5 Two-dimensional example: The embedded protoplanet problem

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*_{0}=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
is equally spaced, with
;
The grid
has
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
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}=0.5.
We plot in Fig. 2 the quantity
after 2.86 orbits. This quantity
represents the effective CFL ratio.
With the standard transport algorithm this ratio
is bounded by *C*_{0}.

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
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*_{i}in the inner part, have shown to allow a speed-up by a factor 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.

Figure 4:
Disk density
;
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 , , and , and with a geometric sequence for (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.

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
and
).
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).

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

We see from the curve obtained in the corotating frame that the accretion rate is about after 400 orbits. This is in relatively good agreement with Kley's results, who gets slightly more than after 400 orbits in a similar run, but with a different grid resolution and a slightly different accretion protocol. We see from Fig. 6 that the accretion rate in the non-rotating frame, with a standard method, is smaller than in the rotating frame run, by a factor . The fact that the accretion rate is slower in this case was to be expected from the curves of Fig. 5. Now the run with the FARGO algorithm leads to an accretion rate which is between the rotating frame results and the non-rotating frame standard transport results, and which are closer to the rotating frame results. From these considerations again we see that the FARGO transport leads to a smaller error wrt the rotating frame results. The point here is that the FARGO transport algorithm is about one order of magnitude or more faster than the corotating frame standard transport run, and that the corotating frame is suitable only to the study of a protoplanet on a fixed circular orbit. From these remarks it clearly appears that the FARGO transport algorithm is particularly well suited to the study of the protoplanet orbit long-term evolution. FARGO has already been used to study the migration and mass accretion of a Jupiter sized protoplanet in a protoplanetary disk. It has been extensively tested against existing independent codes, which use the standard transport algorithm. It has proven to give very similar results, and the slight differences which remain between these codes and FARGO can all be understood in terms of FARGO's lower numerical diffusivity (Nelson et al.1999).

Copyright The European Southern Observatory (ESO)