4 Run results

4.1 Generic disk response

Figure 1: This plot shows the disk surface density response to the planet potential, at the end of run R1744. The rotation is anti-clockwise. The one-armed shock triggered by the planet in the disk is clearly visible, as well as the shallow dip surrounding the orbit, the relative depth of which amounts to $\sim$15% of the disk surface density.

The disk response to the perturbing potential of a planet on a uniform circular orbit is illustrated at Fig. 1. The planet excites waves in the disk at all wavenumbers (Goldreich & Tremaine 1979), which corotate with the planet, and the superposition of which leads to a one-armed spiral shock (Goodman & Rafikov 2001). This latter is in advance w.r.t. the planet in the inner disk, which corresponds to the positive inner Lindblad torque, while it is delayed w.r.t. the planet in the outer disk, which corresponds to the negative outer Lindblad torque. The balance between the one-sided Lindblad torque and the viscous torque, for a steady flow in the planet frame, enables one to evaluate the depth of the dip which surrounds the orbit (Masset 2001), which, as long as it is shallow, scales proportionally to the planet mass square and inversely proportionally to the disk viscosity. The dynamics of the co-orbital region, i.e. the horseshoe region, is better seen in the $(\theta, r-r_{\rm p})$ plane, as shown in the following sections.

   
4.2 An illustrative example

Let one consider the runs R17₅¹ and R33₅¹, which correspond respectively to: $q=5\times 10^{-5}$, h=0.05 and $\hat\nu=1.93\times 10^{-7}$ and to: q=10^-4, h=0.05 and $\hat\nu=1.93\times 10^{-7}$. The aspect of the streamlines at time $t=754~\Omega_{\rm p}^{-1}$ in the frame corotating with the planet in either case is represented at Figs. 2 and 3.

Figure 2: Streamlines aspect in the co-orbital region for run R1751. The thick line represents the outer downstream separatrix.

Figure 3: Same as Fig. 3 for run R3350. The separatrix lies further from the planet orbit than in run R1751.

The streamlines are found integrating the velocity field, which is bilinearly interpolated from the values at the interfaces between adjacent zones. As most of the motion in the co-orbital region can be accounted for by a linear shear v_y=2Axwhere A is the Oort's constant, this bilinear interpolation gives satisfactory results, and in particular enables one to search for the precise position of the separatrix between the librating and circulating streamlines. The separatrix position is defined by its distance $x_{\rm s}$ to the orbit in opposition with the planet, i.e. at azimuth  $\theta =\pi$, as indicated in Figs. 2 and 3. The value of $x_{\rm s}$ is determined with a dichotomic search. If a guess of $x_{\rm s}$ gives a circulating (resp. librating) streamline, a smaller (resp. higher) value is tried. The search is iterated until a sufficient precision is achieved. Because the bilinear interpolation of the velocity field is well adapted to the kinematics of the co-orbital region, the precision that can be achieved on the separatrix position is much smaller than the radial zone size, and in that sense such a precision is illusory. However, this allows to investigate the behavior of the librating zone width when varying any parameter of the problem, the grid resolution being fixed. An estimate of the error on the corotation torque and on the separatrix position due to finite resolution effect can be found in Appendix A. For the run R17₅¹, the separatrix position is found to be:  $x_{\rm s}= 0.0424\pm10^{-4}$, while for the run R33₅¹, the separatrix position is found to be:  $x_{\rm s}= 0.0597\pm10^{-4}$, as can be estimated from Figs. 2 and 3. Therefore the viscous time-scale across the horseshoe zone half-width is:

$$\tau_{{\rm visc}}=\frac{x_{\rm s}^2}{3\nu} \simeq 3.1\times 10^3~\Omega_{\rm p}^{-1} \mbox{~for R}17_5^1$$ (16)

$$\simeq 6.2\times 10^3~\Omega_{\rm p}^{-1}\mbox{~for R}33_5^1$$ (17)

while the outermost horseshoe turnover time is:

$$\tau_{{\rm HS}}=\frac{8\pi}{3\hat x_{\rm s}}\Omega_{\rm p}^{-1} \simeq 200~\Omega_{\rm p}^{-1} \mbox{~for R}17_5^1$$ (18)

$$\simeq 140~\Omega_{\rm p}^{-1} \mbox{~for R}33_5^1.$$ (19)

As the viscous time-scale across the horseshoe region is much higher than the horseshoe turnover time in either case, one expects the corotation torque to saturate. The total torque exerted by the disk on the planet as a function of time is represented at Fig. 4.

