Appendix A: Finite resolution effects on the co-orbital corotation torque estimate

As the radial width of the co-orbital region in the simulations presented in this work amounts to a small number of zone widths (typically 3 to 18), it is of interest to investigate how bad or good is the code used here in giving an estimate of the co-orbital corotation torque. For this purpose, one can consider a simplified situation, such as the one depicted at Fig. A.1. This situation neglects the grid curvature and describes the gas flow around the planet in a shearing sheet. In this situation, the corotation torque is given by

$$\Gamma_{\rm ss}=\frac 34\Omega^2x_{\rm s}^4r_{\rm p}\frac{\partial\Sigma}{\partial r}\cdot$$ (A.1)

Figure A.1: Sketch of a shearing sheet case. A radial set of zones is represented. The separatrices are assumed to be horizontal and are symmetric w.r.t. the orbit. The azimuthal velocity is set on a staggered mesh (centered in radius, staggered in azimuth), and is assumed in this simplified example to be independent of azimuth.

The code used here enforces angular momentum conservation (to the computer accuracy) hence the balance of angular momentum gained/lost by the orbit crossing fluid elements on horseshoe streamlines can be used as an estimate of the corotation torque. The angular momentum is a zone centered variable. In this simplified situation the velocity field is assumed to be the unperturbed one. The zone radial index is denoted i, and vanishes at the orbit (as stated in Sect. 3, the planet lies in the center of a zone). The index of the outermost zone fully embedded in the co-orbital region is denoted $i_{\rm M}$. It is assumed that the co-orbital region is wide enough so that such a zone exists. The expression for $i_{\rm M}$ is therefore:

$$i_{\rm M}=E\left(\frac{x_{\rm s}}{\Delta}-\frac 12\right)$$ (A.2)

where E(X) stands for the integer part of X, and where $\Delta$ is the radial zone width. The co-orbital corotation torque in this situation can therefore be expressed as:

$$\Gamma^{\rm {num}}_{\rm ss}=\Gamma^{\rm {num}^+}_{\rm ss}+\Gamma^{\rm {num}^-}_{\rm ss}$$ (A.3)

where:

 

$$\Gamma^{\rm {num}^+}_{\rm ss} \sum_{i=0}^{i_{\rm M}}\frac 32\Omega_{\rm p} i\Delta^2\cdot \Omega_{\rm p}r_{\rm p}\Delta i\Sigma_i$$ = (A.4)

$$+\left(\frac{x_{\rm s}}{\Delta}-i_{\rm M}-\frac 12\right)\frac 32\Omega_{\rm p}^2 (i_{\rm M}+1)^2\Delta^3r_{\rm p}\Sigma_{i_{\rm M}+1}$$

is the torque due to the horseshoe fluid elements flowing from $i\Delta>0$ to $-i\Delta$, and where:

 

$$\Gamma^{\rm {num}^-}_{\rm ss} -\sum_{i=0}^{-i_{\rm M}}\frac 32\Omega_{\rm p} i\Delta^2\cdot \Omega_{\rm p}r_{\rm p}\Delta i\Sigma_i$$ = (A.5)

$$-\left(\frac{x_{\rm s}}{\Delta}-i_{\rm M}-\frac 12\right)\frac 32\Omega_{\rm p}^2 (i_{\rm M}+1)^2\Delta^3r_{\rm p}\Sigma_{-i_{\rm M}-1}$$

is the torque due to the horseshoe fluid elements flowing from $i\Delta<0$ to $-i\Delta$. If one writes:

$$\Sigma_i=\Sigma_0+i\Delta\frac{\partial\Sigma}{\partial r}+O[(i\Delta/r_{\rm p})^2]$$ (A.6)

then combining Eqs. (A.4) and (A.5) one is led to:

 

$$\Gamma^{\rm {num}}_{\rm ss} \simeq 3\Omega_{\rm p}^2\Delta^3r_{\rm p}\frac{\partial\Sigma}{\partial r}$$

$$\times\left[\sum_{i=0}^{i_{\rm M}}i^3+(i_{\rm M}+1)^3\left(\frac{x_{\rm s}}{\Delta}-i_{\rm M}-\frac 12\right)\right]\cdot$$ (A.7)

