A&A 385, 647-670 (2002)
DOI: 10.1051/0004-6361:20020173

Nested-grid calculations of disk-planet interaction

G. D'Angelo^1,2 - T. Henning¹ - W. Kley²

1 - Astrophysikalisches Institut und Universitäts-Sternwarte, Schillergäßchen 2-3, 07745 Jena, Germany
2 - Computational Physics, Auf der Morgenstelle 10, 72076 Tübingen, Germany

Received 16 May 2001 / Accepted 24 January 2002

Abstract
We study the evolution of embedded protoplanets in a protostellar disk using very high resolution nested-grid computations. This method allows us to perform global simulations of planets orbiting in disks and, at the same time, to resolve in detail the dynamics of the flow inside the Roche lobe of the planet. The primary interest of this work lies in the analysis of the gravitational torque balance acting on the planet. For this purpose we study planets of different masses, ranging from one Earth-mass up to one Jupiter-mass, assuming typical parameters of the protostellar disk. The high resolution supplied by the nested-grid technique permits an evaluation of the torques, resulting from short and very short range disk-planet interactions, more reliable than the one previously estimated with the aid of numerical methods. Likewise, the mass flow onto the planet is computed in a more accurate fashion. The obtained migration time scales are in the range from few times 10⁴ years, for intermediate mass planets, to 10⁶ years, for very low and high mass planets. These are longer than earlier assessments due to the action of circumplanetary material. Typical growth time scales depend strongly on the planetary mass. Below 64 Earth-masses, we find this time scale to increase as the 2/3-power of the planet's mass; otherwise it rises as the 4/3-power. In the case of Jupiter-size planets, the growth time scale is several times ten thousand years.

Key words: accretion, accretion disks - hydrodynamics - methods: numerical - stars: planetary systems

1 Introduction

During the past five years radial velocity studies have allowed the detection of planetary companions around other main-sequence stars. Until now about sixty so-called "extrasolar planets'' have been discovered, that orbit their stars within a distance of a few AU. A recent catalog of extrasolar planets, including their orbital characteristics, is provided by Butler et al. (2001) and up-to-date versions can be found at http://www.obspm.fr/encycl/encycl.html and http://exoplanets.org/, maintained by Jean Schneider and the Department of Astronomy at UC Berkeley, respectively.

In contrast to the solar system, these new planets display quite different orbital properties that challenge the accepted formation scenario for solar planets. The major differences are their high minimum masses (up to 17 Jupiter-masses), their proximity to the central star (a fraction of the Sun-Mercury distance) and their high eccentricities (up to 0.7).

One of the main problems to deal with is the very close distance of massive planets to their parent star. The formation of Jupiter-type planets at these locations is, on theoretical grounds, very unlikely. First of all, from purely geometrical arguments, the matter reservoir of the surrounding disk is too small, so that a planet could hardly accrete its mass. Second, the temperatures within the disk are too high for a rocky core to condense easily.

For these reasons it is generally believed that planets have formed from disk material further out, at distances of several AU from the star, and have then migrated to their present positions. This radial motion of the planet through the disk is primarily caused by gravitational torques acting on the planet. The presence of the planet in the disk disturbs the disk gravitationally, creating spiral density wave perturbations, which emanate from the planet through the disk. Hence, the disk is no longer axisymmetric which results in a net torque on the planet. The sign and magnitude of the vertical component of the torque determines the direction and efficiency of the radial migration.

While initial fully non-linear hydrodynamical numerical computations of embedded planets assumed a fixed circular orbit of the planet (Kley 1999; Bryden et al. 1999; Lubow et al. 1999), more recent simulations took into account the back reaction of the disk and allowed for a change in the parameters of the planetary orbit (Kley 2000; Nelson et al. 2000). For a Jupiter-mass planet and typical parameter values for the disk, the obtained orbital decay time is about 10⁵ years, which agrees reasonably well with previous estimates based on analytic linear theories (Goldreich & Tremaine 1980; Ward 1997).

The majority of the computations, performed so far, have used a single grid which resolves the Roche lobe of a Jupiter-mass planet only with very few grid cells. Recently, Ciecielag et al. (2000a, 2000b) used an Adaptive Mesh Refinement method to resolve the immediate surroundings of the planet, but they didn't give any estimate of the mass accretion rate and magnitude of the gravitational torque. On the other hand, Armitage (2001) reduced the overall simulated region achieving a better resolution. However, also in this case the Roche lobe is only resolved by a few grid cells because of the low mass of the investigated planet.

In this paper we aim to discuss the structure and dynamics of the gas flow in the close vicinity of the planet, while performing global disk simulations. In order to obtain the necessary high spatial and temporal resolution, we use a nested-grid formalism which allows an accurate computation of the mass flow onto the planet and the acting torques.

In the next section we lay out the physical model followed by a description of the numerical method (Sect. 3). We describe the setup of the various numerical models in Sect. 4. The main results are presented in Sect. 5 and our conclusions are given in Sect. 6.

   
2 Physical model

For the purpose of this study, we assume that the opening angle of the protostellar accretion disk is very small. We describe the disk structure by means of a two-dimensional, infinitesimally thin model using vertically averaged quantities, such as the surface mass density

$$\Sigma = \int^\infty_{-\infty} \rho {\rm d}z,$$

where $\rho$ is the regular density. We work in a cylindrical coordinate system $(r, \varphi, z)$ whose origin is fixed at the center of mass of the star and the planet, and where the plane of the disk coincides with the z=0 plane.

The gas in the disk is non-self-gravitating and is orbiting a protostar having a mass $\mbox{$M_\bigstar$ }= 1$ $M_{\odot}$. The total mass of the disk $M_{\rm D}$, within the simulated region, which extends from 2.08 to 13 AU, is $3.5 \times 10^{-3}$ $M_{\odot}$. Embedded in this disk there is a massive protoplanet with a mass $M_{\rm p}$, which ranges from one Earth-mass ( $M_{\oplus}$) to one Jupiter-mass ( $M_{\jupiter}$), depending on the considered model. The planet is assumed to be on a fixed circular orbit throughout the evolution. We employ a rotating coordinate system, corotating with the planet, whose azimuthal position is kept constant at $\varphi_{\rm p}=\pi$. The angular velocity $\Omega$ of the rotating frame is then given by

$$\Omega = \Omega_{\rm p} = \sqrt{\frac{G\,(\mbox{$M_\bigstar$ }+ \mbox{$M_{\rm p}$ })}{a^3}},$$ (1)

where G is the gravitational constant and a is the semi-major axis of the planet's orbit.

The evolution of the disk is given by the two-dimensional ( $r, \varphi$) continuity equation for $\Sigma$ and the Navier-Stokes equations for each of the two components of the velocity field $\vec{u}\equiv (u_r, u_\varphi)$. Thus, the set of equations read

$$\frac{\partial \Sigma}{\partial t} + \nabla \cdot (\Sigma\,\vec{u}) = 0,$$ (2)

$$\frac{\partial (\Sigma\,u_r)}{\partial t} + \nabla \cdot (\Si... ...artial r} - \Sigma\, \frac{\partial \Phi}{\partial r} + f_{r},$$ (3)