Figure 4: Torque on the planet of run R1751 (solid line) and run R3050 (dashed line) as a function of time. The outermost horseshoe turnover time is indicated in either case. The torque is normalized by $\Gamma _0=\pi r_{\rm p}^4\Omega _{\rm p}^2\Sigma q^2/h^3$. The glitch around t=10 orbits corresponds to an edge effect of the smoothing window (these results are smoothed over a temporal window of about 10 orbits).

The total torque tends towards a constant value in time on a time-scale of the order of the outermost horseshoe turnover time. The initial torque value corresponds to the sum of the differential Lindblad and corotation torque, this later being positive, since the surface density is uniform on the co-orbital region. Incidently one can notice that in both cases the initial corotation torque value almost cancels out exactly the limit value at large time, and that this limit value at large time is not the same for both planets.

Figure 5: Torque on the planet of runs R175i as a function of time, for $1\leq i\leq 12$. The steady state regime is reached faster as the viscosity increases, and tends towards a partially saturated state function of the viscosity. The glitch at t=4 orbits is an edge effect of the smoothing window. For this series of runs only the integration time was 150 orbits in order to check that the total torque value does not vary significantly after 120 orbits.

Figure 5 shows the behavior of the total normalized torques as a function of time for the case $q=5\times 10^{-5}$, h=0.05, for different viscosities. The bottom curves correspond to low viscosity, and the limit value at large time increases with viscosity. This limit value reaches a maximum which corresponds roughly to the value at t=0, which corresponds to a totally unsaturated torque, and then decreases beyond the cut-off viscosity, as the topology of the co-orbital region undergoes a change that will be discussed further (see Masset 2001).

4.3 Dependence of corotation torque on the viscosity

Figure 6 shows the value of the normalized total torque acting on the protoplanet as a function of the reduced viscosity, evaluated from runs R17₅ⁱ. In each case the value of the torque is averaged on the time interval [90,120] orbits, in order to get rid of the transitory regime of the beginning, which corresponds, as we mentioned in Sect. 4.2, to the time needed for the corotation torque to (possibly partially) saturate. This figure also shows the best fits by the functional dependences given by Masset (2001), which are respectively:

$$\Gamma_{\rm C}=\Gamma_{\rm C}^{{\rm max}}\cdot \frac 14\left... ...\rm s}^3}-\frac{g(z_{\rm s})}{z_{\rm s}^4g'(z_{\rm s})}\right]$$ (20)

for the solid line,

$$\Gamma_{\rm C}=\Gamma_{\rm C}^{{\rm max}}\frac{g(z_{\rm s})}{z_{\rm s}g'(z_{\rm s})}$$ (21)

for the dashed line, and

$$\Gamma_{\rm C}=\Gamma_{\rm C}^{{\rm max}}\frac{g'(0)}{g'(z_{\rm s})}$$ (22)

for the dotted line, where $z_{\rm s}=\hat x_{\rm s}(2\pi\hat\nu)^{-1/3}$, where $\Gamma_{\rm C}^{\rm max}= (9/8)\Sigma_0\Omega_{\rm p}^2x_{\rm s}^4$, and where g(z) is a linear combination of the Airy functions Ai and Bi which cancels out for z=0, e.g.:

$$g(z)={\rm Bi}(z)-\sqrt 3{\rm Ai}(z).$$ (23)

Figure 6: Torque on the planet of runs R175i (diamonds). The value of the separatrix distance used for the fits is $x_{\rm s}=0.0388$. The vertical three-dot-dashed line shows the cut-off viscosity. Note that the fit maximum value at -0.01 is also the torque value at t=0 (cf. Fig. 4), and corresponds to a fully unsaturated corotation torque.

The fits presented at Fig. 6 have the following functional form:

$$\Gamma(z)=\Gamma_{{\rm min}}+\Gamma_{\rm C}(z)$$ (24)

where $\Gamma_{\rm C}$ is one of the three functions given by Eqs. (20), (21) or (22) and therefore has one free parameter, the floor level $\Gamma_{{\rm min}}$, which corresponds to the case of a saturated corotation torque and therefore to the sum of the differential Lindblad torque and viscosity independent additional term of Eq. (3), which scales with the one-sided Lindblad torque, and which turns out to be $\Gamma_{{\rm min}}=0.050$. Although $x_{\rm s}$ is not a fit free parameter, the strong dependence of $\Gamma_{\rm C}^{\rm max}$ upon $x_{\rm s}$ and the relatively high variation with azimuth of the separatrix distance to the orbit enables one to consider it as a free parameter within a narrow range (e.g. here in the range $0.038<\hat x_{\rm s}<0.043$, which can be gathered from examination of Fig. 2). The value used for the fit of Fig. 6 falls within this range of values, and it corresponds to the separatrix distance at $\theta\simeq 1.6$ rad. It should be noted that if one discards the lowest viscosity points for the fit (for which residual numerical viscosity effects may affect the results), then one can get a good fit by Eq. (20) as well, with a slightly larger value of $x_{\rm s}$ which still falls in the range $0.038<\hat x_{\rm s}<0.043$. Despite of the uncertainties introduced by the low viscosity points, the fit correctly accounts for a number of properties of the corotation torque main term (its absolute value, the half saturation for $\hat\nu\sim 0.05\hat x_{\rm s}^3$, and the fact that the partial saturation regime spans about two orders of magnitude of viscosity).

The behavior at high viscosities differs markedly from the fit by Eq. (24). Indeed a cut-off of the corotation torque is expected when the viscous drift time of a fluid element across the outer half of the horseshoe region is shorter than half of its libration time. In that case the region of fluid elements trapped in the co-orbital region has an azimuthal extent shorter than the azimuthal extent of the horseshoe region. Figures 7 and 8 represent the streamlines for two cases of high viscosity (resp. $\hat\nu=3.7\times 10^{-5}$ and $\hat\nu=5.2\times 10^{-4}$).

Figure 7: Streamline aspect for the run R1759, which corresponds to the peak value of the total torque of Fig. 6. The four thin solid lines represent librating streamlines. One is almost exactly closed, the other ones not quite, as a strictly steady state in the planet frame has not been reached yet. The set of trapped fluid elements still have an azimuthal extent which is just smaller than $2\pi$. The dashed line corresponds to a unique streamline, which originates in the outer disk and eventually reaches the inner disk. The thick solid line represents the outer downstream separatrix between circulating and librating streamlines. The shaded area shows the set of librating fluid elements, the boundary of which (the thick dot dashed line) is also the outer upstream separatrix and inner downstream separatrix.

Figure 8: Streamline aspect for the run R17513. The three thin solid lines represent librating streamlines. They are approximately closed. The azimuthal extent of the trapped material (shaded area) is only about 2 rad, and its boundary roughly coincides with the outer closed streamline. The dotted and dashed lines represent two open streamlines originating in the outer disk. The thick solid line represents the outer downstream separatrix.

The outer downstream separatrix of Fig. 7 roughly coincides with that of Fig. 2. The material lying interior to it is not librating however, as is illustrated by the dashed streamline. This material undergoes one horseshoe like close encounter with the planet and goes to the inner disk. The distance between the outer downstream separatrix and the boundary of the trapped region (given for $r>r_{\rm p}$ by the outer upstream separatrix) increases with the viscosity and eventually becomes as large as the horseshoe region. In that case, the material can go from the outer disk to the inner disk without any close encounter with the planet. Figure 8 shows the example of two open streamlines (dotted and dashed lines) which correspond to two fluid elements which reach the inner disk without any orbit crossing close encounter with the planet, a situation that is impossible at low viscosities. Note that the distance of the outer downstream separatrix to the orbit, measured for $\theta\simeq 0.75$ (instead of $\theta =\pi$ in order for the viscous drift not to affect sizably the measure at these high viscosities), is in the range [0.035; 0.04]in Figs. 2, 7 and 8. Therefore this value is almost constant over three decades of viscosity, which justifies a posteriori the choice of a unique value of $\hat x_{\rm s}$ for the fits of Fig. 6.

4.4 Other runs results

 Figures 9 to 12 show the value of the steady state normalized total torque as a function of viscosity for different planet masses, respectively for disk thickness h=0.03, h=0.04, h=0.05, and h=0.06.

Figure 9: Total torque exerted on the planet as a function of viscosity, for a disk aspect ratio h=0.03 and five different values of the planet mass.

It can checked on these figures that the viscosity at which the torque begins to decrease increases with the planet mass. This is in agreement with the fact that the cut-off at high viscosity occurs for $\hat\nu\simeq(1/4\pi)\hat x_{\rm s}^2$, and that $\hat x_{\rm s}$ increases with the planet mass. It can also be checked that the torque saturation occurs at lower viscosities for smaller masses, as the saturation occurs for $\hat \nu \sim \hat x_{\rm s}^3$.

Figure 10: Total torque exerted on the planet as a function of viscosity, for a disk aspect ratio h=0.04 and five different values of the planet mass.

Figure 11: Total torque exerted on the planet as a function of viscosity, for a disk aspect ratio h=0.05 and five different values of the planet mass.

Figure 12: Total torque exerted on the planet as a function of viscosity, for a disk aspect ratio h=0.06 and four different values of the planet mass.

 

These results can be used to try identify the additional corotation torque term of Eq. (3). The main problem is to evaluate the residual torque, i.e. to correct the measured torque from the main corotation torque term. This latter involves the separatrix distance to the orbit $x_{\rm s}$ to the fourth power, and it is therefore very sensitive to any error on the estimate of this separatrix distance. Two solutions can be adopted: the first one consists in considering that the measured torque value at low viscosity is the sum of the differential Lindblad torque and the additional corotation torque term (in which case there is no need to evaluate the main corotation torque term) and the second one consists in correcting the torque measured at large viscosity - although much before the cut-off - from this main term and attributing the residual to the sum of the differential Lindblad term and the additional corotation torque term. This latter solution turns out to be easier and more reliable, despite of the difficulty in measuring $x_{\rm s}$. Indeed the disk response at low viscosity may involve a strongly non-linear disk response (the dip around the orbit tends to become a fully qualified gap), which modifies the differential Lindblad torque in a way that is not trivial to correct. Furthermore, finite numerical viscosity effects (which are difficult to quantify) and the long viscous time-scale are two additional factors which may lead to a significant difference between the ideal low viscosity steady state situation and the run results after 120 orbits.

Figure 13: Separatrix shape for a low mass, low aspect ratio case. The distance of the separatrix to the orbit varies by about 23% over the azimuth and therefore its fourth power varies by a factor 2.3 with azimuth. $x_{\rm s}^+$ and $x_{\rm s}^-$ stand for the distance to the orbit of the outer and inner separatrices, respectively. They are functions of the azimuth $\theta$. In this example $\Delta \theta =1.5$ rad.

Evaluating the corotation torque main term involves a precise measurement of the value of $x_{\rm s}$. It can be seen in Fig. 13 how the distance between the separatrix and the orbit varies with azimuth. The value for $x_{\rm s}$ is chosen to be the arithmetical average of the outer downstream separatrix position at a given angle $\Delta\theta$ from outer conjunction and the inner downstream separatrix position at the corresponding angle $2\pi-\Delta\theta$ from inner conjunction:

$$x_{\rm s}=\frac 12[x_{\rm s}^+(\Delta\theta)+x_{\rm s}^-(2\pi-\Delta\theta)].$$ (25)

The value of $\Delta\theta$ is then varied between 1 and $\pi$ rad, which leads to a dispersion in the corotation torque main term evaluation, and this leads therefore to a dispersion of the residual torque estimate. The corresponding residual torque estimates are shown for the four different disk thicknesses in Figs. 14 to 17.

Figure 14: Residual normalized torque as a function of the separatrix distance, for the 3% aspect ratio disk. The dotted line shows the linear regression fit. The error bars are obtained by taking respectively the maximum and minimum value of the average of $x_{\rm s}^+(\Delta \theta )$ and $x_{\rm s}^-(2\pi -\Delta \theta )$, when $\Delta$ varies.

Figure 15: Residual normalized torque as a function of the separatrix distance, for the 4% aspect ratio disk. The dotted line shows the linear regression fit.

Figure 16: Residual normalized torque as a function of the separatrix distance, for the 5% aspect ratio disk. The dotted line shows the linear regression fit. The white square shows a measurement for which the separatrix distance to the orbit is smaller than the radial zone width, which means that the librating region is narrower than two zones in radius. It is therefore discarded for the linear regression fit.

Figure 17: Residual normalized torque as a function of the separatrix distance, for the 6% aspect ratio disk. The dotted line shows the linear regression fit.

The results of the linear regression fits performed on these figures are shown in Table 3. A number of comments can be made from these figures and from the table:

1.

The correction is performed at the peak value of the total term, while the separatrix position is searched for at a slightly lower viscosity. The correction takes into account a possibly partial saturation of the corotation torque at its peak value.

2.

The function which is used to correct the main corotation term from partial saturation is the function given by Eq. (20). Inspection of Fig. 6 shows that using either Eq. (21) or Eq. (22) would not change significantly the results, given the uncertainty on $x_{\rm s}$, and given the fact that at its peak value the corotation torque is almost fully unsaturated.

3.

The value at $x_{\rm s}=0$ of the linear regression fit corresponds roughly to the limit value that can be deduced by eye from the limit value at $\hat \nu=0$ in Figs. 9 and 10. In the case of h=0.05, the uncertainties are considerable, and in the h=0.06 case, no curve reaches the fully saturated regime at low $\hat\nu$.

4.

The mean slope $\overline {\cal G}$ distribution is compatible with the expectation. It should be mentioned however that the function ${\cal G}(x_{\rm s})$ is not expected to be a constant, and that its value, which is such that: $1\leq{\cal G}(x_{\rm s})\leq7/3$, depends on the location $x_{\rm s}$ at which the separatrix samples the axisymmetric profiles perturbations.

5.

The variation with mass of the peak value of the total normalized torque for a given aspect ratio can be accounted for by the slope of the fit, i.e. this variation is roughly equal to $\overline{\cal G}(\hat x_{\rm s}^{\rm max}-\hat x_{\rm s}^{\rm min})\Gamma_{\rm LR}$. This therefore accounts for the increasing peak value with the planet mass, whereas one would expect a constant maximum in the linear regime ( $x_{\rm s}^4\propto q^2$) and a decreasing peak value beyond ( $x_{\rm s}^4 \propto q^{4/3}$).

6.

Not all the behavior observed is entirely due to the additional corotation torque term. Indeed there is no reason to assume that the normalized differential Lindblad torque is a constant for a given aspect ratio. Although this is certainly true in the linear regime, it is no longer true for higher masses. Miyoshi et al. 1999 note that in an infinitesimally thin disk, in the non-linear regime (the threshold of which they find at $R_{\rm H} = H/2$, where $R_{\rm H}$ is the planet Hill radius), the one-sided Lindblad torque is smaller than its linear estimate. They mention that part of the cut-off comes from material inside the horseshoe region, which exerts an opposite torque on the perturber, and weakens the linear estimate of the Lindblad torque. As the material in the horseshoe region is librating in the planet frame, this material in average corotates with the planet, and therefore participates in the corotation torque and no in the Lindblad torque. If one defines the Lindblad torque as the torque exerted by the circulating material on the planet, which is, in the steady state regime, the only definition which ensures that the total torque exerted on the planet is the sum of the corotation and Lindblad torques (as the co-orbital dynamics partitions the disk into librating and circulating fluid elements) then it is not clear whether the behavior observed by Miyoshi et al. is due to a Lindblad torque non-linear cut-off or to an avatar of the additional corotation torque term. If one performs however a local scattering calculation (Lin & Papaloizou 1979; Papaloizou & Lin 1984) over the circulating fluid elements, then one expects a q^-1 cut-off of the one-sided Lindblad torque when the separatrix lies further than (2/3)H from the orbit, i.e. when the horseshoe region invades what would be otherwise the place of disk-perturber angular momentum exchange through Lindblad torques. It should be noted that on the present data set, the conditions $x_{\rm s} > (2/3)H$ and $R_{\rm H} >(1/2)H$ are almost strictly equivalent. In the h=3% case (Fig. 14), all points but the bottom one correspond to situation beyond the non-linear threshold, while in the h=6% case (Fig. 17), all points correspond to a situation where the one-sided Lindblad torque is still in the linear regime (the Hill radius of the most massive planet is precisely $0.03r_{\rm p}$). Therefore the behavior observed for the thinner disks is likely to be contaminated by a differential Lindblad torque non-linear cut-off (which also conspires to lift the normalized curves as the planet mass increases), whereas the behavior observed for the h=0.06 disk is exempt of this contamination and should be due entirely to the additional corotation torque term of Eq. (3).

 

 

Table 3: Linear fit regression values for Figs. 14 to 17. The slope $\overline {\cal G}$ is defined by $\Gamma(\hat x_{\rm s})=\Gamma_{\hat x_{\rm s}\rightarrow 0} + \overline{\cal G}\hat x_{\rm s}\Gamma_{\rm LR}$, where the normalized value of $\Gamma _LR$ for all the disks is 0.21.

${h (\%)}$	slope $\overline{\cal G}{\phantom{\int^W}}$	$\Gamma_{\hat x_{\rm s}\rightarrow 0}$	correlation
3	1.61	-0.026	0.978
4	1.61	-0.046	0.90
5	1.89	-0.060	0.44
6	2.15	-0.077	0.997

Despite of the difficulty of these measures, these results confirm the existence of a term in the total torque which scales roughly as $\hat x_{\rm s}\Gamma_{\rm LR}$, in agreement with the expectation that ${\cal G}(x_{\rm s})=O(1)$, and ${\cal G}(x_{\rm s}) < 7/3$.