3 Numerical aspects

The code that was used is an Eulerian ZEUS-like 2D code based on a polar grid centered on the primary (Stone & Norman 1992; Nelson et al. 2000), and corotating with the planet (Kley 1998). The runs described here were performed using a modified azimuthal Courant condition (Masset 2000) in order to increase the timestep, and some of these runs were checked against the standard azimuthal transport procedure, and found to give almost identical results. The grid consists of $N_\theta=450$ by N_r=143 zones, uniformly spaced in azimuth and radius. The planet lies at radius $r_{\rm p}=1$ and azimuth $\theta=0$. The central star mass and gravitational constant are respectively M_*=1 and G=1, and the time unit is chosen to be $(r_{\rm p}/GM_*)^{3/2}=\Omega_{\rm p}^{-1}$, so that the planet orbital period in our unit system is $2\pi$. The grid outer boundary is chosen to be at r=2.5, while the inner boundary is at r=0.504651, so that the planet is located just at the center of a zone, in order to limit any bias in the torque and in the co-orbital dynamics due to an uneven placement of the planet w.r.t. the grid.

3.1 Basic equations

The equations solved by the code are:

• The continuity equation, which reads:
  

  $$\frac{\partial\Sigma}{\partial t}+\frac 1r\frac{\partial(rv_r... ... r}+\frac 1r\frac{\partial(v_\theta\Sigma)}{\partial\theta}=0.$$ (4)

  
  

• The Navier-Stokes equations, which reads respectively for v_r and $v_\theta$:
  

  $$\frac{\partial v_r}{\partial t} +v_r\frac{\partial v_r}{\part... ...\partial r} -\frac{\partial\phi}{\partial r} +\frac{f_r}\Sigma$$ (5)

  
  
  and:
  

  $$\frac{\partial v_\theta}{\partial t} +v_r\frac{\partial v_\th... ... 1r\frac{\partial\phi}{\partial \theta} +\frac{f_\theta}\Sigma$$ (6)

  
  
  where p is the vertically integrated pressure, f_r and $f_\theta$ respectively the radial and azimuthal component of the viscous force per unit surface, and $\phi$ is the gravitational potential.
• The equation of state is that of an isothermal gas with sound speed $c_{\rm s}$:
  

  $$p=c_{\rm s}^2\Sigma.$$ (7)

  
  

• The gravitational potential includes the potential of the primary, the potential of the protoplanet, and an indirect term which arises from the fact that the non-rotating frame centered on the primary is not inertial:
  
   

  $$\phi(r,\theta) -\frac{GM_*}{r}-\frac{Gm_{\rm p}}{(r^2+r_{\rm p}^2-2rr_{\rm p}\cos\theta+ \epsilon^2)^{1/2}}$$

  $$+\frac{Gm_{\rm p}}{r_{\rm p}^2}{r\cos\theta}$$

  $$+r\int_{{\rm disk}}\frac{G\Sigma(r',\theta')}{r'^2} \cos(\theta-\theta')r'{\rm d}r'{\rm d}\theta'.$$ (8)

  
  
  In this expression, $\epsilon$ is the smoothing length of the protoplanet potential, and otherwise stated it amounts to 60% of the disk "thickness'' $H=c_{\rm s}/\Omega$.
• The viscous force is derived from the viscous stress tensor as follows:
  
    

  $$\frac{1}{r} \frac{\partial (r \tau_{rr})}{\partial r} + \frac{1}{r} \frac{\partial \tau_{r \phi}}{\partial \phi} - \frac{\tau_{\phi \phi}}{r}$$ (9)

  $$f_{\phi} \frac{1}{r} \frac{\partial (r \tau_{\phi r})}{\partial r} + \frac{1}{r} \frac{\partial \tau_{\phi \phi}}{\partial \phi} + \frac{\tau_{r \phi}}{r},$$ (10)

  
  
  where the components of the viscous stress tensor are:
  
   

  $$\tau_{rr} 2 \eta D_{rr} - \frac{2}{3} \eta \nabla . \vec{v}$$ (11)

  $$\tau_{\phi \phi} 2 \eta D_{\phi \phi} - \frac{2}{3} \eta \nabla . \vec{v}$$

  $$\tau_{r \phi} \tau_{\phi r} = 2 \eta D_{r \phi},$$

  
  
  where
  
   

  $$\frac{\partial v_r}{\partial r}, D_{\phi \phi} = \frac{1}{r} \frac{\partial v_{\phi}}{\partial \phi} + \frac{v_r}{r}$$ (12)

  $$D_{r \phi} \frac{1}{2} \left[ r \frac{\partial}{\partial r} \left( \frac{v_{\phi}}{r} \right) + \frac{1}{r} \frac{\partial v_r}{\partial \phi} \right],$$

  
  
  and $\eta=\Sigma \nu$ is the vertically integrated dynamical viscosity coefficient.

3.2 Boundary conditions

The boundary conditions are $2\pi$-periodic in $\theta$ in order to account for the disk geometry. As the viscosity can be large in some runs, it is important to take adequate boundary conditions, otherwise the fast radial redistribution of disk material can affect the slope of the surface density profile, which in turn can significantly alter the corotation torque magnitude. For this reason, it is important :

• to have a source of material with surface density $\Sigma_0$ at large radius in order to avoid outer disk depletion;
• to allow the outflow (inwards) of disk material at the inner boundary only if the disk surface density is larger or equal to $\Sigma_0$ (otherwise a positive radial gradient of surface density develops, and one may overestimate the corotation torque in these conditions).

3.3 Initial setup

In all the runs presented here, the unperturbed disk surface density is uniform and is, in our unit system, $\Sigma_0=10^{-4}$. This translates into:

$$\Sigma=\Sigma_0\cdot 8.9\times 10^6\cdot\left(\frac{M_*}{M_\o... ...ft(\frac{r_{\rm p}}{1\mbox{~AU}}\right)^{-2}\mbox{~g~cm}^{-2}.$$ (13)

Note that since the results presented here are scaled by $\Gamma_0$, they are independent of the actual value of $\Sigma_0$ (except for the disk potential indirect term of Eq. (8), which for the value chosen above is negligible).

The disk temperature profile is constant in time and is chosen such that the disk aspect ratio h=H/r be uniform. The initial azimuthal velocity is computed accordingly and is slightly sub-Keplerian due to the partial support of gravity by the radial pressure gradient:

$$v_\theta=\left[\frac{GM_*}{r}(1-h^2)\right]^{1/2}.$$ (14)

The initial radial velocity is set to zero.

3.4 Initial parameters

A number of runs have been performed varying the planet mass, the disk aspect ratio, the viscosity, and in some cases the smoothing prescription or the resolution (which was then chosen twice higher). Each run consisted of 120 planet orbits (which was assumed to be sufficient to reach a steady state in the planet frame). For a given choice of the planet mass and disk aspect ratio, usually 13 runs were performed for different values of the viscosity, logarithmicly spaced:

$$\hat\nu_i=10^{-7+2i/7}\mbox{~~~~~}(1\leq i\leq13).$$ (15)

The set of main runs is described at Table 1. The last cell of the h=6% line is empty: the corresponding sets of runs have not been performed, as the corresponding run would correspond to a horseshoe region radially not wider than one zone, and even if one disregards finite resolution effects, the mismatch between corotation and orbit is not negligible compared to the horseshoe zone half-width, which is one of the condition of validity of Eqs. (2) and (3). The first cell of the h=3% line is also empty, as the situation it would describe would rather correspond to type II migration, i.e. the planet opens a gap around its orbit and the corotation torque is switched off, except for the largest viscosities.

 

 

Table 1: The set of 247 main runs. A run with mass ratio q, aspect ratio hand viscosity $\hat\nu_i$ is labeled Rm_hⁱ, where h is in percent and where m is the round-off of the planet mass, in earth masses, if $M_*=M_\odot$. For instance, run R6₃⁸ stands for the run with $q=1.67\times 10^{-5}$, aspect ratio h=0.03, and viscosity $\hat\nu_8$. The bulk of the runs are Rm_hⁱ, with m=2,6,11,17,30, h=3,4,5,6% and i=1 to 13. Runs with m=30 have been performed for extrapolation purposes only, since for such a mass a protoplanet is expected to rapidly accrete the surrounding gas.

$_q~ ^{h (\%)}$	3	4	5	6
$6.67\times 10^{-6}$	$\nu_{1\rightarrow 13}$	$\nu_{1\rightarrow 13}$	$\nu_{1\rightarrow 13}$	n/a
$1.67\times 10^{-5}$	$\nu_{1\rightarrow 13}$	$\nu_{1\rightarrow 13}$	$\nu_{1\rightarrow 13}$	$\nu_{1\rightarrow 13}$
$3.33\times 10^{-5}$	$\nu_{1\rightarrow 13}$	$\nu_{1\rightarrow 13}$	$\nu_{1\rightarrow 13}$	$\nu_{1\rightarrow 13}$
$5\times 10^{-5}$	$\nu_{1\rightarrow 13}$	$\nu_{1\rightarrow 13}$	$\nu_{1\rightarrow 13}$	$\nu_{1\rightarrow 13}$
10-4	n/a	$\nu_{1\rightarrow 13}$	$\nu_{1\rightarrow 13}$	$\nu_{1\rightarrow 13}$

Subsidiary runs have been performed in which either the smoothing or the resolution has been changed with respect to the main runs. These runs are presented at Table 2.

 

 

Table 2: The 26 subsidiary runs. The smoothing length of the potential has been set to a smaller value (resp. 20% and 40% of the disk thickness) than its usual value (60% of the disk thickness).

Run name	corresponding main run	parameter changed
S20R155i	R155i, $i=1\rightarrow 13$	Smoothing:  $\epsilon=0.2H$
S40R155i	R155i, $i=1\rightarrow 13$	Smoothing:  $\epsilon=0.4H$

Subsidiary runs have been performed in which either the smoothing or the resolution has been changed with respect to the main runs. The runs with a modified smoothing length are presented at Table 2, whereas the high resolution runs are presented at Sect. 6.