$$\frac{\partial [\Sigma\, r^2\, (\omega + \Omega)]}{\partial t... ...Sigma\,\frac{\partial \Phi}{\partial \varphi} + f_{\varphi}.$$ (4)

Here $\omega=u_\varphi/\,r$ is the angular velocity and P is the vertically integrated (two-dimensional) pressure. The gravitational potential $\Phi$, generated by the protostar and the planet, is given by

$$\Phi = \Phi_{\bigstar} + \Phi_{\rm p} = - \frac{G\,\mbox{$M_... ...\,\mbox{$M_{\rm p}$ }}{\vert \vec{r} - \vec{r}_{\rm p} \vert},$$ (5)

where $\vec{r}_{\bigstar}$ and $\vec{r}_{\rm p}$ are the radius vectors to the star and the planet, respectively. The effects of viscosity are contained in the terms f_r and $f_{\varphi}$ which give the viscous force per unit area acting in the radial and azimuthal ( $f_{\varphi}/r$) direction:

$$f_{r} = \frac{1}{r}\,\frac{\partial (r\,S_{rr})}{\partial r} ... ...{r\varphi}}{\partial \varphi} - \frac{S_{\varphi\varphi}}{r},$$

$$f_{\varphi} = \frac{1}{r}\,\frac{\partial (r^2\,S_{r\varphi})... ...} + \frac{\partial S_{\varphi\varphi}}{\partial \varphi}\cdot$$

Since we assume a zero bulk viscosity $\zeta$, a constant kinematic viscosity $\nu$ and we do not include any artificial viscosity in our numerical models, the relevant non-zero components of the three-dimensional viscous stress tensor $\tens{S}$ are

$$S_{rr} = 2\,\nu\,\Sigma\,\left(\frac{\partial u_r}{\partial r} - \frac{1}{3}\,\nabla \cdot \vec{u}\right),$$ (6)

$$S_{\varphi\varphi}= 2\,\nu\,\Sigma\,\left(\frac{\partial \ome... ... + \frac{u_r}{r} - \frac{1}{3}\,\nabla \cdot \vec{u}\right),$$ (7)

$$S_{r\varphi} = \nu\,\Sigma\,\left(\frac{1}{r}\,\frac{\partial... ...ial \varphi} + r\,\frac{\partial \omega}{\partial r}\right),$$ (8)

where the divergence of the velocity field can be written as

$$\nabla \cdot \vec{u} = \frac{1}{r}\,\frac{\partial (r\,u_r)}{\partial r} + \frac{\partial \omega}{\partial \varphi}\cdot$$

A more general form of the relations (6), (7) and (8), within the two-dimensional cylindrical approximation, is given for example in Kley (1999).

In the set of equations above we have omitted the energy equation because in this study we will be concerned only with a relatively simple equation of state which does not require the solution of an energy equation. We shall use an isothermal equation of state where the surface pressure P is related to the density $\Sigma$ through

$$P = c_{\rm s}^2 \, \Sigma.$$ (9)

The local isothermal sound speed $c_{\rm s}$ is given here by

$$c_{\rm s} = \frac{H}{r} \, v_{\rm Kep},$$ (10)

where $v_{\rm Kep} = \sqrt{G\,\mbox{$M_\bigstar$ }/r}$ denotes the Keplerian orbital velocity of the unperturbed disk. Equation (10) follows from vertical hydrostatic equilibrium. The ratio h of the vertical height H to the radial distance r is taken as a fixed input parameter. With this choice, the temperature T is proportional to 1/r. Here we use a standard value

   $$h = \frac{H}{r} = 0.05 ,$$

which is typical for protostellar accretion disks having a mass inflow rate of $\dot{M} \approx 10^{-7}$  $\mbox{$M_{\odot}$ }\,{\rm yr}^{\rm -1}$. With this value of h, our kinematic viscosity coefficient is equivalent to $\alpha =4\times 10^{-3}$ at the radial position of the planet ( $\nu=\alpha\,c_{\rm s}\,H$).

Since the mass of the planet is very small in comparison to the mass of the star, because we always use here a ratio $q = \mbox{$M_{\rm p}$ }/\,\mbox{$M_\bigstar$ }$ smaller than 10^-3, the center of mass is located very close to the position of the star. In the following we will often identify the radial distance from the origin of the coordinate system with the distance from the central star.

   
3 Numerical method

In order to study the planet-disk interaction, we utilize a finite difference method to solve the hydrodynamic equations outlined in the previous section. As we intend to achieve a very high resolution around the planet, we use a nested-grid technique. Both requirements are provided by an early FORTRAN-Version of NIRVANA (Ziegler 1997; Ziegler 1998), which is a 3D, nested-grid, MHD code, based on a covariant Eulerian formalism. The relevant equations are solved on a mesh structure having a constant spacing in each direction. For the present purposes, the code is used in a pure hydrodynamic mode, adopting a cylindrical reference frame where the z-dimension is switched off.

NIRVANA uses a spatially, second-order accurate, explicit method. The advection is computed by means of the second-order monotonic transport algorithm, introduced by van Leer (1977), which guarantees global conservation of mass and angular momentum. It is first-order accurate in time. The viscosity part was added and is treated explicitly. The stress tensor has been implemented for a Newtonian fluid according to the Stokes hypothesis, i.e. with a bulk viscosity $\zeta=0$. This code has been already employed in similar computations, either in 2D (Nelson et al. 2000) and 3D (Kley et al. 2001), but always in a single-grid mode.

   
3.1 Nested-grid technique

This technique is particularly useful when very high local resolution is required at specific and predefined points of the computational domain. In our situation, this allows us to simulate both the overall behavior of the disk and the immediate surroundings of the planet. Since this kind of numerical approach is quite new for the calculation of disk-planet interaction, we describe the method to some extent, but only referring to our particular and specific case.

A similar numerical scheme has been adopted, for astrophysical simulations, by a number of authors. Ruffert (1992) used this approach to investigate the collision between a white dwarf and a main sequence star. In his paper the numerical method is explained in detail. Yorke et al. (1993) and Burkert & Bodenheimer (1993) simulated the collapse of a protostellar cloud. An application to flux-limited radiation hydrodynamics can be found in Yorke & Kaisig (1995).

The method relies on the basic idea that, whenever a greater resolution is needed in a designated region, a finer subgrid is located inside the main grid (the one covering the whole computational domain). If the resolution is not high enough yet, another subgrid may be placed on the underlying one. Since any subgrid can host a finer subgrid structure, a grid hierarchy is generated, also called "system of nested grids''. In principle there is no limit to the degree of nesting. A three-level hierarchy is shown in Fig. 1.

The necessary equations are then integrated, independently, on every grid level. However, two neighbor subgrids must exchange the necessary information whenever the integration proceeds from one grid level to another. Restrictions are imposed on the time step only because, for stability reasons, the Courant-Friedrichs-Lewy (CFL) condition must be fulfilled during each integration, on each level.

	Figure 1: Face-up projection of a three-level grid system in cylindrical coordinates. On the finest subgrid (l=3) the linear spatial resolution is four times as large as it is on the main grid (l=1).
Open with DEXTER

   
3.1.1 Basic integration cycle

In our calculations we use the smallest possible refinement ratio:

$$\frac{\Delta r(l+1)}{\Delta r(l)}= \frac{\Delta \varphi(l+1)}{\Delta \varphi(l)}=\frac{1}{2},$$ (11)

where $\Delta r(l)$ and $\Delta \varphi (l)$ represent the mesh discretization, along each direction, on the grid level l (here l=1 identifies the main grid). In order to analyze a complete integration cycle, let's suppose we have a three-level hierarchy at an evolutionary time t. The cycle will be completed when on each grid the system has evolved for the same time:

1.

We always start the integration from the finest level. During the first step of the cycle, this grid is evolved for a time interval

$$\Delta t_1(3) = \min\left[\Delta t^{\rm CFL}_1(3), \frac{1}{... ...^{\rm CFL}_1(2), \frac{1}{4}\,\Delta t^{\rm CFL}_1(1)\right],$$

where $\Delta t^{\rm CFL}_1(l)$ represents the time step resulting from the CFL criterion applied to the level l, after its latest integration. Thus, $\Delta t_1(3)$ accounts for the CFL stability criterion on the whole set of grids.

2.

The third grid has to be integrated once more because of Eq. (11). Since after the first step this was the only level to evolve, we only have to check the new CFL time step for this grid, $\Delta t^{\rm CFL}_2(3)$. Then it can move further in time for an interval

$$\Delta t_2(3) = \min\left[\Delta t^{\rm CFL}_2(3), \frac{1}{... ...^{\rm CFL}_1(2), \frac{1}{4}\,\Delta t^{\rm CFL}_1(1)\right].$$

3.

Now level 2 can be integrated for a time

$$\Delta t_3(2)=\Delta t_1(3) + \Delta t_2(3) \leq \Delta t^{\rm CFL}_1(2),$$

so that numerical stability is automatically assured. At this point of the cycle the first information exchange takes place: the solution on the level 2 is corrected via the more accurate solution of the level 3; the boundary values of the level 3 are updated by using the solution of the level 2, which covers a larger domain. These fundamental interactions will be described later.

We have just seen that a level l+1 has to be visited two times as often as the level l. Then the next two cycle steps will be similar to the first two, provided that the appropriate CFL time steps are employed to compute $\Delta t_4(3)$ and $\Delta t_5(3)$. During the sixth step, the level 2 evolves for $\Delta t_6(2)=\Delta t_4(3) + \Delta t_5(3)$. Eventually, during the seventh step, the integration of the level 1 is performed using a time step $\Delta t_7(1)=\Delta t_3(2) + \Delta t_6(2)$ which is, by construction, smaller than $\Delta t^{\rm CFL}_1(1)$. The cycle is now complete and each grid level has evolved for the same amount of time $\Delta t_7(1)$. The entire cycle sequence is schematically sketched in Fig. 2.

In general, within this kind of cycle, a level l is integrated 2^l-1 times.

	Figure 2: Scheme of a complete integration cycle for a three-level grid system. Arrows indicate the direction of information transfer when integration proceeds from a level to the next lower one. Straight arrows stand for the solution updating process on the levels 2 and 1. Bow-arrows indicate the data transfer for setting the boundary quantities on the levels 3 and 2.
Open with DEXTER

3.1.2 Downward information transfer

We already mentioned that after the third step of the iteration cycle, grids 3 and 2 have to exchange some information. In general, this exchange must occur every time the grid level l evolves to the same time as the level l+1. Because of the higher resolution, we assume the solution of the level l+1to be more accurate than that of the level l. Therefore, the fine-grid solution replaces the coarse one on the common computational domain. Whatever the level in the hierarchy is, the frame formed by the first and the last two grid cells are ghost cells (see Fig. 3). This indicates that they contain the boundary values necessary to perform the algorithm integration. Ghost cells of level l+1 do not contribute to the updating process of the solution of level l.

Figure 3: Interface between a subgrid and its host. The light-colored zone marks those cells containing boundary values needed for the subgrid integration, the "ghost cells''. The darker region refers to the so-called "active zone'', where values are effectively computed on the subgrid. The thick line, which separates the previous regions, encloses the grid cells on the coarse grid whose content is replaced by the more accurate one coming from the subgrid.

Open with DEXTER

The replacement procedure is straightforward: a surface weighted average, using the nearest fine values, substitutes the coarse quantity. For example, referring to Fig. 4, the averaged coarse density ( $\Sigma ^{\rm C}$) is

$$\Sigma^{\rm C} =\frac{\sum_i \Sigma^{\rm F}_i\,A_i}{\sum_i A_... ...gma^{\rm F}_2 + \Sigma^{\rm F}_3)\,A_2}{2\,(A_1 + A_2)}\cdot$$ (12)

In our case, since the advected quantities are the linear radial momentum density $\Sigma\,u_r$, and the angular momentum density $\Sigma\,r^2\,\omega$, these are the interpolated quantities, along with $\Sigma$. Since velocities are centered at the sides of a cell (see Figs. 5 and 6), this average is a little more complex than the previous one and requires six terms. Indicating with $u^{\rm C}$ the coarse value of the linear momentum density to be interpolated and with $u^{\rm F}_i$ the surrounding fine grid values, we have:

$$u^{\rm C} =\frac{(u^{\rm F}_1 + u^{\rm F}_6)\,A_1 + (u^{\rm ... ...u^{\rm F}_3 + u^{\rm F}_4)\,A_3}% {2\,(A_1 + A_2 + A_3)}\cdot$$ (13)

In a successive step, velocities are retrieved. No other quantity needs to be interpolated because the potential is fixed and no energy equation is solved.

Figure 4: Surface weighted average of the surface density. A coarse cell is shown along with the four fine cells it comprises. As a scalar, the surface density is cell-centered within the appropriate refinement (dots). $\Sigma ^{\rm C}$ represents the new, interpolated, value of the coarse cell; $\Sigma ^{\rm F}_i$ are the fine interpolating quantities. Because of the fixed value of $\Delta \varphi (l)$, A1=A4 and A2=A3.

Open with DEXTER

	Figure 5: Surface weighted average of the radial momentum density $u=\Sigma \,u_r$. Two coarse cells are drawn (thick lines); Fine cells are delimited by thin lines. Because of the staggered structure of the grid, vectors are face-centered within the cell.
Open with DEXTER

In order to guarantee global mass and momentum conservation for the whole hierarchy, we have to make sure that the momentum flux components $[\Sigma\,u_r/(r\,\Delta \varphi\,\Delta t),\: \Sigma\,r^2\,\omega/(\Delta r\,\Delta t)]$, across the border between a subgrid and its host grid (indicated in Fig. 3), are equal in both level solutions whenever the grid has evolved for the same time as the subgrid. However, each grid evolves independently and for a time interval different from that of the lower one. Thus, even after the solution updating process described above, the amount of momentum flowed across the borders might not coincide in the respective solutions. To remove this possible discrepancy, at the coarse-fine grid border, these quantities are taken from the fine grid integration. In Fig. 6 the situation for the azimuthal momentum flux is depicted. Two fine cells participate in this process. Referring to the integration cycle traced in Sect. 3.1.1, f^j_k(l) represents the value of the quantity $\Sigma\,r^2\,\omega/(\Delta r\,\Delta t)$, at the grid-grid interface location, as computed during the kth cycle step on level l. An additional index (j=1,2) is needed to identify the radial position of the two fine cells involved (for example, on level 3), but it does not concern the coarse grid quantity to be replaced (on level 2).

Suppose we are at the end of the third cycle, when the first interaction, between levels 2 and 3, occurs (first straight arrow in Fig. 2). Because of the refinement ratio established by Eq. (11), quantity f₃(2) will be reset as:

$$f_3(2) = \frac{1}{2}\,\frac{\Delta t_1(3)\,\sum_jf^j_1(3) + \Delta t_2(3)\,\sum_jf^j_2(3)}{\Delta t_3(2)}\cdot$$ (14)

This correction is accounted for directly while performing the advection of radial and angular momenta.

Figure 6: Momentum flux correction scheme. Momentum flux components are centered as velocity field components, whose locations are shown in the upper left corner. $\tilde{f}(2)\,\Delta r(2)= f_3(2)\,\Delta t_3(2)\,\Delta r(2)$ represents the amount of angular momentum, flowed across the coarse cell border, during the third cycle step (see Sect. 3.1.1). $\tilde{f}^1(3)\,\Delta r(3) = f^1_1(3)\,\Delta t_1(3)\,\Delta r(3) + f^1_2(3)\,\Delta t_2(3)\,\Delta r(3)$ represents the same quantity transported, during the first two cycle steps, across the j=1 fine cell border. The coarse quantity is replaced by $\left[\tilde{f}^1(3)+ \tilde{f}^2(3)\right]\,\Delta r(3)$.

Open with DEXTER

3.1.3 Upward information transfer

The boundary conditions on the main grid are usually imposed depending on the physics and geometry of the problem: symmetry, periodicity etc. In the case of a subgrid, boundary values must be attached, in some way, to the values of the underlying grid. This point turned out to be extremely delicate for our calculations. In NIRVANA, the boundary values of a certain level l were just set by means of a linear interpolation of the quantities on the lower level l-1. Because of the strong variations in density and velocity, due to the formation of shock fronts in our simulations, this method fails and produces numerical inconsistencies.

Therefore, we raised the order of the interpolation. However, this introduces another potential trap. In fact, a high-order interpolation (higher than the first order) is not monotonic and can produce a new minimum. This is not acceptable since, for example, the density is a non-negative quantity. Then the interpolating function should be monotonised. In order to handle this problem we have used the same approach as described in Ruffert (1992), that is by employing the monotonised harmonic mean (van Leer 1977).

Then, if we have a function sampled at $x-\Delta x$, x and $x+\Delta x$, with values $g_{\rm L}$, g and $g_{\rm R}$, respectively (as shown in Fig. 7), the averaged value at $x+\varepsilon$ is

$$g_\varepsilon = g + \frac{2\,\varepsilon}{\Delta x}\, \max\l... ...(g-g_{\rm L})\,(g_{\rm R}-g)} {g_{\rm R}-g_{\rm L}},0\right],$$ (15)

provided that $-\Delta x/2 \leq \varepsilon \leq \Delta x/2$.

If we adopt this kind of average on a 2D-mesh, each interpolation generally involves $3 \times 3$ coarse quantities. It proceeds by averaging the selected coarse values along a certain direction, three at a time. This results in three new quantities. A further harmonic average of these, along the other direction, generates the subgrid boundary value at the correct position. Figure 8 shows how the procedure works in the case of the surface density.

As density is cell-centered, $\varepsilon _r$ and $\varepsilon_\varphi$ are always one fourth of the coarse grid linear size. Since the averaging process is performed for each direction separately, it is not affected by the metric of the mesh, i.e. it is performed as if a Cartesian grid were used.

For the velocity components we have to distinguish two different cases (see Fig. 9): whether the boundary value lies on a coarse cell border or whether it does not. In the first case, either $\varepsilon _r$ or $\varepsilon_\varphi$ is zero. Then only three coarse values participate in the average, along either the azimuthal or the radial direction. In the second case the interpolation proceeds exactly as explained for the surface density.

Figure 7: Behavior of an harmonic mean against a geometric mean. Grey filled circles indicate the harmonic values assumed by $g_\varepsilon$. Open circles indicate the corresponding values when the geometric mean is applied. The plot intends to simulate a shock front. If we keep the left value $g_{\rm L}$ unchanged and make the right one $g_{\rm R}$ four times as large, the geometric mean value will increase proportionally. In turn, the harmonic mean will change so slightly that its variation is not visible on this plot.

Open with DEXTER

Figure 8: Harmonic average: surface density boundary values. The light-colored zone indicates subgrid ghost cells. Darker region belongs to the active subgrid zone. Nine coarse values (big black circles) are engaged in the average. The value to be interpolated is shown as a small black circle. During a first step, Eq. (15) is applied three times along the $\varphi$-direction. Three new values are generated (grey small circles), having the correct, final azimuthal coordinate. These are averaged to obtain the final value at the correct radial location. Small open circles refer to the other boundary quantities whose value depends on the same coarse quantities.

Open with DEXTER

   
4 General model design

The main goal of this study is the investigation of the characteristic features of the hydrodynamic flow within the Roche lobe of the planet. This means that we must be able to resolve a characteristic length, the Hill radius:

$$\mbox{$R_{\rm H}$ }=a\,\left(\frac{\mbox{$M_{\rm p}$ }}{3\,\mbox{$M_\bigstar$ }}\right)^{\frac{1}{3}}.$$ (16)

Moreover, we intend to do that for a variety of planet-to-star mass ratios. In order to reach such resolutions, we build a series of grid systems having the planet located approximately in the middle, at each level of the hierarchy. Smaller planets require higher degrees of nesting.

From now on, we refer to non-dimensional units. All the lengths are expressed in units of the semi-major axis a. This is constant because the planet moves on a fixed circular orbit, with radius

$$r_{\rm p}=\frac{a}{(1 + q)}\cdot$$ (17)

Masses are in units of the central stellar mass and time is given in units of the planet's orbital period. However, in order to convert them into conventional physical units, we assume that a=5.2 AU and, as already mentioned, $\mbox{$M_\bigstar$ }= 1$ $M_{\odot}$. This implies that one planet orbit takes 11.8 years.

The whole azimuthal range of the disk is taken into account by considering a computational domain represented by $2\,\pi \times [r_{\rm in}, r_{\rm out}]$, where $r_{\rm in}=0.4$ and $r_{\rm out}=2.5$. This is covered by a $142 \times 422$ mesh (main grid), allowing a resolution such that $\Delta r(1) = \Delta \varphi(1) = 0.015$ and constraining the resolution on each other grid level, according to Eq. (11). The size of any higher grid level is $64 \times 64$.

Figure 9: Harmonic average of radial velocity component. Two cases are shown. In one case the position of the velocity component lies on a coarse grid border (thick lines), requiring only three coarse quantities (one single average). In the other case nine coarse quantities are needed as for the density interpolation. The only difference is that $\varepsilon _r$ is one half of the radial coarse size.

Open with DEXTER

   
4.1 Smoothing of the potential

The perturbing action of the planet is exerted via its gravitational potential $\Phi _{\rm p}$. From a numerical point of view, it is usually smoothed in order to prevent numerical problems near the planet. Thus, we re-write the denominator of $\Phi _{\rm p}$ in Eq. (5) as $\sqrt{\vert\vec{r}-\vec{r}_{\rm p}\vert^2 + \delta^2}$. However, the smoothing length $\delta$ cannot be the same all over the grid system because of the different grid size involved at each level. Then we use the following grid-dependent length:

$$\delta(l) = \min\left[\frac{\mbox{$R_{\rm H}$ }}{5}, \lambda(l)\right],$$ (18)

where

$$\lambda(l) = \sqrt{\Delta r^2(l) + r_{\rm p}^2\,\Delta \varphi^2(l)} \simeq \sqrt{2}\,\Delta r(l).$$ (19)

The value of constant part, in Eq. (18), has become a kind of standard in single grid simulations (Kley 1999; Lubow et al. 1999; Kley 2000). This is always used on the main grid whereas the grid dependent part always prevails on the highest grid level. This choice results in a very deep potential in the immediate vicinity of the planet.

   
4.2 Accretion onto the planet

The presence of the planet affects the nearby disk density also because it can accrete matter. Planet accretion is accounted for by removing some mass from the region defined by $\vert\vec{r}-\vec{r}_{\rm p}\vert\leq \kappa_{\rm ac}$. Since mass is removed from the system after each integration step, the evacuation rate depends on an input parameter $\kappa_{\rm ev}$ as well as on the integration time step. The details of the accretion process are given in Fig. 10. This removal is accomplished only on the highest (finest) hierarchy level and the removed mass is not added to the dynamical mass of the planet, but just monitored. A standard value $\kappa _{\rm ev}=5$ is used, while we set $\kappa_{\rm ac}$ between $8.0\times10^{-2}$ and $9.4\times10^{-2}$ $R_{\rm H}$. These values are such that only few grid cells, on each side of the planet, are involved in the gas accretion, making it a locally confined process. For the Jupiter-mass planet, ${\sim } 12\times 12$ cells participate. The lowest number of employed cells is ${\sim} 8\times8$, which is used for the smallest planet ( $\mbox{$M_{\rm p}$ }=1$  $M_{\oplus}$). The largest is ${\sim} 18 \times 18$ and is adopted for a planet with $\mbox{$M_{\rm p}$ }=0.5$  $M_{\jupiter}$ (160  $M_{\oplus}$). This circumstance is due to the high numerical resolution of the model (see Table 1).

Note that the sphere of influence of the accretion process consists of a region with a typical radius of 10^-1 $R_{\rm H}$, quite smaller than the protoplanetary radius, whereas the planet is usually assumed to fill its Roche lobe during the growth phase (Bodenheimer & Pollack 1986; Tajima & Nakagawa 1997). As the present study does not take the energy equation into account, it precludes a detailed treatment of the internal structure of the protoplanet (see Wuchterl et al. 2000 and references therein; see also Fig. 16 below), hence the inferred accretion rates may still be unreliable.

Figure 10: Planet mass accretion scheme. Wherever a cell center lies in the dark-colored zone, after each integration its density is lowered by an amount $\Delta \Sigma$ such that $\Delta \Sigma /\Sigma = 2\,\pi \, \kappa _{\rm ev}\,\Delta t$, where $\Delta t$ is the integration time step interval. If the cell center falls into the light-colored zone $\Delta \Sigma$ is only one third of that value. For a typical main grid time interval and $\kappa _{\rm ev}=5$, we would get $\Delta \Sigma /\Sigma \approx 1.2\times 10^{-2}$. Anyway, this evacuation process is performed only on the highest grid level for better accuracy.

Open with DEXTER

4.3 Initial and boundary conditions

The initial density distribution is proportional to r^-1/2. However, we superimpose to this an axisymmetric gap around the planet, obtained by an approximate balance of the viscous torque and the gravitational torque due to the planet (P. Artymowicz, private communication). Figure 11 shows the surface density, at t=0, for some selected planet masses. The initial velocity field is that of a Keplerian disk.

As boundary conditions, periodicity is imposed at $\varphi=0$ and $\varphi=2\,\pi$. We allow matter to flow out of the computational domain at the inner radial border ( $r_{\rm in}$) whereas we set reflective boundary conditions at the outer radial border ( $r_{\rm out}$). The angular velocity is set equal to the unperturbed Keplerian value $\Omega_{\rm Kep} = \sqrt{G\,\mbox{$M_\bigstar$ }/\,r^3}$, both at $r_{\rm in}$ and $r_{\rm out}$.

For low mass planets ( $\mbox{$M_{\rm p}$ }\lesssim 10$  $M_{\oplus}$), boundary conditions should not affect much the system evolution because density waves damp before reaching $r=r_{\rm out}$ and are very weak when they reach $r=r_{\rm in}$. For more massive planets, some reflection is seen at the outer radial border. Further, the torque exerted by the planet pushes the inside-orbit material inwards. As a result, because of the open boundary at $r=r_{\rm in}$, the inner disk is partially cleared and consequently positive torques originating there reduce. The higher the planet mass is the stronger these two effects appear.

	Figure 11: Initial surface density distribution for q=10-3, q=10-4, and q=10-5. $\Sigma (r,\varphi )=\Sigma (r)$ at t=0. The gap greatly reduces when q gets small enough.
Open with DEXTER

   
4.4 Model specifications

In this paper, we are mainly interested in investigating how disk-planet interactions change by varying q, the planet-to-star mass ratio. Table 1 summarizes the values of q used in different models along with the number of grid levels employed. For reference, a measure of the linear resolution $\xi$, in units of the Hill radius $R_{\rm H}$, is given for the highest level as well. Some models have different prescriptions than the ones outlined above, as specified in the Table 1.

Few models may deserve some comments. Model ELEN2 and WPRO1 aim at checking whether results from CIRO1 and CIRO2 (respectively) are resolution-dependent. This test is negative, as we show in the next sections. Since the planet position $(r_{\rm p},\varphi_{\rm p})$can fall anywhere within a grid cell according to the value of q in Eq. (17) and to the definition of the grid, some asymmetries could arise. These might have some effects on the finest levels, due to the small value of the smoothing factor $\delta$. In order to achieve a complete symmetry, in the model GINO1 the planet is placed at the corner of a main grid cell (i.e., the intersecting point of four grid cells). This property is such that, on every other gird level, the planet always sits on the cross-point of four grid cells. CIRO3 and its counterpart, GINO1, give almost identical results. Since in the Jupiter-mass case the inner-disk is greatly depleted, model ELEN1 was run to evaluate the influence of its presence on the gravitational torque. For this reason, in such model we prevent matter from draining out of the inner radial border by setting reflective boundary conditions.

 

 

Table 1: Model specific parameters: $q = \mbox{$M_{\rm p}$ }/\,\mbox{$M_\bigstar$ }$; ng is the number of grid levels; $\xi= \Delta r(ng)/\mbox{$R_{\rm H}$ }$ is the normalized grid resolution on the finest level. Each model is let evolve, at least, till 200 orbits.

Model	q	ng	$\xi$	Notes
CIRO1	$1.0 \times 10^{-3}$	5	$1.3 \times 10^{-2}$
CIRO2	$1.0 \times 10^{-4}$	6	$1.5 \times 10^{-2}$
CIRO3	$1.0 \times 10^{-5}$	7	$1.6 \times 10^{-2}$
PEPP1	$3.0 \times 10^{-6}$	7	$2.3 \times 10^{-2}$
PEPP2	$1.5 \times 10^{-5}$	7	$1.4 \times 10^{-2}$
PEPP3	$3.0 \times 10^{-5}$	7	$1.1 \times 10^{-2}$
PEPP4	$6.0 \times 10^{-5}$	7	$8.6 \times 10^{-3}$
WPRO1	$1.0 \times 10^{-4}$	7	$7.3 \times 10^{-3}$
WPRO2	$1.0 \times 10^{-4}$	6	$1.5 \times 10^{-2}$	(a)
GINO1	$1.0 \times 10^{-5}$	7	$1.6 \times 10^{-2}$	(b)
GINO2	$2.0 \times 10^{-4}$	6	$1.2 \times 10^{-2}$
GINO3	$5.0 \times 10^{-4}$	6	$8.5 \times 10^{-3}$
ELEN1	$1.0 \times 10^{-3}$	5	$1.3 \times 10^{-2}$	(c)
ELEN2	$1.0 \times 10^{-3}$	6	$6.8 \times 10^{-3}$
ELEN3	$5.0 \times 10^{-4}$	6	$8.5 \times 10^{-3}$	(d)
(a) Same as CIRO2, but $\kappa_{\rm ev}=10$.
(b) Same as CIRO3, but the planet is symmetrically
placed with respect to the grid cells.
(c) Same as CIRO1, but reflective boundary conditions
are set at $r=r_{\rm in}$.
(d) Same as GINO3, but $\kappa_{\rm ac}$ is 40% smaller.

   
5 Results

Hereafter we mainly discuss three models, namely CIRO1, CIRO2, and CIRO3. We concentrate on them because they cover a mass range from 1  $M_{\jupiter}$ down to 3.2  $M_{\oplus}$. Nevertheless, whenever required by our discussion, we mention other particular models. Some results, concerning the whole set of models given in Table 1, are presented as well.

   
5.1 Overall flow structure

Large-scale interactions (whose effects extend over $2\,\pi$ in azimuth and cover a large radial extent) of a Jupiter-mass planet with the surrounding environment, have already been treated numerically in a number of papers (Artymowicz 1992; Kley 1999; Bryden et al. 1999; Kley 2001). An example of large-scale features can be seen in Fig. 12, where the comprehensive result of a nested-grid computation is displayed.

Figure 12: Overview of the surface density $\Sigma$ on the six-grid system of model GINO3 ( $\mbox{$M_{\rm p}$ }=0.5$  $M_{\jupiter}$), after 210 orbits. The top-left panel shows the main grid (l=1) which encloses the whole simulated disk region. From left to right and top to bottom, the other panels zoom in onto the planet. The highest grid level lies over a surface, roughly $0.5 \times 0.5$ $\mbox{$R_{\rm H}$ }^2$, which comprises just four grid-cells on the main grid. The over-plotted curve represents the Roche lobe of the restricted three-body problem. Plus signs indicate the positions of L1 and L2 Lagrangian points, respectively on the left and the right of the planet.

Open with DEXTER

Planets with a lower value of q should perturb the disk less and have a weaker large-scale impact on it. For these reasons we discuss only the medium ( $\vert\vec{r}-\vec{r}_{\rm p}\vert \sim \mbox{$R_{\rm H}$ }$) and small ( $\vert\vec{r}-\vec{r}_{\rm p}\vert \ll \mbox{$R_{\rm H}$ }$) scale effects of such interactions.

Nonetheless, it's worthwhile noticing that large-scale structures are clearly visible in our smallest mass models. In CIRO3 ( $\mbox{$M_{\rm p}$ }=3.2$  $M_{\oplus}$), a trailing density wave, emanating from the planet, spirals around the star for about $4\,\pi$, vanishing approximately at r=2. In PEPP1 ( $\mbox{$M_{\rm p}$ }=1$  $M_{\oplus}$), a similar feature spirals for almost $2\,\pi$, disappearing at r=1.5. Although we did not investigate this issue any deeper, it may happen that results from local simulations could be influenced by not taking into account entirely such global features.

In Fig. 13 the surface density is shown along with the velocity field for the three CIRO-models. As a reference, the Roche lobe (of the relative three-body problem) is over-plotted. From the upper row of this figure ( $\mbox{$M_{\rm p}$ }=1$  $M_{\jupiter}$) we can see the patterns of the two main spirals (left panel). They reach the Roche lobe, but do not enter it. In fact, they are replaced by two streams of material which start from two points (located at r=0.94, $\varphi=3.12$and r=1.07, $\varphi=3.07$, respectively), where the flow is nearly at rest with respect to the planet. Each of them enters the Roche lobe, encircles the planet and hits the other one on the opposite side. As a result of the collision, the material is shocked and the flow is redirected towards the planet. Hence, these streams assume the form of two spirals, winding around the planet (right panel) for $2\,\pi$. That such smaller scale spirals are detached from the main ones can be inferred from the direction of the flow. Along the main spirals the material follows the disk rotation around the star, moving away from the planet. Along the small ones the gas orbits the planet. In fact, they represent the outstanding features of a circumplanetary disk. A more detailed description of the flow regions around and inside the Roche lobe, concerning a Jupiter-mass planet, can be found in Lubow et al. (1999).

Figure 13: Surface density and velocity field for three selected models. Top row. CIRO1: grid level l=3 (left panel) and l=4 (right panel) at t=375 orbits. Middle row. CIRO2: grid level l=4 (left panel) and l=5 (right panel) at t=325 orbits. Bottom row. CIRO3: grid level l=5 (left panel) and l=6 (right panel) at t=225 orbits. To avoid too much confusion, only 40% of the velocity field vectors are drawn. The over-plotted curve and plus signs have the same meaning as in Fig. 12. Close to the planet, the velocity field is roughly keplerian (see Sect. 5.5).

Open with DEXTER

The case $\mbox{$M_{\rm p}$ }= 32$  $M_{\oplus}$ (Fig. 13, middle row) has many analogies to the previous one. This planet is able to open a gap in the disk, as a permanent feature. However, it is neither so wide (the base width is 0.15 against 0.4) nor so deep (40% against 0.5% of the maximum surface density) as it is for a Jupiter-mass planet. The overall behavior of the matter entering the Roche lobe is very similar (left panel). The up-stream disk material, relative to the nearest X-point, reaches it, inverts partially the direction of its motion and flows into the Roche lobe. The gas stream penetrating from the left X-point turns about the planet, at $\varphi < \varphi_{\rm p}$, and collides with the stream incoming from the other X-point, generating the upper spiral arm (at $\varphi > \varphi_{\rm p}$). However, here the locations, from which these gas streams depart (r=0.97, $\varphi=3.14$and r=1.03, $\varphi=3.135$, respectively), lie closer to the L1 and L2 points. The circumplanetary spirals are less twisted around the planet than before. They wrap around it for an angle $\pi$ (right panel).

For the less massive planet, $\mbox{$M_{\rm p}$ }=3.2$  $M_{\oplus}$, the situation is somewhat different (bottom row of Fig. 13). In fact, within the Roche lobe, the signs of spiral fronts are very feeble, though some traces can still be seen. They assume the shape of a bar-like structure which extends, for 0.3 $R_{\rm H}$, from side to side of the planet at $\varphi\simeq \varphi_{\rm p}$. As indicated by the velocity field, the circumplanetary disk roughly occupies the entire Roche lobe (right panel).

Taking into account the other models as well, the following scenario can be proposed: the lower the value of q, the shorter and straighter the circumplanetary spirals become. For example, in model PEPP2, they track a tilde-like pattern, extending for a total length of 0.4 $R_{\rm H}$.

5.2 Density in the planet's environment

So far we have described qualitatively the main hydrodynamic structures which are present, near the planet, on the length scale of a Hill radius. We shall now discuss, more quantitatively, two aspects of the surface density, observable on shorter length scales ($\lesssim$0.5 $R_{\rm H}$). Both of them have repercussions on the torque exerted on the planet by the close-by matter, as explained later. First of all we should notice that $\Sigma(r,\varphi)$ is not perfectly symmetric with respect to the planet, in none of our reference models. Such property can be checked by means of an accurate examination of Fig. 14. In this figure the contour lines of the surface density (i.e., lines on which $\Sigma$ is constant) are plotted, for the finest grid levels l=ng. In case of CIRO1 (top panel, Fig. 14), one can see that the right-hand arm ( $\varphi > \varphi_{\rm p}$) twists towards the planet more than the other one does. The curvature ratio of the more external spiral parts (indicated by thick lines in the figure) is 0.9. Furthermore, the right-hand arm lies closer to the planet than the left-hand one, as can be evaluated from the positions of the arc centers (marked with plus signs). Contour line asymmetries are fainter in the case of CIRO2 (middle panel, Fig. 14). However, a quantitative analysis shows that the external ridges of the spirals (thick lines) have slightly different curvatures. Moreover, their centers (plus signs) are not aligned with the planet position (represented by a big dot). In particular, if we consider the contours between 1.03 and 2.02 (arbitrary units), the right-hand ridge is a little straighter than the left-hand one. In fact, the ratio between their curvatures is 0.93. Also in this case, the arc centers indicate that external parts of the right-hand arm are a little nearer to the planet than those of the left-hand arm. Indeed, surface density asymmetry is evident, in the case of CIRO3 (bottom panel, Fig. 14). At $r\simeq 1.001$ and $\varphi < \varphi_{\rm p}$, the density is systematically lower than it is on the opposite side ( $r\simeq 0.999$, $\varphi > \varphi_{\rm p}$). Following the shape of density contours, it is possible to track what remains of the spiral arms. The pre-shock material conveys the more convex form to these lines, at $r\simeq 0.999$ and $r\simeq 1.001$

	Figure 14: Surface density contours on the finest grid level. The evolutionary time is the same as in Fig. 13. Top: CIRO1. Middle: CIRO2. Bottom: CIRO3. See the text for an explanation of the thick bow-lines and the plus signs.
Open with DEXTER

   
5.2.1 The core

Finally, we would like to mention what happens, in our models, at very short distances from the planet. We described, in Sect. 4.2, how material is removed from the neighborhood of the planet. This is usually a small fraction of the available matter, during each integration time step of the main grid. If $\Delta M^{\rm ev}$ is the evacuated amount of mass, then typically

$$\frac{\Delta M_{\rm ev}}{M}\approx \frac{\Delta \Sigma}{\Sigma}\approx 10^{-2}$$

(see the caption of Fig. 10). This choice accounts for the fact that a planet should not be able to accept, very rapidly, all the material the surrounding environment can offer (Wuchterl 1993).

Since not all of the matter is taken away, it should pile up at the location of the planet, eventually forming a very dense core. Indeed, this is what we find, as already visible in Fig. 14, where density contour lines crowd around the planet at $(r_{\rm p},\varphi_{\rm p})$. Figure 15 displays better how this feature looks like for models CIRO1 and CIRO2. In order to make a comparison of the linear extension of such cores with the Hill radius, we introduce the Hill coordinates. These are defined as a Cartesian reference frame with origin on the planet and coordinates normalized to $R_{\rm H}$ (see the caption of Fig. 15 for details).

Figure 15: Surface density plot of the highest grid level. Top: CIRO1 at t=375 orbits. Bottom: CIRO2 at t=325 orbits. The Hill coordinates, $(x_{\rm H}, y_{\rm H})$, are so defined: $x_{\rm H}= [r_{\rm p}\,\cos(\varphi_{\rm p}) - r\,\cos(\varphi)]/\mbox{$R_{\rm H}$ }$, $y_{\rm H}= [r_{\rm p}\,\sin(\varphi_{\rm p}) - r\,\sin(\varphi)]/\mbox{$R_{\rm H}$ }$. The star is situated along the negative $x_{\rm H}$-axis. The azimuth $\varphi$ increases with $y_{\rm H}$.

Open with DEXTER

In these two cases, the core width, $2\,\eta$ (exactly defined below), can be estimated to be approximatively 0.1 $R_{\rm H}$.

One reason for the sharpness of these peaks is the very small length scale we adopt to smooth the potential of the planet (Sect. 4.1). On the finest grid level the smoothing length $\delta$ is equal to $\lambda(ng) \simeq 1.4\,\Delta r(ng)$ (see Eq. (19)). From Table 1 one can deduce that $\delta \sim 10^{-2}$ $R_{\rm H}$ (though it changes a little for the different values of q). Despite that, the core width always equals at least $6\,\Delta r(ng)$. Other two hints suggest that such features are not just a numeric product. According to models CIRO2 and WPRO1 (which differ only in the number of grid levels), the core width is very similar. If we cut the peak at $\Sigma=2$ on the finest level of each model, its extent, at $r=r_{\rm p}$, is 6 and 11 grid cells, respectively. Models CIRO1 and ELEN2 agree in a very similar way. Furthermore, the structure of these peaks does not depend on the exact placement of the planet within a grid cell, as models CIRO3 and GINO1 prove. Since the amount of core material may be related to the accretion process, the parameters $\kappa_{\rm ev}$ and $\kappa_{\rm ac}$ are likely to modify the local density. The core widths (as introduced before) of models CIRO2 and WPRO2 differ by about 20% while the difference is somewhat less in case of GINO3 and ELEN3.

Figure 16: The core width is generally much larger than that expected for an isothermal, hydrostatic, configuration. In fact the centrifugal force, due to axial rotation of the gas around the planet, may play an important role in supporting the structure. This is not the case for the core around an Earth-mass planet (model PEPP1). Thus, the structure equation reads: ${\rm d}\,\Phi _{\rm p} = -c_{\rm s}^2\,{\rm d}\,\Sigma /\Sigma$, where $\Phi _{\rm p}$ is the smoothed potential of the planet. The constancy of the sound speed ( $c_{\rm s} = h\, \sqrt{G\,\mbox{$M_\bigstar$ }/a}$) implies that the gas is isothermal. In the figure, the surface density of PEPP1, at $\varphi =\varphi _{\rm p}$ ( $y_{\rm H}=0$), is represented by "$\times$'' signs. The over-plotted curve is the solution of the previous structure equation. Clearly, the peak density $\Sigma _0$ is somehow connected to the accretion procedure.

Open with DEXTER

However, in order to study the exact properties of such features a more detailed physical treatment is required in proximity of the planet. This may be part of future computations which include a more detailed treatment of the internal constitution of the protoplanet. Anyway, properties such as the local surface density profile and the velocity field indicate that the core structure could be approximated to that of a rotating and isothermal gas, in hydrostatic equilibrium. In model PEPP1 (1  $M_{\oplus}$), the rotation of the core is very slow. Hence, its density profile can be well reproduced by that of an isothermal and hydrostatic gas, as demonstrated in Fig. 16, where we compare the numerical results (crosses) with an analytic solution for a hydrostatic isothermal core (solid line).

As pointed out, the material within the core region is strongly coupled to the planet, due to the small distances. In some way, it may be considered as part of the protoplanet itself, whose structure we may not resolve well enough in the present paper. Whatever its nature, it is very likely that the angular momentum transferred, by the core material to the planet, may influence the planet's spin angular momentum rather than its orbital one. As we are treating the planet as a point mass we cannot estimate its spin. Therefore we decided to exclude this region from the torque computation. To do that, we need a quantitative estimate of the core radius $\eta$, for every model in Table 1. We adopt the following procedure: the average density $\bar{\Sigma}$ is computed over the region $0.2\ \mbox{$R_{\rm H}$ }\leq \vert\vec{r}-\vec{r}_{\rm p}\vert \leq 0.5\ \mbox{$R_{\rm H}$ }$. Then we define the core width ($2\,\eta$) to be that where $\Sigma = 2\,\bar{\Sigma}$. In Fig. 17 (left panel) the dependence of $2\,\eta$versus q is displayed. The ratio $\eta/\mbox{$R_{\rm H}$ }$ decreases for increasing values of q. However, between 3.2 and 32  $M_{\oplus}$, it seems to vary very little. Our measure of the core sizes is performed at a particular time (for CIRO-models, they are indicated in Fig. 13). Anyway, such estimates are not affected much by this choice because the cores reach a steady state, early during the system evolution. As an example, Fig. 17 (right panel) shows how the core mass $M(\eta )$ assumes a static value very soon.

As already mentioned, the spiral features vanish for low values of q, and for Earth-mass planets the core becomes the most prominent feature within the Roche lobe.

	Figure 17: Left panel. Core width, $2\,\eta$, in units of the Hill radius as a function of the mass ratio q. CIRO and PEPP-models are considered along with GINO2 and GINO3. The measure errors are assumed to be equal to $\xi$ (see Table 1). Righ panel. Mass within a distance $\eta$ from the planet versus time. The core mass $M(\eta )$ is normalized to the mass of the planet.
Open with DEXTER

	Figure 18: Radial distribution and two-dimensional contour map of the gravitational torque exerted by the disk on the planet. Left panels. CIRO1: levels l=4 and l=5. Center panels. CIRO2: levels l=5 and l=6. Right panels. CIRO3: levels l=6 and l=7. The two-dimensional torque distribution (bottom panels) is reported for the finest grid level and is normalized to its maximum, absolute, value.
Open with DEXTER

   
5.3 Torque exerted on the planet

Any protoplanet embedded in a protostellar environment suffers gravitational torques, exerted by the surrounding disk material. If the density distribution were symmetric with respect to the planet, or with respect to the line connecting the star with the planet, the resultant total torque would be zero. However, we have just seen that this is not the case, and even around a planet as small as 3.2  $M_{\oplus}$, the matter distribution is not axially symmetric. Thus, we expect a non-zero total torque from the disk, acting on the planet. In response to it, because of the conservation of the orbital angular momentum, the planet has to adjust its semi-major axis, which leads to a migration phenomenon.

In the present computations, we evaluate the torque exerted on the planet in a straightforward way. First the gravitational force acting on the planet $\vec{f}_{\rm g}(i,j)$, due to each grid cell (i,j), is calculated. The torque, with respect to the star, is then

$$\vec{t}(i,j)= \vec{r}_{\rm p} \times \vec{f}_{\rm g}(i,j).$$ (20)

Since we are interested in the z-component of this vector, we perform the scalar product $\hat{\vec{z}}\cdot\vec{t}(i,j) = t_z(i,j)$, where $\hat{\vec{z}}$ is the unit vector of the vertical direction. Finally, by summing over all i and j, we obtain the total disk torque $\mathcal{T}_{\rm D}$, whereas the simple contraction over the azimuthal index j gives the radial distribution of the torque.

The quantity t_z(i,j) is computed on each grid level. Where the computational domain is covered by more than one level, the torque t_z(i,j) on the finest grid is considered for the evaluation of $\mathcal{T}_{\rm D}$. We don't work with time averaged torques because they were checked not to vary considerably with time during the last 100 orbits of the computation. Thereupon we use their values at the end of the simulations.

In order to avoid the region dominated by the core, we exclude a certain area from the computation of $\mathcal{T}_{\rm D}$. Because of the way we defined the core radius $\eta$, some core material may still lie outside $\vert\vec{r}-\vec{r}_{\rm p}\vert=\eta$. Therefore, for safety reasons, we choose not to take into account the planet neighborhood defined by $\vert\vec{r}-\vec{r}_{\rm p}\vert\leq \beta = 2\,\eta$. The only level, upon which this operation is relevant, is the highest since it provides the gravitational torque form the regions closest to the planet. On coarser levels such operation is meaningless since inner torques are taken from elsewhere. But it can be useful to confer a more regular look to the radial torque profile. We generally adopt the value $\beta = 5\,\delta(l)$ for l<ng.

However, as this might be somewhat arbitrary, we discuss in a separate section how the choice of $\beta$, on the finest level, affects the value of the total torque $\mathcal{T}_{\rm D}$.

Because of the global angular momentum transfer, the disk material (orbiting the star), at $r>r_{\rm p}$, exerts a negative torque on the planet whereas the inside-orbit gas tends to increase its angular momentum. This tendency changes, as material closer to the planet is accounted for, and it may reverse eventually, once in the circumplanetary disk (for $\mbox{$M_{\rm p}$ }=1$  $M_{\jupiter}$, such behavior was also found by Lubow et al. 1999). The radial torque distribution, from CIRO-models, is shown in Fig. 18. The two sets of profiles belong to the grid levels l=ng-1(top row) and l=ng (middle row). In the case of CIRO1 (top row, left panel), the sign reversal of the torque is not completed yet. However some negative torques are exerted from regions inside the planet orbit and some positive torques are exerted from the opposite side. On the domain covered by this grid (fourth level), the torque contribution coming from $x_{\rm H} < 0$ is positive while the one coming from $x_{\rm H} > 0$ is negative. The resulting net torque is positive, as the magnitude of the latter contribution is 2.4 times smaller than that of the former. The torque behavior gets more complex if we restrict to a region closer to the planet (middle row, left panel). Though not evident at a first glance, the signs are reversed if compared to the preceding grid level. The torque exerted by the region $x_{\rm H} > 0$ is definitely positive and, in magnitude, almost 30 times as large as that arising from the region $x_{\rm H} < 0$. Thus, this region exerts a strong, positive, net torque. Indeed, the phenomenon of the torque sign reversal is clear in the case of CIRO2 and CIRO3 (Fig. 18, center and right panels). For both grid levels, inside-orbit material lowers the angular momentum of the planet while outside-orbit material acts in the opposite direction. On the finest level of CIRO2, the ratio of the negative to the positive torque contribution is just 0.96 (in absolute value), whereas it is 0.3 for CIRO3.

In order to check which is the overwhelming contribution, between positive and negative torques, on the various grid levels, we can define the partial torque $\mathcal{T}(l)$. This represents the total torque computed over the entire domain but the part covered by the grid level l+1 and such that $\mathcal{T}(ng)=\mathcal{T}_{\rm D}$. Figure 19 illustrates the sign of the partial torque and its relative strength, on each grid level, for CIRO-models. We can see that the total torque is negative in all of the three models. In CIRO1 and CIRO2 all levels, but the highest, contribute to lower the planet angular momentum. On the contrary, the matter inside the finest level raises the overall torque (of a fair amount in case of CIRO1). In CIRO3, a positive torque is exerted by the material outside a region, around the planet, with a linear extension ${\sim} 2\,h\,r_{\rm p}$. Instead, levels 3, 4, 5, and 6 provide negative torques, which are then weakened considerably by the positive torque coming from the finest level.

	Figure 19: Partial torque $\mathcal{T}(l)$, normalized to $\vert\mathcal{T}_{\rm D}\vert$, plotted against the grid level l for: CIRO1, CIRO2, and CIRO3. $\mathcal{T}(ng)/ \vert\mathcal{T}_{\rm D}\vert = -1$ because, by definition, $\mathcal{T}(ng)=\mathcal{T}_{\rm D}$ and $\mathcal{T}_{\rm D}$ is negative for these models.
Open with DEXTER

An overview of the torques exerted by different portions of the disk, for all the relevant models, is given in Table 2. First of all we notice that, within 0.5 $R_{\rm H}$, the phenomenon of the sign reversal is observed in all of the models, but PEPP1 and GINO2. Anyway, in the latter case, the torque inside this region is positive.

In CIRO-models and in model GINO3, the reversal of the sign is such to produce a positive torque over the domain $\vert\vec{r}-\vec{r}_{\rm p}\vert < 0.5$ $R_{\rm H}$. This is because the surface density of the leading material is slightly higher than that of the trailing matter (Sect. 5.1 and Fig. 12). In contrast, the torque is negative in case of PEPP2, PEPP3, and PEPP4. We note, however, that the positive torques exerted by the outside-orbit material, within this region, strongly attenuate the magnitude of the negative net torque exerted by the rest of the disk.

Clearly, since neighboring material may tend to reduce the magnitude of negative torques acting on the planet, it could either slow down its inward migration or reverse the direction of its motion.

As anticipated, the sign reversal of the radial torque distribution is not observed in model PEPP1 ( $\mbox{$M_{\rm p}$ }=1$  $M_{\oplus}$). Due to the very low mass of the perturber, the only structure present inside the Roche lobe is the density core. Whatever level is considered, the inside-orbit gas always provides positive torques while negative ones come from the outside (see also Table 2). Negative torques are somewhat stronger, on any grid level. Almost the 50% of the total torque is generated between ${\sim} 0.5\,h\,r_{\rm p}$ and ${\sim} h\,r_{\rm p}$, at the starting positions the disk spirals.

  

Table 2: Gravitational torque exerted on the planet arising from different disk regions, for the relevant models of Table 1. The entire domain is divided into three region: outside the Hill circle; inside the circle of radius $\mbox{$R_{\rm H}$ }/2$ and the zone in between. Then they are divided further in order to distinguish between inside-orbit (in) and outside-orbit (out) contribution. The torque is computed over each of these regions, by employing $\beta =2\,\eta$ and according to the value of the normalized distance from the planet $d=\vert\vec{r}-\vec{r}_{\rm p}\vert/\mbox{$R_{\rm H}$ }$. Torques are normalized to $\vert\mathcal{T}_{\rm D}\vert$.

$$\begin{tabular}{lrrrrrr} \hline \hline &% \multicolumn{2}{... ...$-2.06$ & $-0.34$ & $1.90$ & $-4.08$ \\ \hline \end{tabular}$$

   
5.3.1 2D-torque distribution

The detailed balance of the torques, arising from very close material, depends on the medium and small-scale density structures around the planet, such as the shape of spirals. Therefore, referring to what was stated in the previous section, we can deduce that the radial torque asymmetry is a direct consequence of the asymmetric distribution of the gas with respect to the planet.

Thus, a more comprehensive description of the torque behavior requires its full two-dimensional distribution. In the bottom panels of Fig. 18, the contour lines of the two-dimensional torque, t_z, are shown for each reference model. The interesting point to notice here is that the largest magnitude torques arise from the corotation locations, i.e. where $r=r_{\rm p}$( $x_{\rm H} \simeq 0$). Here in fact, $\vec{f}_{\rm g}(i,j)$ is perpendicular to $\vec{r}_{\rm p}$ and the cross-product in Eq. (20) achieves its maximum (minimum). The material leading the planet (at $\varphi > \varphi_{\rm p}$ or $y_{\rm H}>0$), pulls it ahead and makes it gain angular momentum. The trailing material brakes the planet making it lose angular momentum.

Let's consider two fluid elements at $(0,\pm\,y_{\rm H})$ and write their mass density as $\Sigma^\pm$. Then we can write $\vert t_z^\pm\vert \propto r_{\rm p}\,\Sigma^\pm/y_{\rm H}^2$, which yields:

$$\vert t_z^+\vert - \vert t_z^-\vert \propto r_{\rm p}\,\frac{\Delta \Sigma}{y_{\rm H}^2},$$ (21)

where $\Delta \Sigma = \Sigma^+ - \Sigma^-$. Any mismatch of the surface density, $\Delta \Sigma$, causes a torque mismatch amplified by an amount equal to $y_{\rm H}^{-2}$. It's worth noticing that, on larger distances $\vert\vec{r}-\vec{r}_{\rm p}\vert \sim r_{\rm p}$, the torque mismatch is amplified less, in fact:

$$\vert t_z^+\vert - \vert t_z^-\vert \propto \frac{\Delta \Sigma}{r_{\rm p}}\cdot$$ (22)

This is the reason why surface density asymmetries near the planet have a very strong impact on $\mathcal{T}_{\rm D}$, and they can easily prevail against the more distant ones.

The region responsible for the maxima and minima of the radial torque distributions in the middle row of Fig. 18 can be identified by means of the 2D-torque maps. For example, in the case of CIRO1, we see that the maximum at $x_{\rm H}=-0.15$ and the minimum at $x_{\rm H}\simeq 0.15$ are produced at $y_{\rm H}\simeq 0.1$ and $y_{\rm H}\simeq -0.1$, respectively. In the other two cases, radial distribution extremes rise from regions where the torque function t_z is steeper than it is on the opposite side of the planet.

   
5.4 Planet migration

If we consider a planet moving on a circular orbit, we find that the rate of change of its semi-major axis a, caused by an external torque $\mathcal{T}_{\rm D}$, is

$$\frac{{\rm d}a}{{\rm d}t} = \frac{2\,\mathcal{T}_{\rm D}}{\mbox{$M_{\rm p}$ }\,a\,\Omega_{\rm p}}\cdot$$ (23)

Analytical estimates of $\mathcal{T}_{\rm D}$ show that two limiting cases exist, depending on whether the planet is massive enough to generate a gap or not. In the first case of small planetary masses, we have the so-called type I migration. Ward (1997) derives the following expression:

$$\left(\frac{{\rm d}a}{{\rm d}t}\right)_{\rm I} \simeq -\frac{... ... p}\left(\frac{\pi\,a^2\,\Sigma}{\mbox{$M_\bigstar$ }}\right).$$ (24)

The direction of the migration is inwards because of the dominating role of the outer Lindblad resonances (Ward 1986,1997). In the second case the planet is more massive, opens up a gap, and the evolution is locked to that of the disk. As a general trend, the disk material drifts inwards on the viscous time scale, and so does the planet. It follows that

$$\left(\frac{{\rm d}a}{{\rm d}t}\right)_{\rm II} = -\frac{3\,\nu}{2\,a} = -\frac{3}{2}\,\alpha \, h^2 \, a\, \,\Omega_{\rm p},$$ (25)

which is known as type II migration. Comparing Eq. (24) with Eq. (25), it turns out that type I drift is faster than type II (i.e., faster than viscous diffusivity) whenever

$$q \gtrsim 3\,\alpha \, h^5 \, \left(\frac{\mbox{$M_\bigstar$ }}{\pi\,a^2\,\Sigma}\right)\cdot$$ (26)

The parameter values employed here yield a right-hand side equal to $2.8\times10^{-6}$, which is just a bit smaller than an Earth-mass (0.93  $M_{\oplus}$). Fast type I migration should continue till the planet grows enough to impose a gap on the disk. By that time, however, it could have already impacted the parent star. Once entered this fast drifting regime, it seems that the planet may survive only if the growth time scale $\tau_{\rm G}\equiv \mbox{$M_{\rm p}$ }/\dot{M}_{\rm p}$is much smaller than the migration time scale $\tau_{\rm M}\equiv a/\vert\dot{a}\vert$.

Figure 20: Migration time scale $\tau _{\rm M}$ versus the mass ratio q. Open triangles represent the results from models CIRO, PEPP, GINO2, and GINO3. The total disk torque $\mathcal{T}_{\rm D}$ is computed assuming $\beta =2\,\eta$, for each model. This means that the region lying inside $\vert\vec{r}-\vec{r}_{\rm p}\vert= 2\,\eta$ is not taken into account. In order to express $\tau _{\rm M}$ into physical units, we suppose that the planet orbits at a=5.12 AU in a disk with mass $\mbox{$M_{\rm D}$ }= 3.5 \times 10^{-3}$  $M_\bigstar$. The solid curve over-plotted represents predictions of the analytical theory as formulated by Ward (1997), for the case of "strong'' viscosity ( $\alpha =4\times 10^{-3}$ in our case) and accounting only for Lindblad torques. It is derived assuming an unperturbed constant surface density and a disk temperature dropping as 1/r. The behavior of this curve reduces to Eq. (24) and Eq. (25) letting $q \rightarrow 0$ and $q \rightarrow \infty$, respectively. The meaning of the dashed line is discussed in the text. The vertical line marks the value given by Eq. (26).

Open with DEXTER

In Fig. 20 the migration time scale $\tau _{\rm M}$ is shown for the main models listed in Table 1. The drifting motion is directed inwards in all cases. The lowest migration velocity belongs to the Earth-mass planet (PEPP1). The second lowest drifting velocity is that of the Jupiter-mass planet (CIRO1). The most rapidly migrating planet is the one having $\mbox{$M_{\rm p}$ }= 20$  $M_{\oplus}$ (PEPP4).

In agreement with predictions of analytical theories, the drift velocity $\vert\dot{a}\vert$ increases for increasing planet mass, just as prescribed by type I migration (Eq. (24)). The fast speed branch has a turning point around $q\simeq 6\times 10^{-5}$, after which migration slows down considerably. Past this point, $\vert\dot{a}\vert$ drops as the planet mass increases. According to the linear theory, this property announces the transition to type II migration (Eq. (25)). As a comparison, the complete theoretical behavior of $\tau _{\rm M}$ (solid line) is also reported in Fig. 20. It was derived by Ward (1997) for viscous disks with $\alpha\gtrsim 10^{-4}$. Equations (24) and (25) represent the asymptotic branches of this curve for very light and very heavy planets, respectively. We also plot (long-dash line) the recent outcome by Tanaka et al. (2002). They derived a general formula which is comprehensive of both Lindblad and corotation torques, in 2D as well as in 3D. Their linear theory, which was developed for (vertically and radially) isothermal and inviscid disks, is applicable as long as the Hill radius of the planet is much smaller than the local disk scale height. Hence we can refer to it only in the Type I migration regime.

In all the models under examination, numerical simulations predict a slower drift than Ward's analytic theory does. Roughly, migration is two times as slow for models GINO2, GINO3, PEPP3, and PEPP4; three times as slow for CIRO1, CIRO2, CIRO3, and PEPP2; planets as massive as 3.2 and 64 Earth-masses migrate on the viscous time scale of the disk (Eq. (25)). Such discrepancies are likely to arise also because the theoretical (solid) curve in Fig. 20 does not include corotation torques. Yet, we have seen that just these torques slow down inward migration. This can be inferred also by comparing the two linear analyses presented in Fig. 20.

It is generally assumed that the smaller the perturbations, the more suitable linear theories are in describing disk-planet interactions. Consequently, the migration time scale we compute for PEPP1 (1  $M_{\oplus}$) should match as close as possible the theoretical prevision. Compared to Ward's prediction, $\tau _{\rm M}$ is six times as large, whereas the agreement with the formula by Tanaka et al. (2002) is complete. Between 3 and 10  $M_{\oplus}$, our simulations predict slightly faster migration velocities. This could be attributed to the inviscid approximation and to non-linearity effects. In this respect, Myioshi et al. (1999) found that, at 5 AU, non-linear effects appear whenever $\mbox{$M_{\rm p}$ }\gtrsim 5$  $M_{\oplus}$. The choice of $\beta$ in the numerical experiments performed here may play some role as well.

It's worthy to note here that the relevance of coorbital torques, to the orbital evolution of a protoplanet, was already pointed out by Ward (1993). A further increment of migration times is to be expected if a full three-dimensional treatment of the system is achieved, especially when low mass planets are dealt with. Local simulations by Miyoshi et al. (1999), global simulations by Kley et al. (2001), and the linear wave theory by Tanaka et al. (2002) indicate that effects due to the disk thickness contribute to reduce the magnitude of the overall torque experienced by the planet. In a flat disk, gas is more efficient in generating torques because it constantly dwells in the planet's orbital plane. In a 3D geometry, material near the vertical of the planet exerts almost a vertical force which does not affect the z-component of the torque in Eq. (20).

   
5.4.1 $\beta$-dependence

Figure 21: Migration time $\tau _{\rm M}$ as a function of $\beta$, the radius of the region excluded from the computation of the torque. The ratio $v=\dot{a}/\vert\dot{a}\vert$ indicates the direction of the planet drift: v<0 for inward migration. The semi-major axis a has the same value as in Fig. 20. The disk mass is cast into the form $\mbox{$M_{\rm D}$ }= k \times 10^{-3}$  $M_\bigstar$. Since in this paper we are using the value k=3.5, the factors in round brackets become $2.9 \times 10^5$ and $2.9 \times 10^4$ years, respectively. The solid, vertical, line indicates the length $\beta =2\,\eta$(see Fig. 17).

Open with DEXTER

In the above discussion, as explained in Sect. 5.3, we do not account for the torque exerted by matter closer than $\beta =2\,\eta$ to the planet. Now we would like to relax this assumption and consider how material, lying even closer, affects the overall torque $\mathcal{T}_{\rm D}$. We mentioned already that the tendency of the nearby gas is to increase the angular momentum of the planet. As we enter the core dominated zone, such tendency may grow stronger and stronger, reducing more and more the magnitude of the negative torque exerted by the rest of the disk.

An overview of the effects, due to nearby matter, on the migration time scale $\tau _{\rm M}$, is given in Fig. 21. Actually, these plots show the dependence of $v/\tau_{\rm M} \propto \mathcal{T}_{\rm D}$ upon $\beta$, where $v=\dot{a}/\vert\dot{a}\vert$ indicates the direction of the planet's migration. The distance where $\beta =2\,\eta$ is marked with a vertical, solid line. $\mathcal{T}_{\rm D}$ is directly proportional to the mass of the disk $M_{\rm D}$. Therefore, to remove this potential restriction, in these plots we let it as a free parameter and write $M_{\rm D}$ as $k \times 10^{-3}$  $M_\bigstar$, although we use k=3.5 for our estimates. In Fig. 21, $\mathcal{T}_{\rm D}=\mathcal{T}_{\rm D}(\beta)$is presented for all of the three reference models.

In CIRO1 (top panel), $\mathcal{T}_{\rm D}$ becomes larger as $\beta$ gets smaller. The sign of the total torque changes around $\beta=\eta$, the threshold of the density core. As a comparison, the behavior of $\mathcal{T}_{\rm D}$, versus $\beta$, is reported also for the model ELEN1. In this case, because of the closed inner radial border, the amount of matter inside the orbit of the planet is five times as large as that of CIRO1. Outside the Hill circle, the torque exerted by the inner-disk, in case of ELEN1, is twice as large as that measured in CIRO1. Instead, torques arising from the outer-disk nearly coincide. Inside the Hill circle, in ELEN1, the contribution to  $\mathcal{T}_{\rm D}$ is relatively small down to $\sim$0.1 $R_{\rm H}$ whereas, in CIRO1, it never appears to be negligible.

CIRO2 (middle panel) behaves somewhat differently from CIRO1. The total torque attains a minimum around $\beta= 0.15$ $R_{\rm H}$, where $\tau_{\rm M}\simeq 3 \times 10^4$ years. Then the positive torques, exerted by close matter, increase the total torque, though it remains negative all the way down to $\beta=0.03$ $R_{\rm H}$. Below such value, $\mathcal{T}_{\rm D}$ diverges positively. Results from the higher resolution model, WPRO1, do not differ significantly (dashed line in Fig. 21).

$\mathcal{T}_{\rm D}$ varies smoothly, as a function of $\beta$, in case of CIRO3. The torques arising from the region enclosed between $\beta\approx 0.2$ $R_{\rm H}$ and $\beta\approx 0.4$ $R_{\rm H}$ almost cancel out, so that the total torque appears nearly constant ( $\tau_{\rm M}\approx 10^5$ years). At shorter distances, positive torques prevail over the negative ones and $\mathcal{T}_{\rm D}$ starts to increase. The sign of the total torque reverses at $\beta\simeq 0.12$ $R_{\rm H}$. For example, at $\beta=\eta$, its value is quite positive, imposing an outward migration rate $\tau_{\rm M}\simeq 3 \times 10^4$ years.

Some comments should be devoted to how the smoothing length affects the total torque. We did not try to reduce further its value, however we ran some models, identical to CIRO1, but multiplying $\lambda$ (see Eq. (19)) by some integer number. A larger smoothing length tends to smear out more the surface density nearby the planet. Besides, it also causes the material to be distributed more symmetrically around it. Both tendencies contribute to reduce the magnitude of the net torque arising from a region with radius ${\approx} \lambda$.

   
5.5 Circumplanetary disk: gas flow

Figure 22: Contour lines of velocity ratios $w_{\rm in}/w_{\rm in}^{\rm D}$ (top panels) and $w_{\rm rot}/w_{\rm rot}^{\rm K}$ (bottom panels). From left to right, we report the above quantities, as computed on the highest level, of models: CIRO1, CIRO2, and CIRO3. $w_{\rm in}>0$ contours indicate places where the flow approaches the planet. At locations where the in-fall velocity is zero, the flow performs a pure rotation around the planet.

Open with DEXTER

In Sect. 5.1 we described in detail how gas flows into the Roche lobe of a planet, for three particular values of its mass.

If a planet is massive enough, say $\mbox{$M_{\rm p}$ }\gtrsim 10$  $M_{\oplus}$, the streams of matter, entering the Roche lobe, produce strong shock waves which then rule the gas flow inside this region. Material passing through the shock fronts is deflected towards the planet, tightening its orbit on it. Less massive planets are not able to cause strong perturbations inside the Roche lobe. As a consequence, the flow pattern appears more uniform around the planet.

Now we would like to investigate quantitatively the rotational regime of the gas inside the circumplanetary disks. In particular, we would like to estimate how much it resembles a Keplerian one. In order to address this issue we decompose the local velocity field $\vec{u}$, in two components, representing the in-fall velocity $w_{\rm in}$ and the rotational velocity $w_{\rm rot}$ of the fluid relative to the planet. The first component is defined as:

$$w_{\rm in} = - \vec{u}\cdot \frac{\vec{r}-\vec{r}_{\rm p}}% {\vert\vec{r}-\vec{r}_{\rm p}\vert},$$

which is positive if a fluid element moves towards the planet. The quantity $w_{\rm rot}$ is the projection of $\vec{u}$ along the direction orthogonal to $\vec{r}-\vec{r}_{\rm p}$ and such that $(\vec{w}_{\rm rot} \times \vec{w}_{\rm in}) \cdot \hat{\vec{z}}$is positive. With this choice $w_{\rm rot}$ is positive for a counter-clockwise rotation.

If circumplanetary disks were regular accretion disks, we should expect them to be in a "Keplerian'' regime. This is characterized by the rotational velocity

$$w_{\rm rot}^{\rm K}= \sqrt{\frac{G\,\mbox{$M_{\rm p}$ }}{\vert\vec{r}-\vec{r}_{\rm p}\vert}},$$ (27)

and the inward viscous diffusion

$$w_{\rm in}^{\rm D}= \frac{3\,\nu}{2\,\vert\vec{r}-\vec{r}_{\rm p}\vert}\cdot$$ (28)

Equations (27) and (28) don't include the smoothing length, $\delta$, because its effects were checked to be unimportant. Comparing $w_{\rm rot}$ and $w_{\rm in}$ with Eq. (27) and Eq. (28), respectively, we can estimate how much the circumplanetary flow is close to be Keplerian, i.e. close that of an unperturbed viscous disk.

Figure 22 shows, for CIRO-models, the contour lines of $w_{\rm in}$ normalized to $w_{\rm in}^{\rm D}$ (top panels) and $w_{\rm rot}$ normalized to $w_{\rm rot}^{\rm K}$ (bottom panels).

As first remark we note that, if we compare Fig. 22 to Fig. 14, lines of equal surface density perturbation are also lines of equal velocity perturbation, as spiral wave theory predicts.

From the top panels of Fig. 22, we can see that material approaches the planet along well defined patterns. Contours $w_{\rm in}=0$ mark locations where the flow rotates around the planet without altering its distance from it. They also separate regions in which material proceeds towards the planet from those where it moves away. One of these contours runs along the spiral ridges. Across it, the in-fall velocity changes abruptly its sign.

The ratio $w_{\rm in}/w_{\rm in}^{\rm D}$ becomes smaller as the gas comes closer to the planet. Since the viscous diffusion $w_{\rm in}^{\rm D}$ is not related to $M_{\rm p}$, it's possible to compare the magnitude of $w_{\rm in}$ for the different cases. Contour level values indicate that it gradually reduces as $M_{\rm p}$ gets smaller.

As regards the rotational component of the velocity field $w_{\rm rot}$ (Fig. 22, bottom panels), we can see that it is generally slightly below $w_{\rm rot}^{\rm K}$. For the Jupiter-mass case, this can be seen also in Fig. 23. Centrifugal over-balance regions are not present around the smallest planet ( $\mbox{$M_{\rm p}$ }=3.2$  $M_{\oplus}$). Instead, they are observed in CIRO1, at $(x_{\rm H},y_{\rm H}) \approx (-0.1,0.3)$ and (0.15,0.1), and in CIRO2, where they are labeled. In both cases, anyway, the ratio $w_{\rm rot}/w_{\rm rot}^{\rm K}$is very close to unity. Centrifugal under-balance is mainly established along spiral ridges. This is in agreement with the idea that spiral arms are zones of compression, hence pressure plays a more active role in supporting the gas.

Figure 23: The plot shows the ratio $w_{\rm rot}/w_{\rm rot}^{\rm K}$at $\varphi =\varphi _{\rm p}$ ( $y_{\rm H}=0$), as computed on the grid level l=6 of the model ELEN2 ( $\mbox{$M_{\rm p}$ }=1$  $M_{\jupiter}$). The core area, in the Jupiter-mass case, extends for $[0.08 \times 0.08]$ $R_{\rm H}$2 (see Fig. 17, left panel). The resolution, in ELEN2, is such that this is covered by ${\sim } 12\times 12$ grid cells.

Open with DEXTER

In all of the three cases shown in Fig. 22, $w_{\rm rot}$ reveals somewhat a circular symmetry only within a distance ${\sim} \eta$from the planet. Yet, we have to notice that this coincides nearly with the region from which matter is removed to simulate the gas accretion. Figure 23 also indicates that the core material (Sect. 5.2.1) rotates slower, when it approaches its center.

   
5.6 Accretion onto the planet

	Figure 24: Mass accretion rate onto the planet versus time, in units of 10-4 disk masses per orbital period of the planet. For $M_{\rm D}=3.5 \times 10^{-3}$  $M_\bigstar$ and a = 5.2 AU, one dimensionless unit corresponds to $2.95 \times 10^{-5}$  $\mbox{$M_{\jupiter}$ }\,{\rm yr}^{\rm -1}$. The initial accretion rate is very small for model CIRO1 because of the imposed initial gap.
Open with DEXTER

Gas matter closely orbiting the planet is eligible to be accreted once its distance, $\vert\vec{r}-\vec{r}_{\rm p}\vert$, is less than $\kappa_{\rm ac}\approx 9 \times10^{-2}$ $R_{\rm H}$. The details of this process are described in Sect. 4.2 (see also Fig. 10). In general, the mass accretion rate of a planet, $\dot{M}_{\rm p}$, becomes relatively constant after the gap (if any) has evolved to a quasi-stationary state.

For a Jupiter-mass planet (CIRO1) this happens around 100 orbital periods, as indicated by the solid line in Fig. 24. The theoretical gap imposed at the beginning of the evolution (see Fig. 11) is deeper and wider, at $\varphi =\varphi _{\rm p}$, than it is later on. For this reason $\dot{M}_{\rm p}$is negligibly small at early evolutionary times. The partial replenishment is related to the formation of the circumplanetary disk which supplies matter for the accretion process.

Smaller planets dig narrower and shallower gaps so this quasi-steady regime is reached even earlier. For both CIRO2 (Fig. 24, short-dash line) and CIRO3 (long-dash line), $\dot{M}_{\rm p}$ reduces a little as the evolution proceeds. This is likely due to the depletion of the inside-orbit disk. In case of CIRO1, most of the inner-orbit material is cleared out during the transitional phase (75% after 100 orbits). Therefore, it does not contribute much during the quasi-steady phase.

In Fig. 25, the accretion rate is plotted against the mass ratio (left panel), in dimensionless units as in Fig. 24. $\dot{M}_{\rm p}$ increases as the planet mass increases and reaches a maximum around 0.5  $M_{\jupiter}$ (model GINO3). Beyond this point the curve starts to drop. In fact, the accretion rate of CIRO1 lies between those of CIRO2 ( $\mbox{$M_{\rm p}$ }= 32$  $M_{\oplus}$) and GINO2 ( $\mbox{$M_{\rm p}$ }= 64$  $M_{\oplus}$). The accretion rates for planets with masses above about 10^-4  $M_\bigstar$ (32  $M_{\oplus}$) are consistent with those obtained by recent models studying the evolution of giant protoplanets (Tajima & Nakagawa 1997). For smaller masses our accretion rates are substantially higher. Two-dimensional approximation may be partially responsible for this since all the gas is confined to the disk midplane. Actually, Kley et al. (2001) showed that accretion rates in disks with thickness are smaller than those in infinitesimally thin disks as soon as the scale height of the disk becomes comparable with the Hill radius of the planet. However, all the more detailed studies of protoplanetary evolution are spherically symmetric (see the review by Wuchterl et al. 2000), and accretion via a planetary accretion disk may allow, anyway, for higher rates.

We did not perform computations involving planets heavier than one Jupiter-mass so we cannot follow the trend of the curve for larger values of q. However, Lubow et al. (1999) found that $\dot{M}_{\rm p}$ decreases in the mass range from 1 to 6  $M_{\jupiter}$.

In the right panel of Fig. 25 the growth time scale $\tau_{\rm G}\equiv \mbox{$M_{\rm p}$ }/\dot{M}_{\rm p}$ is plotted versus q. Logarithmic scaling of the axes shows that the curve decreases almost linearly with respect to the mass of the planet. If we perform a linear least-square fit of the values for which $q \ge 2\times 10^{-4}$, corresponding to $\mbox{$M_{\rm p}$ }= 64$  $M_{\oplus}$, we get:

$$\frac{{\rm d}\log(1/\tau_{\rm G})}{{\rm d}\log(q)} = -0.66.$$

This equation yields:

$$\tau_{\rm G} \propto q^{0.66} \simeq q^{2/3}.$$ (29)

At higher values of q, the curve steepens, decreasing more rapidly. Taking into account the growth time scales relative to the two most massive planets, one finds

$$\tau_{\rm G} \propto q^{1.34} \simeq q^{4/3}.$$ (30)

Therefore, as a planet becomes more massive, the growth time scale grows with an increasing power of its mass. Thus, we can argue that very high mass planets should be extremely rare.

Figure 25: Left panel. Mass accretion rate $\dot{M}_{\rm p}$ as function of the mass ratio q. Quantities are normalized to the disk mass $M_{\rm D}$. Right panel. Dependence of the growth time scale $\tau_{\rm G}\equiv \mbox{$M_{\rm p}$ }/\dot{M}_{\rm p}$ on the planet mass (filled triangles). $\mbox{$M_{\rm D}$ }= 3.5 \times 10^{-3}$  $M_\bigstar$ is employed to express $\tau _{\rm G}$ into years. The two-branch solid line represents the first-order polynomial fitting the data (see Eqs. (29) and (30)).

Open with DEXTER

5.6.1 Influence of numerics

Finally, we comment on the influence of some numerical parameters upon the mass accretion rate. Table 3 lists some results, after the same number of orbits, for various models which can be useful to the goal. The value of $\dot{M}_{\rm p}$ may depend on the mass evacuation rate $\kappa_{\rm ev}$as well as on the radius $\kappa_{\rm ac}$ of the accretion region. By comparing models CIRO2 and WPRO2, it turns out that doubling $\kappa_{\rm ev}$, the accretion rate is only 11% higher. In order to estimate how relevant the extension of the accretion region may be, we ran the model ELEN3 for which we adopted a length $\kappa_{\rm ac}$ equal to 0.6 times that used for GINO3 (the rest of the model set-up being identical). We obtain an accretion rate 8% smaller.

Further, $\dot{M}_{\rm p}$ could depend on the numerical resolution as well. Indeed, this dependency turns out to be very weak, as indicated by two sets of models. ELEN2 and WPRO1 have, in the accretion region, a resolution two times as high as that of CIRO1 and CIRO2, respectively (see Table 1). Despite this fact, in the first case, the accretion rates differ by just 2%, whereas the difference amounts to 4% in the second.

As mentioned in Sect. 4.4, the planet in not symmetrically centered within a grid cell. Model GINO1 was designed to accomplish a fully symmetric configuration (as already explained in Sect. 4.4). However, the planet accretion rate is not affected by this position shift, as can be seen by comparing the values in Table 3, belonging to CIRO3 and GINO1.

The accretion procedure we use in the present work is the same as in Kley et al. (2001). Their evaluations of $\dot{M}_{\rm p}$, for two-dimensional models with $\mbox{$M_{\rm p}$ }=0.5$ and 1  $M_{\jupiter}$, are roughly twice as high as the estimates given in the left panel of Fig. 25. Yet, in that paper, though $\kappa_{\rm ev}$ was set to 5, $\kappa_{\rm ac}$ is equal to 0.5 $R_{\rm H}$, thus the accretion region exceeds by 40 times that used in this paper for models GINO3 and CIRO1.

 

 

Table 3: Mass accretion rate onto the planet for different numerical parameters. This quantity is reported at the same evolutionary time for similar models. $\dot{M}_{\rm p}$ is given in units of disk masses per orbit.

Model	q	$\dot{M}_{\rm p}$	Orbits
CIRO1	$1.0 \times 10^{-3}$	$8.7 \times 10^{-5}$	200
ELEN2	$1.0 \times 10^{-3}$	$8.9 \times 10^{-5}$	200
GINO3	$5.0 \times 10^{-4}$	$1.2 \times 10^{-4}$	200
ELEN3	$5.0 \times 10^{-4}$	$1.1 \times 10^{-4}$	200
CIRO2	$1.0 \times 10^{-4}$	$9.0 \times 10^{-5}$	280
WPRO1	$1.0 \times 10^{-4}$	$8.6 \times 10^{-5}$	280
WPRO2	$1.0 \times 10^{-4}$	$1.0 \times 10^{-4}$	280
CIRO3	$1.0 \times 10^{-5}$	$4.5 \times 10^{-5}$	120
GINO1	$1.0 \times 10^{-5}$	$4.5 \times 10^{-5}$	120

   
6 Conclusions

A number of numerical simulations concerning disk-planet interactions have been performed to get new insights into the scenario of the joint evolution of protoplanets and their environment. They have confirmed analytical theories for gap formation and planet migration.

However, many open questions still remain. The most important unsolved issue is the influence of the ambient gas on the dynamical evolution of a planet. Another one is the way disk-planet interaction changes when small planets, in the mass range of Neptune and Earth, are considered.

We began to investigate in both directions by means of a nested-grid technique, which is particularly suitable for treating these problems. The main asset of this numerical scheme is the possibility of achieving, locally, a very high spatial and temporal resolution. With such a method we are able to resolve very accurately both the inner parts of the Hill sphere of the planet and the global structure of the disk; consequently we can treat the whole azimuthal extent of the disk.

As regards the issues mentioned above, with the present paper we tackled some outstanding problems concerning the growth and migration of protoplanets, covering a range from one Jupiter-mass down to one Earth-mass. Thus, even though we do not include the detailed energetic balance of the planetary structure, which is still beyond present day computer facilities, this study represents a definite improvement in the determination of the torque balance on protoplanets.

Our main achievements can be summarized as follows.

1.

Inside the Hill sphere of the planet, a circumplanetary disk forms. Within it, strong spiral density waves develop, if the planet mass is greater than $\sim$5  $M_{\oplus}$. These waves assume the shape of a two-arm pattern. The two spiral arms are slightly asymmetric with respect to the planet. For decreasing planet masses, they stretch and shorten. Matter is observed to pile up at the location of the planet, generating a very high density zone (we named as density "core''), which might represent its primordial gaseous envelope.

2.

Nearby material exerts positive torque on the planet, slowing down, considerably in some cases, its inward migration. Most of these torques arise from corotation regions, i.e. from gas lying on the planet's orbit. Not all analytical models about migration include them. This can be one of the reasons why our estimates of the migration time scales give somewhat higher values than those predicted by Ward's 1997-theory. On the other hand, Tanaka et al. (2002) account also for corotation torques, but non-linearity and viscosity may play some role in the little discrepancies between their formula and our results.

3.

Within a distance of $\sim$0.1 $R_{\rm H}$ from the planet, the point-mass approximation becomes too restrictive and maybe not appropriate. Therefore, the structure of the planet should be also taken into account over such a length scale. This is absolutely necessary if one wants to evaluate how much of the angular momentum, transferred by closely orbiting matter, is conveyed to the spin of the planet rather than to its orbital angular momentum.

4.

The Keplerian rotational regime of circumplanetary disks is affected by spiral perturbations. Just as for the mass density, the more massive the planet is, the stronger such perturbations are. Gas material, passing through the spiral fronts, is deflected towards the planet. Instead, in the inter-arm regions it moves away from it. This is in analogy with the spiral wave theory in galaxies. Around Earth-mass planets, the rotation of the gas is very slow if compared to the Keplerian rotation. In fact, in this particular case, the density core has nearly a hydrostatic structure.

5.

The mass accretion rate, as a function of the mass of the planet, has a maximum around $\mbox{$M_{\rm p}$ }=0.5$  $M_{\jupiter}$. As long as $\mbox{$M_{\rm p}$ }\lesssim 0.2$  $M_{\jupiter}$, the growth time scale of a planet increases, approximatively, as $\mbox{$M_{\rm p}$ }^{2/3}$. For more massive planets, it increases roughly as $\mbox{$M_{\rm p}$ }^{4/3}$. Such a dependence may contribute to limit the size of a massive planet.

As we have just started to explore these new grounds, each of the items above may deserve a more specific and dedicated study.

Further efforts should be devoted to refining the physical model, especially in the vicinity of the planet, and to approach high resolution three-dimensional computations. Thereupon, some of the future developments could be:

• including an energy equation, by implementing an approximate treatment of radiative transfer and viscous dissipation;
• improving the equation of state by using an alternative form which accounts for the planet's structure;
• evaluating possible effects due to the two-dimensional approximation of the disk, via nested-grid three-dimensional simulations.

Acknowledgements

We would like to thank Udo Ziegler for making the FORTRAN Version of his code NIRVANA available to us. We would like to thank also Dr. Ewald Müller for stimulating discussions on nested-grid calculations. This work was supported by the German Science Foundation (DFG) under grant KL 650/1-1. The numerical computations have been carried out at the Computer Center of the University of Jena.

References

 
• Armitage, J. P. 2001, MNRAS, submitted In the text  
• Artymowicz, P. 1992, PASP, 104, 769 In the text NASA ADS  
• Bodenheimer, P., & Pollack, J. B. 1986, Icarus, 67, 391 In the text NASA ADS  
• Bryden, G., Chen, X., Lin, D. N. C, Nelson, R. P., & Papaloizou, J. C. B. 1999, ApJ, 514, 344 In the text NASA ADS  
• Burkert, A., & Bodenheimer, P. 1993, MNRAS, 264, 798 In the text NASA ADS  
• Butler, R. P., Marcy, G. W., Fischer, D. A., et al. 2001, to appear in Planetary Systems in the Universe: Observations, Formation and Evolution, ed. A. Penny, P. Artymowicz, A.-M. Lagrange, & S. Russell, ASP Conf. Ser. In the text  
• Ciecielag, P., Plewa, T., & Rózyczka, M. 2000, in Disks, Planetesimals and Planets, ed. F. Garzon, C. Eiroa, D. de Winter, & T. J. Mahoney, ASP Conf. Ser., 219, 45 In the text  
• Ciecielag, P., Plewa, T., & Rózyczka, 2000, Astronomische Nachrichten, 321, 171 In the text NASA ADS  
• Goldreich, P., & Tremaine, S. 1980, ApJ, 241, 425 In the text NASA ADS  
• Kley, W. 1998, A&A, 338, L37 NASA ADS  
• Kley, W. 1999, MNRAS, 303, 696 In the text NASA ADS  
• Kley, W. 2000, MNRAS, 313, L47 In the text NASA ADS  
• Kley, W., D'Angelo, G., & Henning, Th. 2001, ApJ, 547, 457 In the text NASA ADS  
• Lubow, S. H., Seibert, M., & Artymowicz, P. 1999, ApJ, 526, 1001 In the text NASA ADS  
• Miyoshi, K., Takeuchi, T., Tanaka, H., & Ida, S. 1999, ApJ, 516, 451 In the text NASA ADS  
• Nelson, R., Papaloizou, J. C. B., Masset, F., & Kley, W. 2000, MNRAS, 318, 18 In the text NASA ADS  
• Ruffert, M. 1992, A&A, 265, 82 In the text NASA ADS  
• Tajima, N., & Nakagawa, Y. 1997, Icarus, 126, 282 In the text NASA ADS  
• Tanaka, H., Takeuchi, T., & Ward, W. 2002, ApJ, 565, 1257 In the text NASA ADS  
• van Leer, B. 1977, J. Comput. Phys., 23, 276 In the text NASA ADS  
• Ward, W. R. 1986, Icarus, 67, 164 In the text NASA ADS  
• Ward, W. R. 1993, Anal. NY Acad. Sci., 675, 314 In the text  
• Ward, W. R. 1997, Icarus, 126, 261 In the text NASA ADS  
• Wuchterl, G. 1993, Icarus, 106, 323 In the text NASA ADS  
• Wuchterl, G., Guillot, T., & Lissauer, J. L. 2000, Protostars and Planets IV, 1081 In the text  
• Yorke, H. W., Bodenheimer, P., & Laughlin, G. 1993, ApJ, 411, 274 In the text NASA ADS  
• Yorke, H. W., & Kaisig, M. 1995, Comput. Phys. Commun., 89, 29 In the text  
• Ziegler, U., & Yorke, H. W. 1997, Comput. Phys. Commun., 101, 54 In the text  
• Ziegler, U. 1998, Comput. Phys. Commun., 109, 111 In the text