text:
# Nested-grid calculations of disk-planet interaction

## 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, which 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. (butler2001 (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 little 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 kley1999 (1999); Bryden et al. bryden1999 (1999); Lubow et al. lubow1999 (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 kley2000 (2000); Nelson et al. rnelson2000 (2000)). For a Jupiter-mass planet and typical parameter values for the disk, the obtained orbital decay time is about $`10^5`$ years, which agrees reasonably well with previous estimates based on analytic linear theories (Goldreich & Tremaine gt1980 (1980); Ward ward1997 (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, Cieciela̧g 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 (armitage2001 (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 at 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 layout 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, infinitesimal thin, model using vertically averaged quantities, such as the surface mass density

$$\mathrm{\Sigma }=_{\mathrm{}}^{\mathrm{}}\rho 𝑑z,$$

where $`\rho `$ is the regular density. We work in a cylindrical coordinate system $`(r,\phi ,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 $`M_{\mathrm{}}=1`$ $`M_{}`$. The total mass of the disk $`M_\mathrm{D}`$, within the simulated region, which extends from $`2.08`$ to $`13`$ AU, is $`3.5\times 10^3`$ $`M_{}`$. Embedded in this disk there is a massive protoplanet with a mass $`M_\mathrm{p}`$, which ranges from one Earth-mass ($`M_{\mathrm{}}`$) to one Jupiter-mass ($`M_{\mathrm{}}`$), 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 $`\phi _\mathrm{p}=\pi `$. The angular velocity $`\mathrm{\Omega }`$ of the rotating frame is then given by

$$\mathrm{\Omega }=\mathrm{\Omega }_p=\sqrt{\frac{G(M_{\mathrm{}}+M_\mathrm{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,\phi `$) continuity equation for $`\mathrm{\Sigma }`$ and the Navier-Stokes equations for each of the two components of the velocity field $`u(u_r,u_\phi )`$. Thus, the set of equations read

$$\frac{\mathrm{\Sigma }}{t}+(\mathrm{\Sigma }u)=0,$$
(2)
$$\frac{(\mathrm{\Sigma }u_r)}{t}+(\mathrm{\Sigma }u_ru)=\mathrm{\Sigma }r(\omega +\mathrm{\Omega })^2\frac{P}{r}\mathrm{\Sigma }\frac{\mathrm{\Phi }}{r}+f_r,$$
(3)
$$\frac{[\mathrm{\Sigma }r^2(\omega +\mathrm{\Omega })]}{t}+[\mathrm{\Sigma }r^2(\omega +\mathrm{\Omega })u]=\frac{P}{\phi }\mathrm{\Sigma }\frac{\mathrm{\Phi }}{\phi }+f_\phi .$$
(4)
Here $`\omega =u_\phi /r`$ is the angular velocity and $`P`$ is the vertically integrated (two-dimensional) pressure. The gravitational potential $`\mathrm{\Phi }`$, generated by the protostar and the planet, is given by

$$\mathrm{\Phi }=\mathrm{\Phi }_{\mathrm{}}+\mathrm{\Phi }_\mathrm{p}=\frac{GM_{\mathrm{}}}{|rr_{\mathrm{}}|}\frac{GM_\mathrm{p}}{|rr_\mathrm{p}|},$$
(5)
where $`r_{\mathrm{}}`$ and $`r_\mathrm{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_\phi `$ which give the viscous force per unit area acting in the radial and azimuthal ($`f_\phi /r`$) direction:

$$f_r=\frac{1}{r}\frac{(rS_{rr})}{r}+\frac{1}{r}\frac{S_{r\phi }}{\phi }\frac{S_{\phi \phi }}{r},$$

$$f_\phi =\frac{1}{r}\frac{(r^2S_{r\phi })}{r}+\frac{S_{\phi \phi }}{\phi }.$$

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 $`𝖲`$ are

$$S_{rr}=2\nu \mathrm{\Sigma }\left(\frac{u_r}{r}\frac{1}{3}u\right),$$
(6)
$$S_{\phi \phi }=2\nu \mathrm{\Sigma }\left(\frac{\omega }{\phi }+\frac{u_r}{r}\frac{1}{3}u\right),$$
(7)
$$S_{r\phi }=\nu \mathrm{\Sigma }\left(\frac{1}{r}\frac{u_r}{\phi }+r\frac{\omega }{r}\right),$$
(8)
where the divergence of the velocity field can be written as

$$u=\frac{1}{r}\frac{(ru_r)}{r}+\frac{\omega }{\phi }.$$

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 $`\mathrm{\Sigma }`$ through

$$P=c_\mathrm{s}^2\mathrm{\Sigma }.$$
(9)
The local isothermal sound speed $`c_\mathrm{s}`$ is given here by