The relative error on the corotation torque evaluation, ${\cal E}=(\Gamma^{\rm {num}}_{\rm ss}-\Gamma_{\rm ss})/\Gamma_{\rm ss}$, is therefore given by:

$${\cal E}(Y) = \frac{4}{Y^4}\left[\sum_{i=0}^{i_{\rm M}}i^3+(i_{\rm M}+1)^3\left(Y-i_{\rm M}-\frac 12\right)\right]$$ (A.8)

where $Y=x_{\rm s}/\Delta$ is the horseshoe region half-width, expressed in zone radial widths. The graph of the function ${\cal E}(Y)$ is the solid line of Fig. A.3.

The second step of this error estimate consists of evaluating the discrepancy between the actual and effective separatrix position. Indeed, the brackets of the second line of Eqs. (A.4) and (A.5) stand for the mass fraction of the outermost horseshoe zone included within the separatrix (i.e. it is also the ratio of the shaded area surface to the zone surface in Fig. A.1). There is no warranty however that the separatrix position which has to be used in this evaluation coincides with the one provided by the streamline analysis. In other words, one has to evaluate how much separatrix crossing numerical effects can lead to, in order to infer how distant the actual and effective separatrices can be. One approximate way to do that is to consider a 1D problem (along the radial dimension). Disregarding azimuthal advection should not lead to significant errors, as this latter is treated apart during the advection procedure based on the operator splitting technique, and as azimuthal motion should not lead to separatrix crossing provided this latter is horizontal (or weakly tilted) in the $(\theta, r-r_{\rm p})$ plane.

Figure A.2: Sketch of the 1D mesh used to evaluate the impact of numerical effects on separatrix crossing. The linear density $\Sigma$ is zone-centered, while the radial velocity u is defined at zone interfaces (on a staggered mesh).

Following the notations of Fig. A.2, one can write the total mass M_- contained below the separatrix (or below any arbitrary position $x_{\rm s}$) at time t as:

$$M_- = \Delta\sum_{j=0}^{i-1}\Sigma_j+\Sigma_i(x_{\rm s}-i\Delta)$$ (A.9)

where $\Delta$ is the zone spacing. This mass will be referred to hereafter as the inner mass. At time $t+\delta t$, the inner mass becomes:

$$\Delta\sum_{j=0}^{i-1}\Sigma_j'+\Sigma_i'(x_{\rm s}'-i\Delta)$$ (A.10)

where the primes denote the new quantities, and where one has the following relationship:

$$\Sigma_j'=\Sigma_j+\frac{\delta t}{\Delta}(F_j-F_{j+1})$$ (A.11)

where F_j is the mass flux (oriented upwards) across the interface between the zones j-1 and j. One also needs the following relationship:

$$x_{\rm s}'=x_{\rm s}+\delta t\left[u_i\left(i+1-\frac{x_{\rm ... ...}\right)+u_{i+1}\left(\frac{x_{\rm s}}{\Delta}-i\right)\right]$$ (A.12)

which comes from the fact that the streamlines (and therefore the separatrix) are found integrating the linearly interpolated velocity field. One can therefore write:

$$\Delta\sum_{j=0}^{i-1}\Sigma_j+\delta t\sum_{j=0}^{i-1}(F_j-F_{j+1})$$ M_-' =

$$+\left[\Sigma_i+\frac{\delta t}{\Delta}(F_i-F_{i+1})\right](x_{\rm s}'-i\Delta).$$ (A.13)

The mass fluxes can be written as:

$$F_j=\Sigma_j^*u_j$$ (A.14)

where $\Sigma_j^*$ is an estimate of the linear density at the interface, the expression of which depends on the numerical method (Stone & Norman 1992). To lowest order in $u_j\delta t/\Delta$, the mass variation can be expressed, after a few transformations, and assuming F₀=0, as:

$$\delta M_- \delta t\left[u_i(\Sigma_i-\Sigma_i^*)\left(i+1-\frac{x_{\rm s}}{\Delta}\right)\right.$$ =

$$\left.+u_{i+1}(\Sigma_i-\Sigma^*_{i+1})\left(\frac{x_{\rm s}}{\Delta}-i\right)\right].$$ (A.15)

Figure A.3: Relative error on the co-orbital corotation torque in the shearing sheet, as a function of resolution (number of zones in the horseshoe region half-width). The dotted and dashed lines are respectively for the worst case and a conservative case for the runs presented here (see text for details).

This mass variation necessarily corresponds to the separatrix crossing mass, since F₀=0. One can get an estimate of this mass variation assuming $u_i=u_{i+1}\equiv u$ (a reasonable assumption in the case of a radial viscous drift) and writing $\Sigma_i-\Sigma^*_i=\Sigma_{i+1}^*-\Sigma_i\simeq\partial_r\Sigma\Delta/2$ (in which case one has to assume that $\vert u\delta t\vert\ll\Delta$, a reasonable assumption for the case of interest here). One is then led to:

$$\delta M_-\simeq u~\delta t~\partial_r\Sigma(x_i^c-x_{\rm s})$$ (A.16)

where $x_i^c=(i+1/2)\Delta$ is the x value of the zone center. One can therefore write:

$$\frac{{\rm d}M_-}{{\rm d}t} = u(x_i^c-x_{\rm s})\partial_r\Sigma = \dot{x_{\rm s}}(x_i^c-x_{\rm s})\partial_r\Sigma$$ (A.17)

from which one can infer:

$$M_-=M_-^{\rm c}-\frac 12(x_{\rm c}^i-x_{\rm s})^2\partial_r\Sigma$$ (A.18)

where $M_-^{\rm c}$ is the inner mass when the separatrix crosses the zone center. Equation (A.18) shows that the inner mass at the top ( $x_{\rm s}=x_{\rm c}^i+\Delta/2$) and the bottom ( $x_{\rm s}=x_{\rm c}^i-\Delta/2$) of the zone are the same, from which one can conclude that, although some mass can cross the separatrix as its sweeps the zone, there is no net effect as it sweeps the entire zone from one interface to the other. As a consequence, the maximum amount of mass that can cross the separatrix as it sweeps an arbitrary number of zones is $M_{\rm sc}=\partial_r\Sigma\cdot\Delta^2/8$. Denoting with $\delta x_{\rm s}$ the error on the separatrix position, which is defined as: $\delta x_{\rm s} \cdot\Sigma=M_{\rm sc}$, one is led to the following relation:

$$\delta x_{\rm s} = \frac{\Delta^2}{8}\frac{\partial_r\Sigma}{\Sigma}\cdot$$ (A.19)

The worst case is obtained for $\partial_r\Sigma\sim\Sigma/\Delta$, i.e. when the separatrix lies on the edge of a fully emptied gap, and when this edge is not resolved (i.e. when the density falls abruptly from $\Sigma$ to 0 from one zone to its neighbor). It should be stressed that this situation is far from being met in the simulations presented here. Even in this worst case, the maximal error on the separatrix position will be: $\delta x_{\rm s}=\Delta/8$, whereas in a more reasonable case for which $\partial_r\Sigma$will typically amount to a small fraction of $\Sigma/H$, the maximal error on the separatrix position will be a fraction of $(\Delta/H)\cdot \Delta/8$. The worst case is displayed in Fig. A.3 with dotted lines, while the conservative case where $\delta x_{\rm s}=\Delta /16$, obtained for $\partial_r\Sigma=\Sigma/H$ and H=0.03 in the normal resolution runs, is displayed with dashed lines. It is worth noting that the separatrix error is much smaller than this conservative estimate, in particular for small planet masses, which hardly perturb the disk density profile, and for which $Y=x_{\rm s}/\Delta$ is the smallest. One can conclude from Fig. A.3 that the corotation torque is described with an accuracy better than 15% as soon as $x_{\rm s}/\Delta\ga 2.3$. Although this condition is not met for the smallest mass planets in the normal resolution runs presented here, the corotation torque, which on average for $1\leq Y\leq 2$ numerical effects tend to overestimate, is never found to be large enough to cancel out the differential Lindblad torque for these objects. Therefore, one can conclude that even with a relatively small number of mesh rings involved in the co-orbital region, the code used here can describe with a reasonable precision the dynamics of the co-orbital corotation torque, and that whenever the planet mass is too small for the dynamics of this region to be correctly described (for the resolution and parameter set used here), the corotation torque is too small to significantly affect the migration anyway (one recovers the linear regime, see Sect. 1 for details).