We show in this section the validity of the modified transport algorithm when applied to the interaction of a Jupiter sized protoplanet with a minimum mass protoplanetary disk in which it is embedded. The perturbed potential associated with the planet excites spiral density waves in the disk, which propagate away both inwards and outwards, with a pattern frequency equal to the planet orbital frequency. The spiral waves interact with the disk and give it the angular momentum they removed from the planet, and eventually open a gap centered on the planet orbit, provided the planet mass is high enough (Papaloizou & Lin 1984). We present a run with a one solar mass primary, one Jupiter mass protoplanet initially on a fixed circular orbit at r₀=5 AU embedded in a standard protoplanetary nebula whose parameters have been mentioned above. The grid has an inner radius at 2 AU and an outer radius at 12.5 AU. The sequence $(R_i)_{i\in[0,N_{\rm r}]}$ is equally spaced, with $N_{\rm r}=49$ ; The grid has $N_{\rm s}=143$ sectors, it is fixed in a non-Galilean non-rotating frame centered on the primary. Its outer boundary is rigid and its inner boundary allows outflow but no inflow. The disk aspect ratio is set to $4\ 10^{-2}$ everywhere. The planet perturbed potential is smoothed on a length scale which amounts to 40% of the Roche radius. In Eq. (14) we choose C₀=0.5. We plot in Fig. 2 the quantity $e_{ij}=\frac{v^{\theta a}_{ij}\Delta t}{\Delta y_i}$ after 2.86 orbits. This quantity represents the effective CFL ratio. With the standard transport algorithm this ratio is bounded by C₀.

$\begin{figure}\par\psfig{file=elong.ps,width=8.8cm} \par\end{figure}$

Figure 2: Number of cells crossed during one timestep. See text for parameters. The inner boundary is at the left (high values) and the outer boundary at the right (low values)

We see that the innermost ring sweeps almost four cells on one timestep, hence the use of the FARGO transport algorithm in this case results in a speed-up by a factor $\sim 8$ of the computation. One can note that the difference in e_ij between the innermost ring and its immediate neighbor is 0.5, which is the maximum allowed by Eq. (14). Indeed the timestep in this run is shear-limited, and the constraint on the residual velocities only would lead to an even bigger timestep, since as one can see the residuals of the distance swept over one timestep amounts to far less than 1/2, even in the vicinity of the planet. Indeed, runs performed with a logarithmic polar grid (i.e. with R_i+1/R_i constant), which have a smaller value of R_i+1-R_iin the inner part, have shown to allow a speed-up by a factor $\sim 30$ wrt the standard method.

In order to see how numerical viscosity affects the disk response in both cases, we plot in the Fig. 3 the disk density after 28.6 orbits, obtained from different algorithms. The left plot corresponds to a non-rotating frame standard transport run, while the middle plot represents a non-rotating frame FARGO transport run, and the right plot represents a standard transport run in a frame corotating with the planet (hence the planet is fixed with respect to the grid, so we expect from the results of Sect. 4 the density response in the vicinity of the planet to be given with a high accuracy). Note that special care has to be devoted to the treatment of the Coriolis force in that case in order to conserve exactly the angular momentum and then to avoid a spurious outwards transport in the disk (Kley 1998). We clearly see that the global spiral pattern excited by the protoplanet in the disk is identical in the three cases, though the response in the immediate vicinity of the planet is much more spread out in the non-rotating frame standard accretion case (left plot), and that the most sharply peaked response is achieved through the use of a corotating frame (right plot), as expected. Indeed, we plot in Fig. 5 a cut of the disk density at the planet radius in the three cases. The solid line represents the FARGO transport result, and the dot-dashed line the corotating frame result. They both have the same width, though the maximum of the density in the corotating case is higher. The dashed curve represents the result of the standard transport in a non-rotating frame. Its width is about twice as large as the other curves' width, and we also see that numerical effects in that case lead to additional leading and trailing material (near cells number 65 and 77), and to a smaller density peak value.

$\begin{figure}\par\psfig{file=threedens.ps,width=13.5cm} \par\end{figure}$

Figure 3: Disk density $\Sigma _{ij}$ ; j is in abscissa and i in ordinate. The left plot has been obtained by a non-rotating frame standard method, the middle one by a non-rotating frame FARGO transport method and the right one by a corotating frame standard method. Since each of these plots is approximately square, any circular feature in the disk should appear on the plots as a 1:3 vertical ellipse. This is not quite the case of the material surrounding the planet in the left panel, which leads to the conclusion that in a non-rotating frame standard transport method, the matter is artificially elongated along the orbital motion. The FARGO case, in the middle panel, shows much better behavior, and the coorbital material has a distribution which looks very much like the right panel one