$$c_\mathrm{s}=\frac{H}{r}v_{\mathrm{Kep}},$$
(10)
where $`v_{\mathrm{Kep}}=\sqrt{GM_{\mathrm{}}/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}10^7`$ $`M_{}\mathrm{yr}^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_sH`$).

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=M_\mathrm{p}/M_{\mathrm{}}`$ 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 ziegler1997 (1997); Ziegler ziegler1998 (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 (vanleer1977 (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 explicitely. 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. rnelson2000 (2000)) and 3D (Kley et al. kley2001 (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 (ruffert1992 (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. (yorke1993 (1993)) and Burkert & Bodenheimer (burkert1993 (1993)) simulated the collapse of a protostellar cloud. An application to flux-limited radiation hydrodynamics can be found in Yorke & Kaisig (yorke1995 (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.

#### 3.1.1 Basic integration cycle

In our calculations we use the smallest possible refinement ratio:

$$\frac{\mathrm{\Delta }r(l+1)}{\mathrm{\Delta }r(l)}=\frac{\mathrm{\Delta }\phi (l+1)}{\mathrm{\Delta }\phi (l)}=2,$$
(11)
where $`\mathrm{\Delta }r(l)`$ and $`\mathrm{\Delta }\phi (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

$$\mathrm{\Delta }t_1(3)=\mathrm{min}[\mathrm{\Delta }t_1^{\mathrm{CFL}}(3),\frac{1}{2}\mathrm{\Delta }t_1^{\mathrm{CFL}}(2),\frac{1}{4}\mathrm{\Delta }t_1^{\mathrm{CFL}}(1)],$$

where $`\mathrm{\Delta }t_1^{\mathrm{CFL}}(l)`$ represents the time step resulting from the CFL criterion applied to the level $`l`$, after its latest integration. Thus, $`\mathrm{\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, $`\mathrm{\Delta }t_2^{\mathrm{CFL}}(3)`$. Then it can move further in time for an interval

$$\mathrm{\Delta }t_2(3)=\mathrm{min}[\mathrm{\Delta }t_2^{\mathrm{CFL}}(3),\frac{1}{2}\mathrm{\Delta }t_1^{\mathrm{CFL}}(2),\frac{1}{4}\mathrm{\Delta }t_1^{\mathrm{CFL}}(1)].$$
3. Now level $`2`$ can be integrated for a time

$$\mathrm{\Delta }t_3(2)=\mathrm{\Delta }t_1(3)+\mathrm{\Delta }t_2(3)\mathrm{\Delta }t_1^{\mathrm{CFL}}(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 $`\mathrm{\Delta }t_4(3)`$ and $`\mathrm{\Delta }t_5(3)`$. During the sixth step, the level $`2`$ evolves for $`\mathrm{\Delta }t_6(2)=\mathrm{\Delta }t_4(3)+\mathrm{\Delta }t_5(3)`$. Eventually, during the seventh step, the integration of the level $`1`$ is performed using a time step $`\mathrm{\Delta }t_7(1)=\mathrm{\Delta }t_3(2)+\mathrm{\Delta }t_6(2)`$ which is, by construction, smaller than $`\mathrm{\Delta }t_1^{\mathrm{CFL}}(1)`$. The cycle is now complete and each grid level has evolved for the same amount of time $`\mathrm{\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^{l1}`$ times.

#### 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+1`$ to be more accurate than that of the level $`l`$. Therefore, the fine-grid solution replaces the coarse 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`$.

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 ($`\mathrm{\Sigma }^\mathrm{C}`$) is

$$\mathrm{\Sigma }^\mathrm{C}=\frac{_i\mathrm{\Sigma }_i^\mathrm{F}A_i}{_iA_i}=\frac{(\mathrm{\Sigma }_1^\mathrm{F}+\mathrm{\Sigma }_4^\mathrm{F})A_1+(\mathrm{\Sigma }_2^\mathrm{F}+\mathrm{\Sigma }_3^\mathrm{F})A_2}{2(A_1+A_2)}.$$
(12)
In our case, since the advected quantities are the linear radial momentum density $`\mathrm{\Sigma }u_r`$, and the angular momentum density $`\mathrm{\Sigma }r^2\omega `$, these are the interpolated quantities, along with $`\mathrm{\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^\mathrm{C}`$ the coarse value of the linear momentum density to be interpolated and with $`u_i^\mathrm{F}`$ the surrounding fine grid values, we have:

$$u^\mathrm{C}=\frac{(u_1^\mathrm{F}+u_6^\mathrm{F})A_1+(u_2^\mathrm{F}+u_5^\mathrm{F})A_2+(u_3^\mathrm{F}+u_4^\mathrm{F})A_3}{2(A_1+A_2+A_3)}.$$
(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.

In order to guarantee global mass and momentum conservation for the whole hierarchy, we have to make sure that the momentum flux components $`(\mathrm{\Delta }r\mathrm{\Sigma }u_r/\mathrm{\Delta }t,r\mathrm{\Delta }\phi \mathrm{\Sigma }r^2\omega /\mathrm{\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_k^j(l)`$ represents the value of the quantity $`r\mathrm{\Delta }\phi \mathrm{\Sigma }r^2\omega /\mathrm{\Delta }t`$, at the grid-grid interface location, as computed during the $`k`$-th 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_3(2)`$ will be reset as:

$$f_3(2)=\frac{1}{2}\frac{\mathrm{\Delta }t_1(3)_jf_1^j(3)+\mathrm{\Delta }t_2(3)_jf_2^j(3)}{\mathrm{\Delta }t_3(2)}$$
(14)
This correction is accounted for directly while performing the advection of radial and angular momenta.

#### 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 $`l1`$. 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 (ruffert1992 (1992)), that is by employing the monotonised harmonic mean (van Leer vanleer1977 (1977)).

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

$$g_\epsilon =g+\frac{2\epsilon }{\mathrm{\Delta }x}\mathrm{max}[\frac{(gg_\mathrm{L})(g_\mathrm{R}g)}{g_\mathrm{R}g_\mathrm{L}},0],$$
(15)
provided that $`\mathrm{\Delta }x/2\epsilon \mathrm{\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. Fig. 8 shows how the procedure works in the case of the surface density.

As density is cell-centered, $`\epsilon _r`$ and $`\epsilon _\phi `$ 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 $`\epsilon _r`$ or $`\epsilon _\phi `$ 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.

## 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:

$$R_\mathrm{H}=a\left(\frac{M_\mathrm{p}}{3M_{\mathrm{}}}\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_\mathrm{p}=\frac{a}{(1+q)}.$$
(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, $`M_{\mathrm{}}=1`$ $`M_{}`$. 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_{\mathrm{in}},r_{\mathrm{out}}]`$, where $`r_{\mathrm{in}}=0.4`$ and $`r_{\mathrm{out}}=2.5`$. This is covered by a $`142\times 422`$ mesh (main grid), allowing a resolution such that $`\mathrm{\Delta }r(1)=\mathrm{\Delta }\phi (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`$.

### 4.1 Smoothing of the potential

The perturbing action of the planet is exerted via its gravitational potential $`\mathrm{\Phi }_\mathrm{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 $`\mathrm{\Phi }_\mathrm{p}`$ in Eq. (5) as $`\sqrt{|rr_\mathrm{p}|^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)=\mathrm{min}[\frac{R_\mathrm{H}}{5},\lambda (l)],$$
(18)
where

$$\lambda (l)=\sqrt{\mathrm{\Delta }r^2(l)+r_p^2\mathrm{\Delta }\phi ^2(l)}\sqrt{2}\mathrm{\Delta }r(l).$$
(19)
The value of constant part, in Eq. (18), has become a kind of standard in single grid simulations (Kley kley1999 (1999); Lubow et al. lubow1999 (1999); Kley kley2000 (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<sup>1</sup><sup>1</sup>1Using this notation, we refer not only to the radial extent of the accretion region but also to the fact that the region is centered at $`r_\mathrm{p}`$. by $`|rr_\mathrm{p}|\kappa _{\mathrm{ac}}`$. Since mass is removed from the system after each integration step, the evacuation rate depends on an input parameter $`\kappa _{\mathrm{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 _{\mathrm{ev}}=5`$ is used, while we set $`\kappa _{\mathrm{ac}}`$ between $`8.0\times 10^2`$ and $`9.4\times 10^2`$ $`R_\mathrm{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, $`12\times 12`$ cells participate. The lowest number of employed cells is $`8\times 8`$, which is used for the smallest planet ($`M_\mathrm{p}=1`$ $`M_{\mathrm{}}`$). The largest is $`18\times 18`$ and is adopted for a planet with $`M_\mathrm{p}=0.5`$ $`M_{\mathrm{}}`$ (160 $`M_{\mathrm{}}`$). 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 typically $`10^1`$ $`R_\mathrm{H}`$ which is smaller than the radius of the protoplanet, which is assumed to fill its Roche lobe during its growth phase (Bodenheimer & Pollack bodenheimer1986 (1986); Tajima & Nakagawa tajima1997 (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. wuchterl2000 (2000) and references therein, and also Fig. 16 below), hence the inferred accretion rates may still be unreliable.

### 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 $`\phi =0`$ and $`\phi =2\pi `$. We allow matter to flow out of the computational domain at the inner radial border ($`r_{\mathrm{in}}`$) whereas we set reflective boundary conditions at the outer radial border ($`r_{\mathrm{out}}`$). The angular velocity is set equal to the unperturbed Keplerian value $`\mathrm{\Omega }_{\mathrm{Kep}}=\sqrt{GM_{\mathrm{}}/r^3}`$, both at $`r_{\mathrm{in}}`$ and $`r_{\mathrm{out}}`$.

For low-mass planets ($`M_\mathrm{p}10`$ $`M_{\mathrm{}}`$), boundary conditions should not affect much the system evolution because density waves damp before reaching $`r=r_{\mathrm{out}}`$ and are very weak when they reach $`r=r_{\mathrm{in}}`$. For more massive planets, some reflection is seen at the radial outer 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_{\mathrm{in}}`$, the inner disk is partially cleared. The higher the planet mass is the stronger these effects appear.

### 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_\mathrm{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_\mathrm{p},\phi _\mathrm{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.

## 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_{\mathrm{}}`$ down to $`3.2`$ $`M_{\mathrm{}}`$. 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 Artymowicz1992 (1992); Kley kley1999 (1999); Bryden et al. bryden1999 (1999); Kley kley2001 (2001)). An example of large-scale features can be seen in Fig. 12, where the comprehensive result of a nested-grid computation is displayed.

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 ($`|rr_\mathrm{p}|R_\mathrm{H}`$) and small ($`|rr_\mathrm{p}|R_\mathrm{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 ($`M_\mathrm{p}=3.2`$ $`M_{\mathrm{}}`$), a trailing density wave, emanating from the planet, spirals around the star for about $`4\pi `$, vanishing approximately at $`r=2`$. In Pepp1 ($`M_\mathrm{p}=1`$ $`M_{\mathrm{}}`$), 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 ($`M_\mathrm{p}=1`$ $`M_{\mathrm{}}`$) 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`$, $`\phi =3.12`$ and $`r=1.07`$, $`\phi =3.07`$, respectively), where the flow is nearly at rest with respect to the planet<sup>2</sup><sup>2</sup>2They are designated as “X-points” by Lubow et al. (lubow1999 (1999)). 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. (lubow1999 (1999)).

The case $`M_\mathrm{p}=32`$ $`M_{\mathrm{}}`$ (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 $`\phi <\phi _\mathrm{p}`$, and collides with the stream incoming from the other X-point, generating the upper spiral arm (at $`\phi >\phi _\mathrm{p}`$). However, here the locations, from which these gas streams depart ($`r=0.97`$, $`\phi =3.14`$ and $`r=1.03`$, $`\phi =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, $`M_\mathrm{p}=3.2`$ $`M_{\mathrm{}}`$, 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_\mathrm{H}`$, from side to side of the planet at $`\phi \phi _\mathrm{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_\mathrm{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 ($`0.5`$ $`R_\mathrm{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 $`\mathrm{\Sigma }(r,\phi )`$ 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 $`\mathrm{\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 ($`\phi >\phi _\mathrm{p}`$) twists towards the planet more than the other one does. The curvature<sup>3</sup><sup>3</sup>3The curvature of a circle of radius $`R`$ is $`1/R`$. 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 $`r1.001`$ and $`\phi <\phi _\mathrm{p}`$, the density is systematically lower than it is on the opposite side ($`r0.999`$, $`\phi >\phi _\mathrm{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 $`r0.999`$ and $`r1.001`$

#### 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 $`\mathrm{\Delta }M^{\mathrm{ev}}`$ is the evacuated amount of mass, then typically

$$\frac{\mathrm{\Delta }M_{\mathrm{ev}}}{M}\frac{\mathrm{\Delta }\mathrm{\Sigma }}{\mathrm{\Sigma }}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 wuchterl1993 (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_\mathrm{p},\phi _\mathrm{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_\mathrm{H}`$ (see the caption of Fig. 15 for details).

In these two cases, the core width, $`2\eta `$ (exactly defined below), can be estimated to be approximatively $`0.1`$ $`R_\mathrm{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)1.4\mathrm{\Delta }r(ng)`$ (see Eq. 19). From Table 1 one can deduce that $`\delta 10^2`$ $`R_\mathrm{H}`$ (though it changes a little for the different values of $`q`$). Despite that, the core width always equals at least $`6\mathrm{\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 $`\mathrm{\Sigma }=2`$ on the finest level of each model, its extent, at $`r=r_\mathrm{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.

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_{\mathrm{}}`$), 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 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 $`\overline{\mathrm{\Sigma }}`$ is computed over the region $`0.2R_\mathrm{H}|rr_\mathrm{p}|0.5R_\mathrm{H}`$. Then we define the core width ($`2\eta `$) to be that where $`\mathrm{\Sigma }=2\overline{\mathrm{\Sigma }}`$. In Fig. 17 (left panel) the dependence of $`2\eta `$ versus $`q`$ is displayed. The ratio $`\eta /R_\mathrm{H}`$ decreases for increasing values of $`q`$. However, between $`3.2`$ and $`32`$ $`M_{\mathrm{}}`$, 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.

### 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_{\mathrm{}}`$, 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 $`f_g(i,j)`$, due to each grid cell $`(i,j)`$, is calculated. The torque, with respect to the star, is then

$$t(i,j)=r_\mathrm{p}\times f_g(i,j).$$
(20)
Since we are interested in the $`z`$-component of this vector, we perform the scalar product $`\widehat{z}t(i,j)=t_z(i,j)`$, where $`\widehat{z}`$ is the unit vector of the vertical direction. Finally, by summing over all $`i`$ and $`j`$, we obtain the total disk torque $`𝒯_\mathrm{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 $`𝒯_\mathrm{D}`$.

In order to avoid the region dominated by the core, we exclude a certain area from the computation of $`𝒯_\mathrm{D}`$. Because of the way we defined the core radius $`\eta `$, some core material may still lie outside $`|rr_\mathrm{p}|=\eta `$. Therefore, for safety reasons, we choose not to take into account the planet neighborhood defined by $`|rr_\mathrm{p}|\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 $`𝒯_\mathrm{D}`$.

Because of the global angular momentum transfer, the disk material (orbiting the star), at $`r>r_\mathrm{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 $`M_\mathrm{p}=1`$ $`M_{\mathrm{}}`$, such behavior was also found by Lubow et al. lubow1999 (1999)). The radial torque distribution, from Ciro-models, is shown in Fig. 18. The two sets of profiles belong to the grid levels $`l=ng1`$ (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_\mathrm{H}<0`$ is positive while the one coming from $`x_\mathrm{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_\mathrm{H}>0`$ is definitely positive and, in magnitude, almost $`30`$ times as large as that arising from the region $`x_\mathrm{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 $`𝒯(l)`$. This represents the total torque computed over the entire domain but the part covered by the grid level $`l+1`$ and such that $`𝒯(ng)=𝒯_\mathrm{D}`$. Fig. 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 $`2hr_\mathrm{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.

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_\mathrm{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 $`|rr_\mathrm{p}|<0.5`$ $`R_\mathrm{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 ($`M_\mathrm{p}=1`$ $`M_{\mathrm{}}`$). 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 $`0.5hr_\mathrm{p}`$ and $`hr_\mathrm{p}`$, at the starting positions the disk spirals.

#### 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_\mathrm{p}`$ ($`x_\mathrm{H}0`$). Here in fact, $`f_g(i,j)`$ is perpendicular to $`r_\mathrm{p}`$ and the cross-product in Eq. (20) achieves its maximum (minimum). The material leading the planet (at $`\phi >\phi _\mathrm{p}`$ or $`y_\mathrm{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_\mathrm{H})`$ and write their mass density as $`\mathrm{\Sigma }^\pm `$. Then we can write $`|t_z^\pm |r_\mathrm{p}\mathrm{\Sigma }^\pm /y_\mathrm{H}^2`$, which yields:

$$|t_z^+||t_z^{}|r_\mathrm{p}\frac{\mathrm{\Delta }\mathrm{\Sigma }}{y_\mathrm{H}^2},$$
(21)
where $`\mathrm{\Delta }\mathrm{\Sigma }=\mathrm{\Sigma }^+\mathrm{\Sigma }^{}`$. Any mismatch of the surface density, $`\mathrm{\Delta }\mathrm{\Sigma }`$, causes a torque mismatch amplified by an amount equal to $`y_\mathrm{H}^2`$. It’s worth noticing that, on larger distances $`|rr_\mathrm{p}|r_\mathrm{p}`$, the torque mismatch is amplified less, in fact:

$$|t_z^+||t_z^{}|\frac{\mathrm{\Delta }\mathrm{\Sigma }}{r_\mathrm{p}}.$$
(22)
This is the reason why surface density asymmetries near the planet have a very strong impact on $`𝒯_\mathrm{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_\mathrm{H}=0.15`$ and the minimum at $`x_\mathrm{H}0.15`$ are produced at $`y_\mathrm{H}0.1`$ and $`y_\mathrm{H}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 $`𝒯_\mathrm{D}`$, is

$$\frac{da}{dt}=\frac{2𝒯_\mathrm{D}}{M_\mathrm{p}a\mathrm{\Omega }_\mathrm{p}}.$$
(23)
Analytical estimates of $`𝒯_\mathrm{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 (ward1997 (1997)) derives the following expression:

$$\left(\frac{da}{dt}\right)_\mathrm{I}\frac{1}{2}qh^3a\mathrm{\Omega }_\mathrm{p}\left(\frac{\pi a^2\mathrm{\Sigma }}{M_{\mathrm{}}}\right).$$
(24)
The direction of the migration is inwards because of the dominating role of the outer Lindblad resonances<sup>4</sup><sup>4</sup>4This is true as long as the temperature gradient within the disk is negative (Ward ward1986 (1986)). (Ward ward1986 (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{da}{dt}\right)_{\mathrm{II}}=\frac{3\nu }{2a}=\frac{3}{2}\alpha h^2a\mathrm{\Omega }_\mathrm{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

$$q3\alpha h^5\left(\frac{M_{\mathrm{}}}{\pi a^2\mathrm{\Sigma }}\right).$$
(26)
The parameter values employed here<sup>5</sup><sup>5</sup>5If we take into account the dependence of the unperturbed surface density upon the radial distance $`r`$, $`\pi a^2\mathrm{\Sigma }=0.38M_\mathrm{D}`$. yield a right-hand side equal to $`2.8\times 10^6`$, which is just a bit smaller than an Earth mass ($`0.93`$ $`M_{\mathrm{}}`$). 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 _\mathrm{G}M_\mathrm{p}/\dot{M}_\mathrm{p}`$ is much smaller than the migration time scale $`\tau _\mathrm{M}a/|\dot{a}|`$.

In Fig. 20 the migration time scale $`\tau _\mathrm{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 $`M_\mathrm{p}=20`$ $`M_{\mathrm{}}`$ (Pepp4).

In agreement with predictions of analytical linear theories, the drift velocity $`|\dot{a}|`$ increases for increasing planet mass, just as prescribed by type I migration (Eq. 24). The fast speed branch has a turning point around $`q6\times 10^5`$, after which migration slows down considerably. Past this point, $`|\dot{a}|`$ 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 _\mathrm{M}`$ is also reported in Fig. 20. It was derived by Ward (ward1997 (1997)) for viscous disks with $`\alpha 10^4`$. Eqs. (24) and (25) represent the asymptotic branches of this curve for very light and very heavy planets, respectively. The vertical line, in the figure, indicates the mass value above which migration is faster than inward viscous diffusion.

In all the models under examination, numerical simulations predict a slower drift than analytic linear 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). $`\tau _\mathrm{M}`$ is six times as large in case of the Earth-mass planet, when compared to analytic predictions. Such discrepancies are likely to arise also because the theoretical curve in Fig. 20 does not include corotation torques. Yet, we have seen that just these torques slow down inward migration, preventing the total torque $`𝒯_\mathrm{D}`$ from being very negative.

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 (ward1993 (1993)).

#### 5.4.1 $`\beta `$-dependence

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 $`𝒯_\mathrm{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 _\mathrm{M}`$, is given in Fig. 21. Actually, these plots show the dependence of $`v/\tau _\mathrm{M}𝒯_\mathrm{D}`$ upon $`\beta `$, where $`v=\dot{a}/|\dot{a}|`$ indicates the direction of the planet’s migration. The distance where $`\beta =2\eta `$ is marked with a vertical, solid line. $`𝒯_\mathrm{D}`$ is directly proportional to the mass of the disk $`M_\mathrm{D}`$. Therefore, to remove this potential restriction, in these plots we let it as a free parameter and write $`M_\mathrm{D}`$ as $`k\times 10^3`$ $`M_{\mathrm{}}`$, although we use $`k=3.5`$ for our estimates. In Fig. 21, $`𝒯_\mathrm{D}=𝒯_\mathrm{D}(\beta )`$ is presented for all of the three reference models.

In Ciro1 (top panel), $`𝒯_\mathrm{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 $`𝒯_\mathrm{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 $`𝒯_\mathrm{D}`$ is relatively small down to $`0.1`$ $`R_\mathrm{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_\mathrm{H}`$, where $`\tau _\mathrm{M}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_\mathrm{H}`$. Below such value, $`𝒯_\mathrm{D}`$ diverges positively. Results from the higher resolution model, Wpro1, do not differ significantly (dash-line in Fig. 21).

$`𝒯_\mathrm{D}`$ varies smoothly, as a function of $`\beta `$, in case of Ciro3. The torques arising from the region enclosed between $`\beta 0.2`$ $`R_\mathrm{H}`$ and $`\beta 0.4`$ $`R_\mathrm{H}`$ almost cancel out, so that the total torque appears nearly constant ($`\tau _\mathrm{M}10^5`$ years). At shorter distances, positive torques prevail over the negative ones and $`𝒯_\mathrm{D}`$ starts to increase. The sign of the total torque reverses at $`\beta 0.12`$ $`R_\mathrm{H}`$. For example, at $`\beta =\eta `$, its value is quite positive, imposing an outward migration rate $`\tau _\mathrm{M}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 $`\lambda `$.

### 5.5 Circumplanetary disk: gas flow

In Sect. 5 we have 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 $`M_\mathrm{p}10`$ $`M_{\mathrm{}}`$, 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 $`u`$, in two components, representing the in-fall velocity $`w_{\mathrm{in}}`$ and the rotational velocity $`w_{\mathrm{rot}}`$ of the fluid relative to the planet. The first component is defined as:

$$w_{\mathrm{in}}=u\frac{rr_\mathrm{p}}{|rr_\mathrm{p}|},$$

which is positive if a fluid element moves towards the planet. The quantity $`w_{\mathrm{rot}}`$ is the projection of $`u`$ along the direction orthogonal to $`rr_\mathrm{p}`$ and such that $`(w_{\mathrm{rot}}\times w_{\mathrm{in}})\widehat{z}`$ is positive. With this choice $`w_{\mathrm{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_{\mathrm{rot}}^\mathrm{K}=\sqrt{\frac{GM_\mathrm{p}}{|rr_\mathrm{p}|}},$$
(27)
and the inward viscous diffusion

$$w_{\mathrm{in}}^\mathrm{D}=\frac{3\nu }{2|rr_\mathrm{p}|}.$$
(28)
Equations (27) and (28) don’t include the smoothing length, $`\delta `$, because its effects were checked to be unimportant. Comparing $`w_{\mathrm{rot}}`$ and $`w_{\mathrm{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_{\mathrm{in}}`$ normalized to $`w_{\mathrm{in}}^\mathrm{D}`$ (top panels) and $`w_{\mathrm{rot}}`$ normalized to $`w_{\mathrm{rot}}^\mathrm{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_{\mathrm{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_{\mathrm{in}}/w_{\mathrm{in}}^\mathrm{D}`$ becomes smaller as the gas comes closer to the planet. Since the viscous diffusion $`w_{\mathrm{in}}^\mathrm{D}`$ is not related to $`M_\mathrm{p}`$, it’s possible to compare the magnitude of $`w_{\mathrm{in}}`$ for the different cases. Contour level values indicate that it gradually reduces as $`M_\mathrm{p}`$ gets smaller.

As regards the rotational component of the velocity field $`w_{\mathrm{rot}}`$ (Fig. 22, bottom panels), we can see that it is generally slightly below $`w_{\mathrm{rot}}^\mathrm{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 ($`M_\mathrm{p}=3.2`$ $`M_{\mathrm{}}`$). Instead, they are observed in Ciro1, at $`(x_\mathrm{H},y_\mathrm{H})(0.1,0.3)`$ and $`(0.15,0.1)`$, and in Ciro2, where they are labeled. In both cases, anyway, the ratio $`w_{\mathrm{rot}}/w_{\mathrm{rot}}^\mathrm{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.

In all of the three cases shown in Fig. 22, $`w_{\mathrm{rot}}`$ reveals somewhat a circular symmetry only within a distance $`\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

Gas matter closely orbiting the planet is eligible to be accreted once its distance, $`|rr_\mathrm{p}|`$, is less than $`\kappa _{\mathrm{ac}}9\times 10^2`$ $`R_\mathrm{H}`$. The details of this process are described in Sect. 4.2 (see Fig. 10). In general, the mass accretion rate of a planet, $`\dot{M}_\mathrm{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 $`\phi =\phi _\mathrm{p}`$, than it is later on. For this reason $`\dot{M}_\mathrm{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}_\mathrm{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}_\mathrm{p}`$ increases as the planet mass increases and reaches a maximum around $`0.5`$ $`M_{\mathrm{}}`$ (model Gino3). Beyond this point the curve starts to drop. In fact, the accretion rate of Ciro1 lies between those of Ciro2 ($`M_\mathrm{p}=32`$ $`M_{\mathrm{}}`$) and Gino2 ($`M_\mathrm{p}=64`$ $`M_{\mathrm{}}`$). The accretion rates for planets with masses above about $`10^4`$ $`M_{\mathrm{}}`$ ($`30`$ $`M_{\mathrm{}}`$) are consistent with those obtained by recent models studying the evolution of giant proto-planets (Tajima & Nakagawa tajima1997 (1997)). For smaller masses our accretion rates are substantially higher. However, all the more detailed studies of protoplanetary evolution are spherically symmetric (see the review by Wuchterl et al. wuchterl2000 (2000)), and accretion via a planetary accretion disk may allow for higher accretion 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. (lubow1999 (1999)) found that $`\dot{M}_\mathrm{p}`$ decreases in the mass range from 1 to 6 $`M_{\mathrm{}}`$.

In the right panel of Fig. 25 the growth time scale $`\tau _\mathrm{G}M_\mathrm{p}/\dot{M}_\mathrm{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 $`q2\times 10^4`$, corresponding to $`M_\mathrm{p}=64`$ $`M_{\mathrm{}}`$, we get:

$$\frac{d\mathrm{log}(1/\tau _\mathrm{G})}{d\mathrm{log}(q)}=0.66.$$

This equation yields:

$$\tau _\mathrm{G}q^{0.66}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 _\mathrm{G}q^{1.34}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.

#### 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}_\mathrm{p}`$ may depend on either the radius $`\kappa _{\mathrm{ac}}`$ of the accretion region and the mass evacuation rate $`\kappa _{\mathrm{ev}}`$. By comparing models Ciro2 and Wpro2, it turns out that doubling $`\kappa _{\mathrm{ev}}`$, the accretion rate is only 11% higher. In order to estimate how relevant the extension of the accretion region may be, we ran a model having an accretion length scale ($`\kappa _{\mathrm{ac}}`$) $`0.6`$ times as small as that of Gino3 (the rest of the model set-up being identical). We obtain an accretion rate 7% smaller.

Further, $`\dot{M}_\mathrm{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 Wpro2 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 Gino1 and Ciro3.

## 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, because we treat the whole azimuthal extent of the disk.

Concerning 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 larger than $`5`$ $`M_{\mathrm{}}`$. 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 torques 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. Analytical models about migration do not account for them. This can be one of the reasons why our estimates of the migration time scales give somewhat higher values than those predicted by such theories.
3. Within a distance of $`0.1`$ $`R_\mathrm{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 $`M_\mathrm{p}=0.5`$ $`M_{\mathrm{}}`$. As long as $`M_\mathrm{p}0.2`$ $`M_{\mathrm{}}`$, the growth time scale of a planet increases, approximatively, as $`M_\mathrm{p}^{2/3}`$. For more massive planets, it increases roughly as $`M_\mathrm{p}^{4/3}`$. Such a dependence may contribute to limit the size of a massive planet.

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

Next efforts should be devoted to refine the physical model, especially in the vicinity of the planet. Henceforth, 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 three-dimensional simulations.

###### Acknowledgements.

We would like to thank Udo Ziegler for making the FORTRAN Version of his code Nirvana<sup>6</sup><sup>6</sup>6http://www.aip.de/ ziegler/ 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.


format: text

text:
# Additional Stringy Sources of the Electric Dipole Moments

\[

## Abstract

We show that string models with low energy supersymmetry which accommodate the fermion mass hierarchy generally give non-universal soft trilinear couplings (A-terms). In conjunction with the apparently large Cabibbo-Kobayashi-Maskawa (CKM) phase, this results in large fermion EDMs even in the absence of CP violating phases in the supersymmetry breaking auxiliary fields and the $`\mu `$-term. Nonobservation of the EDMs therefore implies that strings select special flavor and/or supersymmetry breaking patterns.

\]

Recent measurements of the CP asymmetry in the $`B\mathrm{\Psi }K_s`$ decay  imply that CP is significantly violated in nature and is not an approximate symmetry. Supersymmetric models with large CP phases generally predict fermion electric dipole moments (EDMs) that far exceed the experimental limits , and this presents a major challenge for low energy supersymmetry, the so called SUSY CP problem.

This problem arises because there is no a priori reason for the SUSY breaking dynamics to conserve CP. In particular, in supergravity models the auxiliary fields that break supersymmetry ($`F_S`$,$`F_T`$, etc.) are in general complex, which leads to complex soft SUSY breaking parameters. These new CP-phases, which are absent in the non-supersymmetric case, induce large EDMs. Conventionally, this is believed to be the source of the CP problem in this class of models.

In this letter, we point out that there is another source of EDMs inherent in more fundamental models such as supersymmetric models derived from strings. It is present even if the SUSY breaking ($`F_S`$,$`F_T`$, etc.) conserves CP, and has its origin in the flavor structures. String theory, as a fundamental theory, has to explain the observed fermion mass hierarchy and mixings. This generally leaves an imprint on the flavor structure of the soft SUSY breaking terms and, when combined with the phenomenological requirement that the Yukawa matrices contain $`𝒪(1)`$ CP-phases, results in unacceptably large EDMs. The source of the EDMs lies in the CP-phases present already in the non-supersymmetric case, i.e. those in the Yukawa matrices, which are rendered observable by supersymmetry.

Non-observation of EDMs provides a strong constraint on the fundamental flavor structure of the Yukawas and/or supersymmetry breaking, in addition to constraints on any $`new`$ CP phases that may occur in the latter, and indicates that the fundamental theory selects quite special patterns of the flavor structures or supersymmetry breaking.

Let us start with the general arguments that imply non-universality in the soft trilinear terms relevant to the EDM calculations. In all supergravity models, the soft SUSY breaking parameters are given in terms of the Kähler potential $`K`$ and the superpotential $`W`$. In particular, the trilinear parameters are written as

$`A_{\alpha \beta \gamma }=F^m\left[\widehat{K}_m+_m\mathrm{log}Y_{\alpha \beta \gamma }_m\mathrm{log}(\stackrel{~}{K}_\alpha \stackrel{~}{K}_\beta \stackrel{~}{K}_\gamma )\right].`$ (1)
Here the Latin indices refer to the hidden sector fields while the Greek indices refer to the observable fields; the Kähler potential is expanded in observable fields as $`K=\widehat{K}+\stackrel{~}{K}_\alpha |C^\alpha |^2+\mathrm{}`$ and $`\widehat{K}_m_m\widehat{K}`$. The sum in $`m`$ runs over all of the SUSY breaking fields. We note that the supergravity notation for the A-terms and Yukawas can be connected to the “usual” one by fixing the order of the indices as follows: the first index is to refer to the Higgs fields, the second to the quark doublets, and the third to the quark singlets, e.g. $`Y_{H_1Q_iD_j}Y_{ij}^d`$.

Since string theory is a fundamental theory, the Yukawa matrices must be generated dynamically. That is, the Yukawa couplings are functions of fields that dynamically get vacuum expectation values leading to the observed fermion masses and mixings (in the same way that the VEV of the dilaton field produces the gauge coupling). Also, since CP is a gauge symmetry in strings  and can be broken only spontaneously, the fields responsible for generating the Yukawa couplings must attain CP-violating VEVs in order to produce the CKM phase. Therefore, the derivatives of the Yukawa couplings in Eq.1 will generally be non-zero and should be taken into account. As we will see, this contribution is often significant even if the corresponding SUSY breaking component $`F^m`$ is small.

To estimate how large these contributions are expected to be, we will consider two possible ways to generate the Yukawa hierarchies: through exponential factors, as in the heterotic string, and through non-renormalizable operators, as for example in type I or other models.

In weakly coupled heterotic orbifolds, the Yukawa couplings are calculable and given in terms of the T-moduli fields. The fermion mass hierarchy has a geometrical origin and appears due to the fact that interaction of the states placed at different orbifold fixed points is exponentially suppressed . The matter fields must belong to the twisted sectors since otherwise the couplings are 1 or 0 and the mass hierarchy cannot be generated at the renormalizable level. The Yukawa coupling of the states at the fixed points $`f_{1,2,3}`$ belonging to the twisted sectors $`\theta _{1,2,3}`$ is given by

$`Y_{f_1f_2f_3}=N{\displaystyle \underset{uZ^n}{}}\mathrm{exp}\left[4\pi T\left(f_{23}+u\right)^TM\left(f_{23}+u\right)\right].`$ (2)
where $`f_{23}f_2f_3`$, $`N`$ is a normalization factor, and the matrix $`M`$ (with fractional entries) is related to the internal metric of the orbifold. Here $`T`$ is normalized such that $`TT+i`$ is a symmetry of the Lagrangian . For a realistic case, Re$`T=𝒪(1)`$, the sum is dominated by one term and

$$_T\mathrm{ln}Y_{f_1f_2f_3}4\pi f_{23}^{}{}_{}{}^{T}Mf_{23}$$
(3)
for some (typically fractional) $`f_{23}`$ depending on the positions of the fixed points. This expession is independent of $`T`$ and can be evaluated for various orbifolds. It vanishes if the fixed points coincide, i.e. there is no suppression of the Yukawa interaction. In other cases it is of order one (or larger, see e.g. ). This creates a significant non-universality in the A-terms unless $`F_T0`$. Indeed, to produce a fermion mass hierarchy we are bound to place different quark generations at different orbifold fixed points. Thus, the expression (3) will be generation dependent. The other contributing terms in Eq.1, $`\widehat{K}_T=3/(T+\overline{T})`$ and $`_T\mathrm{log}(\stackrel{~}{K}_\alpha \stackrel{~}{K}_\beta \stackrel{~}{K}_\gamma )=(n_\alpha +n_\beta +n_\gamma )/(T+\overline{T})`$ are generation independent since the sum of the modular weights $`n_i`$ is fixed by the point group selection rule for the Yukawa interactions. Unless we encounter a special case $`F_T0`$, the non-universality will be present. We thus conclude that the A-term non-universality is a generic feature of realistic heterotic models.

In general, it might be too strict a requirement to demand that all of the required features appear at the renormalizable level. In many string models the renormalizable Yukawa couplings are either unknown or are 0,1 (e.g. Type I models). Thus, to reproduce the observed fermion masses, non-renormalizable operators must be taken into account. The mass hierarchy is then created via powers of a small (in Planck units) VEV of a certain field $`\varphi `$ ,, i.e.

$$Y_{\alpha \beta \gamma }\varphi ^{q_{\alpha \beta \gamma }}.$$
(4)
Here $`q_{\alpha \beta \gamma }`$ is integer and can sometimes be associated with a $`U(1)`$ charge (as in the Froggatt-Nielsen mechanism). Clearly, if this field breaks supersymmetry, $`F_\varphi 0`$, the generated A-terms will be nonuniversal:

$$_\varphi \mathrm{ln}Y_{\alpha \beta \gamma }=\frac{q_{\alpha \beta \gamma }}{\varphi }.$$
(5)
Generally, $`\varphi `$ is expected to give a (small) contribution to supersymmetry breaking. To estimate its natural size, let us recall that the Kähler potential for untwisted fields is given by $`K=\mathrm{ln}(S+\overline{S})3\mathrm{ln}(T+\overline{T}\varphi \overline{\varphi })`$. The resulting auxiliary field is then

$$F_\varphi =m_{3/2}\frac{T+\overline{T}\varphi \overline{\varphi }}{3}\left[\overline{\varphi }\frac{W_T}{W}+\frac{W_\varphi }{W}\right].$$
(6)
A similar result holds for a twisted $`\varphi `$. For small $`\varphi `$, $`F_\varphi `$ is typically of order $`\varphi m_{3/2}`$ (see also e.g. ). As a result, the Yukawa-induced contribution to the A-terms is

$$\mathrm{\Delta }A_{\alpha \beta \gamma }q_{\alpha \beta \gamma }m_{3/2}.$$
(7)
The “charges” $`q_{\alpha \beta \gamma }`$ are generation-dependent and order one, so the resulting non-universality is very significant. Again, we come to the conclusion that realistic models of the fermion masses entail non-universal A-terms. We note that additional nonuniversal contributions may come from the Kähler potential. However, these are not direcly related to the Yukawa structures, so we do not discuss them here.

On general grounds, it is natural to expect that once non-universality is created in the Yukawa interactions, it will also appear in the analogous soft-breaking terms. This leads to large EDMs even if the SUSY breaking dynamics do not generate new CP phases, i.e. $`F_S`$,$`F_T`$, etc. and the $`\mu `$-term are real. To see this let us suppose that initially there are no CP-violating phases in any of the soft breaking parameters. At first sight this seems to avoid overproduction of the EDMs. However this is not the case. Indeed, the relevant trilinear interactions are

$`\mathrm{\Delta }_{\mathrm{soft}}={\displaystyle \frac{1}{6}}A_{\alpha \beta \gamma }Y_{\alpha \beta \gamma }C^\alpha C^\beta C^\gamma .`$ (8)
Denoting $`\widehat{A}_{\alpha \beta \gamma }A_{\alpha \beta \gamma }Y_{\alpha \beta \gamma }`$ and identifying $`\widehat{A}_{H_1Q_iD_j}\widehat{A}_{ij}^d`$, it is easy to see that generally $`\widehat{A}`$ and Y matrices are “misaligned”. In the basis where the Yukawas are diagonal and real (the super-CKM basis), the $`\widehat{A}`$ matrices are not diagonal while their diagonal elements contain complex phases. Under the superfield rotation $`\widehat{U}_{L,R}V_{L,R}^u\widehat{U}_{L,R},\widehat{D}_{L,R}V_{L,R}^d\widehat{D}_{L,R}`$ to the super-CKM basis, i.e. $`Y^uV_L^uY^uV_R^u=\mathrm{diag}(h_u,h_c,h_t)`$, $`Y^dV_L^dY^dV_R^d=\mathrm{diag}(h_d,h_s,h_b)`$, the $`\widehat{A}`$ matrices transform as

$`\widehat{A}^uV_L^u\widehat{A}^uV_R^u,\widehat{A}^dV_L^d\widehat{A}^dV_R^d.`$ (9)
Since the Yukawa matrices contain $`𝒪(1)`$ CP-phases, so do the rotation matrices $`V_{L,R}`$. Therefore, the rotation by itself can induce diagonal CP-phases even if the $`A`$\- or $`\widehat{A}`$\- matrices initially are completely real. Furthermore, due to non-universality the diagonal elements of $`\widehat{A}`$ will have contributions from all three generations, e.g.

$$\widehat{A}_{11}^um_u+ϵm_c+ϵ^{}m_t.$$
(10)
This is to be contrasted with the universal case when the Yukawas and $`\widehat{A}`$ matrices are diagonalized simultaneously and $`\widehat{A}_{ii}`$ are proportional to $`m_i`$. The imaginary parts of the diagonal entries of the $`\widehat{A}`$ matrix in the super-CKM basis induce quark EDMs via the gluino-squark loops with LR mass insertions. Eq.10 implies that the magnitude of the mass insertions contributing to the EDMs, $`(\delta _{LR}^{u,d})_{ii}\widehat{A}_{ii}^{u,d}H_{u,d}/\stackrel{~}{m}^2`$, is significantly enhanced and even a small CP-phase can result in large EDMs.

To illustrate the strength of this effect, we present the scatter plots (Fig.1,2) of the neutron and mercury EDMs versus $`real`$ A-terms. We take representative Yukawa textures (similar to those of Ref. but possessing no symmetry) with $`conservative`$ flavor mixing, i.e. $`(V_{L,R})_{13}V_{td}^{CKM}𝒪(10^2)`$, $`(V_{L,R})_{12}\mathrm{sin}\theta _C`$, and containing order one complex phases. We fix the universal scalar mass $`m_0`$ and $`m_{1/2}`$ to be 200 GeV (i.e. $`m_{\stackrel{~}{q}}m_{\stackrel{~}{g}}`$ 500 GeV at $`M_Z`$) and $`\mathrm{tan}\beta =3`$. We write $`A_{ij}^{u,d}=m_0+\delta Ax_{ij}`$ and vary $`\delta A`$ from 0 to 3$`m_0`$ with $`x_{ij}`$ randomly selected in the range $`[1,1]`$. To calculate the mercury EDM we use the results of Ref. and the bounds on the mass insertions of Ref.. The neutron EDM requires Im$`(\delta _{LR}^{u,d})_{11}𝒪(10^6)`$, whereas the bound from the mercury EDM is Im$`(\delta _{LR}^{u,d})_{11}𝒪(10^7)`$.

We see that the EDMs exceed the experimental bounds by up to 3 orders of magnitude. In order sufficiently to suppress this effect, the sfermion masses must be pushed up above 10 TeV. This would however create a large hierarchy between the electroweak and the SUSY scales leading to a large fine-tuning , so one of the primary motivations for supersymmetry would be lost.

The effect described above seems quite close in spirit to the notorious FCNC problem. Yet, it is different from the latter since even if the FCNC are suppressed sufficiently, the EDM problem will still persist. Indeed, the strongest FCNC constraints come from the kaon physics and, for instance, $`ϵ_K`$ requires $`\sqrt{|\mathrm{Im}(\delta _{LR}^d)_{12}^2|}3.5\times 10^4`$ for the above parameters  . This bound is satisfied automatically in most SUGRA models with $`𝒪(1)`$ non-universality in the A-terms because these mass insertions are suppressed by the quark masses. In contrast, the level of sensitivity of the EDMs is $`𝒪(10^7)`$ and even a 10% nonuniversality can overproduce EDMs by an order of magnitude (Figs.1,2).

To be specific, let us consider an example. In the $`Z_6`$-I orbifold model of Refs., the Higgs fields belong to the $`\theta `$ twisted sector, quark doublets – to the $`\theta ^2`$ sector, and quark singlets – to the $`\theta ^3`$ sector (for the fixed point assignment, see ). The squark masses are generation–independent since different generations belong to the same twisted sector. The resulting A-terms depend on the Goldstino angle $`\mathrm{\Theta }`$ which measures the balance between the dilaton and moduli contributions to SUSY breaking. For $`|T|1`$ and a conservative value $`\mathrm{tan}\mathrm{\Theta }2`$ which corresponds to mostly dilaton SUSY breaking ($`F_S/F_{T_i}10`$), the arising A-term texture is

$`A^d\sqrt{3}m_{3/2}\left(\begin{array}{ccc}0& 0& 2\\ 1& 1& 3\\ 1& 2& 1\end{array}\right).`$ (11)
Here $`m_{3/2}`$ is the gravitino mass and we have assumed real $`F_S`$ and $`F_T`$. The resulting mass insertions are $`(\delta _{LR}^d)_{11}3\times 10^53\times 10^5`$i, $`(\delta _{LR}^d)_{12}2\times 10^5+2\times 10^5`$i, $`|(\delta _{LR}^d)_{23}|𝒪(10^4)`$, $`|(\delta _{LR}^d)_{13}|𝒪(10^5)`$. Thus, the SUSY contribution to $`K\overline{K}`$ mixing is negligible (only the contribution to $`ϵ^{}`$ is significant). The FCNC bounds involving the third generation are much looser and easily satisfied. On the other hand, the neutron and mercury EDMs are overproduced by one-two orders of magnitude. We found that these results are qualitatively stable under variation of the Goldstino angle. Clearly, this situation is quite general and holds in most models with $`𝒪(1)`$ non-universality in the A-terms.

The additional EDM contributions emerge due to the fact that string theory, as a fundamental theory, must explain the fermion mass patterns as well as accommodate a large CKM phase. The underlying mechanism generally affects the soft terms thereby inducing large EDMs. The two sectors, the Standard Model and the supersymmetric sector, which can be assumed to be disconnected in the context of the MSSM, are tightly related in string models and it is the absence of CP violation in one of them and its presence in the other that leads to the conflict.

In this context, EDMs serve as a probe of supersymmetric flavor physics as well as supersymmetry breaking. Non-observation of the EDMs implies that string theory selects special flavor structures or a special pattern of SUSY breaking. The non-universal contributions to the EDMs are absent if the A-terms happen to be $`universal`$. This requires that (1) the fields that generate the Yukawa matrices do not break SUSY,

$$F_T0,F_\varphi 0,$$
(12)
and (2) that there are no generation-dependent contributions from the Kähler potential in Eq.1. This occurs if different moduli are responsible for SUSY breaking and generating the Yukawa couplings, yet it remains to be understood why the latter moduli do not break supersymmetry. Exact universality can be a result of dilaton-dominated SUSY breaking. This scenario is known to lead to charge and color breaking minima which, however, can be avoided by lowering the string scale to intermediate values. On the other hand, dilaton domination would be hard to implement in semi-realistic string models which exhibit CP violation and dilaton stabilization ,.

The problem can also be avoided if the flavor structures are $`hermitian`$,

$$Y^a=Y^a,A^a=A^a.$$
(13)
In this case, the $`\widehat{A}`$ matrices are also hermitian and, since $`V_L^a=V_R^a`$, the hermiticity is preserved by the basis rotation . This eliminates the flavor-diagonal phases in the A-terms. The non-hermitian corrections appear only due to small renormalization group effects and lead to naturally suppressed EDMs. To justify the hermiticity requires the presence of an additional symmetry  and it is a non-trivial task to embed such models in string theory.

The problem becomes milder (yet it exists) if one assumes that the up Yukawa matrix is diagonal while the down one is proportional to the CKM matrix. In this case, the CP-phase appears only in the (13) and (31) entries of the down Yukawa matrix. This suppresses the effect of the super-CKM rotation, yet the mercury EDM often exceeds the limit by about an order of magnitude. Another interesting possibility is to have $`matrix`$-factorizable A-terms: $`\widehat{A}=AY`$ or $`YA`$. This suppresses the magnitude of the mass insertions since now there is a contribution from only one generation, $`(\delta _{LR})_{ii}m_i`$. However, again it is not clear whether these special cases can be obtained in realistic string models. In any case, the above considerations provide a selection rule: any realistic string model must satisfy the requirement of absence of the non-universal contributions to the EDMs.

To summarize, we have argued that string models which accommodate the fermion mass hierarchy and mixings generally lead to large EDMs even in the absence of CP-phases in the SUSY breaking auxiliary fields and the $`\mu `$-term. Non-observation of the EDMs implies that the supersymmetric structures have a special flavor or SUSY breaking pattern, which can be probed in current and future particle experiments. We have benefited from valuable discussions with G. Kane.