$\begin{figure}\par\psfig{file=logplots.ps,width=13.5cm} \par\end{figure}$

Figure 4: Disk density $\Sigma _{ij}$ ; j is in abscissa and i in ordinate, for the log-grid runs described in the text. The left plot has been obtained by a non-rotating frame standard method, the middle one by a non-rotating frame FARGO transport method and the right one by a corotating frame standard method. The same comments as in Fig. 3 apply here. On this specific example, the FARGO run turned out to be 17 times faster than the standard run in the non-rotating frame, and 15 times faster than the standard run in the co-rotating frame

The FARGO plot in Fig. 3 exhibits at its inner boundary an oscillatory behavior which originates from three combined effects. First, this is a shear-limited run -- see Eq. (14) and Fig. 2 --; if we change the 0.5 factor in Eq. (14) to 0.3, this oscillatory behavior disappears (hence in any high resolution run, where the algorithm is most likely to be residual velocity limited rather than shear-limited, it never turns up). Second, the inner grid has strongly radially elongated cells. If we take a log-grid (see e.g. Nelson 1999, or the example below), where the cells are almost "square'' everywhere, this behavior is not observed, even if the run remains shear-limited. And finally, we have a steep density gradient close to the inner boundary. If the inner boundary was closed and hence if we had no density gradient, this oscillatory behavior would never appear. In all the cases where it was observed this behavior always disappeared after a few tens of dynamical times.

It should be noted that the numerical damping observed in the non-stationery frame in Sect. 4 occurs both in the non-rotating frame and corotating frame (far from the coorbital region) standard method runs. Hence the amplitude of the protoplanet triggered density wave is marginally higher in a FARGO run at the inner boundary. Both this reason and the effect we noticed in the previous paragraph lead to a marginally higher mass loss through the inner boundary, at least during the first stages of the evolution of the system, which results in the darker band at the inner boundary in the middle panel of Fig. 3.

We present in Fig. 4 the results of three runs (non-rotating standard and FARGO, and corotating standard), which describe the same physical system as before after the same amount of time, but with a grid for which $N_{\rm r}=70$ , $N_{\rm s}=180$ , $R_{\rm min}=0.25$ and $R_{\rm max}=2.5$ , and with a geometric sequence for $(R_i)_{i\in[0,N_{\rm r}]}$ (hence it is a log-grid, and everywhere its cells are almost "square''). One can check on these plots that there is no oscillatory behavior in the FARGO results (this time the cells are no radially elongated near the inner boundary), whereas the run is still shear-limited. Furthermore, as stated above, a careful look at the inner spiral structure shows that it has a slightly higher amplitude in the FARGO case.

$\begin{figure}\par\psfig{file=threecuts.ps,width=\columnwidth} \par\end{figure}$

Figure 5: Disk density cuts at the planet radius. The solid line represents the FARGO transport case, the dashed line represents the standard case, and the dot-dashed line represents the corotating frame result. The dotted line indicates the unperturbed surface density. Note that the local maxima at $j\simeq 46$ and $j\simeq 94$ correspond to a temporary residual accumulation of material at the L₄ and L₅ Lagrange points of the protoplanet

From the results depicted in Figs. 3, 4 and 5, one can deduce that the FARGO transport algorithm on this particular problem is much closer than the usual standard transport algorithm to the exact solution (which must closely resemble the results given by the corotating frame run, at least in the coorbital region, since in Sect. 4 we have seen that one needs to be in the comoving frame in order to get accurate results even in the limit of a vanishing timestep). Another quantitative evaluation of the FARGO algorithm consists in monitoring the accretion rate onto the planet as a function of time. We present in Fig. 6 the accretion rate onto a one Jupiter mass protoplanet embedded in a minimum mass protostellar disk with no initial gap. The disk parameters are the same as before, as well as the grid resolution (arithmetic radial spacing with $N_{\rm r}=49$ and $N_{\rm s}=143$ ). Three runs are presented with three different schemes: the standard method in the rotating frame, which gives, according to Sect. 4, the most accurate results, the standard method in the non-rotating frame, and the FARGO method in the non-rotating frame. We use a slightly different accretion procedure than the one described by Kley (1999).

$\begin{figure}\par\psfig{file=monitoracc.ps,width=\columnwidth} \par\end{figure}$

Figure 6: Accretion rate as a function of time onto a one Jupiter mass protoplanet with three different methods. See text for details

5 Two-dimensional example: The embedded protoplanet problem