1 - Observatoire de la Côte d'Azur, BP 4229, 06304 Nice Cedex 4, France
2 - UMR AIM, DSM/DAPNIA/SAp, Orme des Merisiers, CE-Saclay, 91191 Gif-sur-Yvette Cedex, France IA-UNAM, Apartado Postal 70-264, Ciudad Universitaria, Mexico City 04510, Mexico

Abstract
Context. Numerical simulations of planet-disk interactions are usually performed with hydro-codes that - because they consider only an annulus of the disk, over a 2D grid - cannot take into account the global evolution of the disk. However, the global evolution governs type II planetary migration, so that the accuracy of the planetary evolution can be questioned.
Aims. To develop an algorithm that models the local planet-disk interactions together with the global viscous evolution of the disk.
Methods. We surround the usual 2D grid with a 1D grid ranging over the real extension of the disk. The 1D and 2D grids are coupled at their common boundaries via ghost rings, paying particular attention to the fluxes at the interface, especially the flux of angular momentum carried by waves. The computation is done in the frame centered on the center of mass to ensure angular momentum conservation.
Results. The global evolution of the disk and the local planet-disk interactions are both well described and the feedback of one on the other can be studied with this algorithm, for a negligible additional computing cost with respect to the usual algorithms.

Key words: methods: numerical - solar system: formation - accretion, accretion disks

Planetary formation occurs in disks of gas and dust around protostars. Giant planets - whose mass is mostly made of hydrogen and helium - must have formed before the dissipation of the gas disk. Consequently, they must have exerted tidal forces on the gas and their orbits must have evolved in response to the gas.

The presence of a planet on a circular orbit leads to the formation of a spiral density wake in the disk. Goldreich & Tremaine (1980) and Lin & Papaloizou (1979) showed that, through the wake, the planet transfers angular momentum to the part of the disk exterior to its orbit, while it receives angular momentum from the inner part. For low mass planets, Ward (1997) showed that the net result is a loss of angular momentum for the planet, which makes its orbit decay on a short timescale. This is usually referred to as type I migration.

As the planet transfers angular momentum to the outer part of the disk, it repels it outward; it symmetrically pushes the inner part inward. If the planet is massive enough (the threshold mass depends on the disk's viscosity and scale height; see Crida et al. 2006, and references therein) a clear gap opens around the planet's orbit, effectively splitting the disk into an internal and an external part. In this situation, the planet is locked in the middle of the gap, because both the outer and the inner part of the disk are pushing it away. Thus, the planet moves in parallel to the migration of the gap (Lin & Papaloizou 1986; Ward 2003), which in turn has to follow the viscous evolution of the disk (characterized by a radial spreading of the disk as the gas is accreted onto the central star; see Lynden-Bell & Pringle 1974). This is usually referred to as type II migration.

Although analytic theories have brought a great deal of understanding of the fundamentals of planet-disk interactions, it has become increasingly evident that numerical simulations are an essential tool of investigation (see Papaloizou et al. 2005, for a review). Numerical simulations, however, are difficult and time-consuming, and they involve simplifications, which may produce artifacts. In this paper we are concerned with the simulation of type II migration.

Proper simulation of type II migration requires not only the correct calculation of planet-disk interactions - which are essentially local - but also of the global evolution of the disk. Unfortunately, planet migration is typically simulated with hydro-dynamical codes using a 2 dimensional polar grid which, for numerical reasons, is truncated at an inner and an outer radius. This enables to describe the local interaction of the planet with the disk but not the global evolution of the disk.

Indeed, the boundary conditions at the extreme rings of the polar grid cannot take into account what happens in the whole disk, outside of the grid. For instance, the use of open boundary conditions allows the gas to leave the region covered by the grid; this makes the studied disk annulus behave as if it were surrounded by vacuum, so that it empties very rapidly, which is not realistic. In the opposite extreme case, boundary conditions that impose that the mean density on the extreme rings remains constant with time disable the accretion and spreading of the gas. Some other prescriptions between these two extremes may be used (for instance allowing inflow, or imposing an outflow given by an analytical model, or setting the flow on the last ring equal to the viscous flow measured in the neighboring rings, etc.). However, it is very difficult to adapt these prescriptions to the changing behavior of the disk with time, in particular when the disk undergoes perturbation from the planets. Also these prescriptions are rather arbitrary and they may introduce artifacts in the planetary evolution.

In this paper we present a novel idea for the correct calculation of the global evolution of the disk. It consists of surrounding the 2D polar grid with a 1D radial grid. The 1D grid extends from the real inner edge of the gas disk (e.g. the X-wind truncation radius at a few tenths of AU) to the real outer edge (e.g. the photo-dissociation radius at hundreds of AU). This 1D grid has open boundaries at the inner and outer edges, and exchanges information with the 2D grid for the definition of realistic, time-dependent boundary conditions of the latter. Our algorithm for the interface between the 1D and 2D grids is driven by the requirement that the angular momentum of the global system (the disk in the 2D section, plus the disk in the 1D section plus the planet-star system) is conserved. We will describe it in detail in Sect. 3. First, however, we revisit the algorithm used to model the gas evolution and planet-disk interactions on the 2D grid, to ensure that it also conserves angular momentum, which is often not the case in standard implementations. We discuss this issue in Sect. 2. In Sect. 4, the results of our new algorithm are discussed, in terms of CPU time and robustness with respect to the positioning of the interfaces between the two grids. In Sect. 5, we describe some interesting astrophysical applications of this new hydrodynamical code.

2 Conservation of angular momentum in 2 dimensional hydrodynamical algorithms

Consider a simulation of the planet-disk interactions, with the disk represented on a 2D polar grid with open boundary conditions. A necessary requirement for the simulation to be correct is that the sum of the angular momenta of the disk, of the star-planet(s) system and of the gas outflowing through the boundaries remains constant over time.

Our simulation accounts for a Jupiter mass planet on an initially circular orbit. Here and in the rest of the paper we adopt the following units: solar mass, the initial semi-major axis of the planet, and its orbital frequency so that the gravitational constant G=1 and an orbit lasts $2\pi$ time units at r=1. The grid used to represent the disk extends from r=0.25 to r=3 with open boundaries. It is equally divided in $N_{\rm r}=165$ elementary rings, and $N_{\rm s}=320$ sectors. The disk aspect ratio (H/r) is set uniformly to 5%. It is assumed to be constant in time, hence the disk is locally isothermal. The equation of state used in FARGO is: $P={c_{\rm s}}^2\Sigma$ , with $c_{\rm s} = H\Omega$ as usual, where P is the pressure, $\Sigma$ the density, $c_{\rm s}$ the sound speed, H the disk scale height, and $\Omega$ the local angular velocity. The gas kinematic viscosity is $\nu=10^{-5.5}$ in our normalized units (which corresponds to the viscosity at the location of the planet for $\alpha=1.25\times 10^{-3}$ in a Shakura & Sunyaev (1973) prescription). Its mean density is $\Sigma=3\times10^{-4}$ , which is a bit less than the Hayashi Minimal Mass Solar Nebula at Jupiter (Hayashi 1981); the initial density profile can be seen in Fig. 5 and corresponds to a disk that evolved for some time under the effect of its own viscosity; it can be very well approximated by $\Sigma(r)=0.000306~\exp~(-r^2/52.8)$ .

The evolution of the total mass and angular momentum of the whole system (including the outflow, e.g. the cumulated mass and momentum advected outside of the region spanned by the 2D grid by the gas outflowing through the boundaries) is presented in Fig. 1. While the total mass is conserved at the level of numerical errors (top panel), the conservation of the angular momentum is quite poor (middle panel).

$\begin{figure} \par\includegraphics[angle=-90,width=8.5cm,clip]{5870_f01.ps}\end{figure}$	Figure 1: Relative variation of mass ( top panel) and angular momentum ( middle panel) of the system {gas + Jupiter + outflow} with time. The absolute variation of angular momentum of the system is compared to that of the planet on the bottom panel (in our units, in which the initial angular momentum of the planet is 10^-3).
Open with DEXTER

The bottom panel shows that the gain of angular momentum of the whole system (plain line), which is the error of the simulation, amounts to 10% of what the planet exchanges with the gas in its migration process (dashed line). Thus we can expect that the migration of the planet measured in the simulation is correct only at the 10% level.

This shows the need to improve the algorithm in order to achieve a much better conservation of angular momentum. We describe a procedure that considerably improves angular momentum conservation below.

2.1 Choice of reference frame

Need for a frame centered on the center of mass.

Implementation.

Caveats.

In the case of an axisymmetric problem (accretion disk with no planet), the choice of such a frame may lead to numerical instabilities because small deviations from axial symmetry in the disk due to numerical errors cause a shift of the star relative to the center of mass, which in turn can enhance the disk's asymmetries. Thus, if the planetary system is made of the star alone, we recommend to keep it fixed at the origin of the frame.

In many applications, the equations of motion are implemented in a rotating frame in order to keep the planet at a constant position. This can still be done in a frame centered on the center of mass.

2.2 Sub-steps sequence

As we anticipated above, the integration algorithm separately treats the pure gravitational part and the hydrodynamical part. To ensure a good preservation of the conserved quantities, it is necessary to respect as much as possible the action-reaction principle in the planet-disk interactions. This requires that the sequence of integration sub-steps is taken in a specific order.

2.3 Action-reaction symmetry

The action-reaction principle has to be perfectly fulfilled in the computation of sub-steps 2 and 3. It is not the case, for instance, if one computes the action of a planet on the disk from the gradient of its potential on each cell, and the action of the gas on the planet from the sum of the elementary gravitational forces that it feels from each cell. In principle, both calculations should be equivalent, giving total forces of the same strength and opposite directions. However, if the gradient is computed by finite differences (typical in grid calculations), this introduces a difference with respect to the force exerted by the corresponding cell on the planet.

In order to impose action-reaction symmetry, we proceed as follows. After sub-step 1, we measure the total change of angular momentum that the gas will have at sub-step 3 due to the potential of each planet. Then, in sub-step 2, we compute the change of azimuthal velocity of each planet according to the change of angular momentum that it causes to the disk. This ensures angular momentum conservation between sub-step 2 and sub-step 3. Similarly, in sub-step 2 the planet's radial velocity change is computed by measuring the total change of linear momentum in this direction caused to the gas disk by the same planet.

2.4 Results and discussion

All these precautions allow us to greatly improve the angular momentum conservation of the whole system. The same simulation as in Fig. 1 is computed with our modified algorithm. We denote by $\delta H$ the variation of angular momentum, in the inertial frame centered on the center of mass, of the whole system (the sum of the momenta of the gas in the 2D grid of the planet-star system and of the outflow which is the cumulated momentum carried by the gas that left the grid). It should be zero if the code were perfectly conservative. Thus $\delta H$ is a measure of the error of the integration scheme. It is compared with three relevant quantities in Fig. 2. The bottom curve (solid line) corresponds to the logarithm of the relative variation of angular momentum of the system: $\delta H/H_0$ , where H₀ is the initial total angular momentum (it corresponds to the middle panel of Fig. 1); on a long timescale (16 000 time units, $\gtrsim$ 2500 orbits at r=1), this normalized error does not exceed 10^-5.5. The middle curve (long dashed line) corresponds to the logarithm of $\delta H/H_{\rm p}$ , where $H_{\rm p}$ is the angular momentum of the planet: the error in total angular momentum is about 4.5orders of magnitude smaller than the angular momentum of the planet. Thus, it should not affect its migration. The top curve (short dashed line) shows the logarithm of the ratio between $\delta H$ and $\delta H_{\rm p}$ , where $\delta H_{\rm p}$ is the variation of angular momentum of the planet: this ratio is smaller than 10^-3.

$\begin{figure} \par\includegraphics[width=8cm,angle=270]{5870_f02.ps}\end{figure}$

Figure 2: Relative angular momentum variation as a function of time. The variation of angular momentum of the system $\{$ gas + planet + star + outflow $\}$ , $\delta H$ , is compared to its initial value H₀ (bottom curve, plain line), to the angular momentum of the planet $H_{\rm p}$ (middle curve, long-dashed line), and to the variation of the latter $\delta H_{\rm p}$ (top curve, short-dashed line).

Open with DEXTER

Although highly satisfactory for our purposes, the conservation is admittedly not perfect, in particular when compared to the mass conservation (which in the new simulation is as good as in the example of Fig. 1). This is due to two sources of error.

The most obvious one is that the algorithm used for the advancement of the planetary system, which is in our case a Runge-Kutta algorithm, is not symplectic. However, we checked that in our simulation, the amount of angular momentum artificially introduced in the system in this way is 3 orders of magnitude smaller than the total global error.

The second source of error is in sub-step 6, due to the translation of the center of mass of the planet-star system. This translation occurs for two reasons. (i) It is required to compensate the center of mass motion due to errors introduced by the discretization of the grid and to the second order (in time) advection of the disk. Thus, the amplitude of this effect decreases with increasing disk resolution and decreasing time-step. (ii) Gas continuously leaves the grid through its inner or outer boundary. This gas is deleted from the simulation, while its total mass and momentum are recorded as outflow. The problem is that the outflowing gas is in general not axisymmetric. Therefore, its elimination from the simulation causes a shift (in position and velocity) of the center of mass of the whole system. The shift is artificial, in the sense that it would not exist if the gas were modeled with an infinite 2D grid, and it requires, for compensation, the translation of the planetary system. This issue is a conceptual one. Thus this source of error cannot be reduced by tuning the grid resolution or time-step. This seems unavoidable: it is the price to pay to work with 2D grids that are not as extended as the physical disk. However, as we have seen in Fig. 2, this error is so small that for all practical purposes we can safely ignore it.

3 Coupling of the 2D grid with a 1D grid

As we discussed in the introduction, for a correct simulation of type II migration the conservation of angular momentum is a necessary, but not sufficient condition. It is also necessary that the global evolution of the disk is correctly reproduced, which is not the case if only a portion of the disk (the annulus represented by the 2D grid) is considered.

For this reason, we consider the whole disk, from its real inner to outer radius, corresponding respectively to the "X'' radius at the inner edge (Shu 2002), and to a photo-evaporation radius at the outer edge (see Hollenbach & Adams 2004b,a). Because the representation of the whole disk with a 2D grid would be computationally prohibitive, outside of the radial distance range covered by the usual 2D grid we study the evolution of the extended disk in one dimension, with the aid of a 1D grid representation. In practise, the 1D grid is made of two parts: one extending from the real inner edge of the disk to the inner edge of the 2D grid, and the other one extending from the outer edge of the 2D grid to the real outer edge of the disk.

In this section, we first present the equations for the 1D disk. Then we detail how we couple the 1D and 2D portions of the disk in a way that conserves mass and angular momentum, and that ensures a proper simulation of the disk's global evolution.

3.1 1D equations

The Navier-Stokes equation and the advection of the gas are computed in the 1D grid exactly as in the 2D grid, but with $N_{\rm s}=1$ , so that each cell can be considered as an elementary ring. Consequently, there are no azimuthal components for the various gradients, nor an advective transport phase in the azimuthal direction. The gravitational potential due to the planets and the star is computed for the 1D grid as a function of r only: the whole planetary system is considered as a single central mass, the mass of which is the sum of all the bodies inside the considered radius. Denoting by $\mathcal{H}$ the Heaviside function equal to 1 for positive arguments and to 0 for negative ones, the gravitational potential felt by a cell of the 1D grid located at a radius r from the center of mass reads:

$\begin{displaymath}\phi^{\rm 1D}(r)=\sum_{\rm p} -\frac{GM_{\rm p}}{r}\mathcal{H}(r-r_{\rm p}) \end{displaymath}$

We assume that the gas disk in the 1D grid is axisymmetric, so that in principle there is no torque in the interaction between the planets and the 1D grid. For more realism, however, one can assume that each 1D ring feels a gravitational torque due to each planet, using the formula of Goldreich & Tremaine (1980) or Lin & Papaloizou (1979):

$\begin{displaymath}\delta T_{\rm g}(r) \approx 0.4~\left(\frac{M_{\rm p}}{M_*}\r... ...{\rm p}}{\Delta}\right)^4 \left(2\pi r \Sigma~\delta r\right) \end{displaymath}$

(1)

Because of the assumption of axial symmetry, what leaves the 1D grid through its boundaries does not impose any re-adjustment of the star-planet(s) center of mass. Thus, the conceptual problem discussed at the end of Sect. 2.4 is not relevant here. Consequently, the evolution on the 1D grid always shows a perfect angular momentum conservation, at the level of numerical errors.

Despite this perfect conservation, however, for the results to be correct it is necessary that the assumption of axial symmetry is a good approximation for the real disk. Moreover, because Eq. (1) is an approximation, it is necessary that the planetary torque is small for the error to be small in absolute value. Both requirements are fulfilled if the interfaces between the 1D and 2D portions of the disk are placed sufficiently far from the planet(s). We will discuss this issue in Sect. 4, with quantitative tests.

3.2 Ghost rings

The 1D grid and the 2D grid communicate with each other through a system of ghost rings. This technique is derived from that used for the parallelization of a hydro-code on a distributed memory architecture. We describe it here for the inner interface between the 1D grid and the 2D grid, that is the interface at the inner edge of the 2D grid; the 1D grid is inside, extending from the real inner edge of the disk (the truncation radius or the surface of the star) to the 2D grid; the 2D grid lies outside this interface, as sketched in Fig. 3.

$\begin{figure} \par\includegraphics[width=8.8cm]{5870_f03.eps}\end{figure}$	Figure 3: Sketch of the coupling via ghost rings at the inner interface between the 1D grid and the 2D grid.
Open with DEXTER

Outside the interface, a number $N_{\rm G}$ of 1D rings are added, which are superimposed on the first $N_{\rm G}$ rings of the 2D grid; they are the ghosts of the inner 1D grid (see Fig. 3). At the beginning of each time-step, the density and the velocity in the 1D ghost rings are set as follows. The density of gas ( $\Sigma$ ), radial momentum ( $\Sigma v_r$ ), and angular momentum ( $\Sigma r v_\theta$ ) in the ghost ring are set to the azimuthal average of the respective quantities in the corresponding ring of the 2D grid. Denoting with the superscript "1D'' the quantities in 1D rings and by h the specific angular momentum ( $h=rv_\theta$ ), this leads to the following equations:

$\displaystyle \Sigma^{\rm 1D}$	=	$\displaystyle \frac{1}{2\pi r}\int_0^{2\pi}\Sigma ~r{\rm d}\theta\ ,$	(2)
$\displaystyle \Sigma^{\rm 1D} v_r^{\rm 1D}$	=	$\displaystyle \frac{1}{2\pi r} \int_0^{2\pi}\Sigma v_r~r{\rm d}\theta\ ,$	(3)
$\displaystyle \Sigma^{\rm 1D}h^{\rm 1D}$	=	$\displaystyle \frac{1}{2\pi r}\int_0^{2\pi} \Sigma h~r{\rm d}\theta\ .$	(4)

Then, during the time step, the computation of all stages is performed in the 1D grid as well as in the ghost rings. The 1D grid "feels'' the ghosts via the Navier-Stokes equations. As the 1D ghost rings have been filled with quantities inherited from the 2D grid, the 1D grid behaves as if it felt the 2D grid outside the interface.

Symmetrically, inside the interface, a $N_{\rm G}$ number of 2D rings are added: the 2D grid ghosts. They are treated as normal 2D rings, except that at the beginning of each time step, the azimuthal means of surface density, of radial momentum and angular momentum density are set equal to the respective values measured in the corresponding 1D rings, as sketched with arrows in Fig. 3. So, the 2D grid effectively "feels'' the 1D grid inside the interface. Filling the 2D ghosts is more elaborate than in the 1D ghost rings case. For a given 2D ghost ring, one has to proceed in four steps: (i) store for each cell the surface density, the radial momentum, and the angular momentum; (ii) compute the azimuthal means of these quantities; (iii) subtract in each cell, for each of the three considered quantities, the difference between its azimuthal mean and its value in the corresponding 1D ring; (iv) find the velocity in each cell by dividing the new momenta by the new surface density.

The minimal number of ghost rings needed, $N_{\rm G}$ , depends on the numerical scheme: it is the number of rings that causally affects a given ring during a time-step, or equivalently the number of rings in which the information contained in a ring propagates during a time-step; it is often called the kernel. For instance in FARGO, $N_{\rm G}=5$ .

This way of coupling the 1D and 2D grids ensures a smooth transition for each of the quantities computed in the code. In particular, if there is no planetary perturbation, the global disk (2D + 1D parts) behaves exactly as predicted by the Lynden-Bell & Pringle (1974) equations. This can be seen in Fig. 4, which corresponds to Fig. 4 in the aforementioned paper: it shows the evolution of the density distribution at four times, for $\Sigma_{t=0}(r)=\delta(r\!-\!1)$ , with $\delta$ the Dirac distribution. As in Lynden-Bell & Pringle (1974), T^*denotes the viscous time: $T^*=6\nu t$ in our units. The little discrepancy between the theoretical profile (thin) and the numerical one (bold) at T^*=0.004 most likely comes from the fact that the initial distribution is not exactly a $\delta$ function; this little difference vanishes with time.

$\begin{figure} \par\includegraphics[width=8cm,angle=270]{5870_f04.ps} \end{figure}$	Figure 4: Viscous spreading of a ring. The theoretical curves (obtained by numerical solution of Eq. (14) of Lynden-Bell & Pringle (1974) with the boundary condition $\Sigma (0.02)=0$ ) are thin, while the numerical ones obtained with our code are bold. The vertical dashed lines show the position of the interfaces, where the solution provided by our code remains perfectly smooth.
Open with DEXTER

3.3 Fluxes at the interface

The fluxes of mass and angular momentum through the interfaces have to be perfectly conservative. What leaves the 2D grid has to enter the surrounding 1D grid, and vice versa.

3.3.1 Mass flux

The definition of $v_r^{\rm 1D}$ in Eq. (3) makes the mass flux through a centered circle of radius r computed in 1D the same as the one computed in 2D. Their theoretical expressions are:

$\begin{eqnarray*}F_m^{\rm 2D} & = & \int_0^{2\pi}\Sigma v_r~r{\rm d}\theta\\ F_m^{\rm 1D} & = & 2\pi r\Sigma^{\rm 1D} v_r^{\rm 1D}\ . \end{eqnarray*}$

$\begin{eqnarray*}{F_m^{\rm 1D}}_{\rm computed} & = & 2\pi r \Sigma^{{\rm 1D}*} v... ...2D}}_{\rm computed} & = & \sum_j \Sigma^*[j]v_r[j]~r\delta\theta \end{eqnarray*}$

To ensure the mass conservation, we solve the following equation instead of (3) for the ring next to the interface:

$\begin{displaymath}{F_m^{\rm 1D}}_{\rm computed} = {F_m^{\rm 2D}}_{\rm computed}. \end{displaymath}$

(5)

3.3.2 Angular momentum flux

Denoting with a prime sign the non axisymmetric components of the quantities ( $h' = h-h^{\rm 1D}$ , $v_r'=v_r-v_r^{\rm 1D}$ ), using Eqs. (2)-(4) the angular momentum fluxes in 1D and in 2D are:

$\begin{eqnarray*}F_h^{\rm 1D} & = & 2\pi r\Sigma^{\rm 1D}h^{\rm 1D}v_r^{\rm 1D}\... ... =F_h^{\rm 1D} + \int_0^{2\pi} \Sigma h' v_r'~r{\rm d}\theta\ . \end{eqnarray*}$

We face a degree of freedom problem. With only 3 free variables ( $\Sigma^{\rm 1D}$ , $v_r^{\rm 1D}$ , and $v_\theta^{\rm 1D}$ ), 4 quantities have to be set in the 1D ghosts: the mass and the angular momentum that is in each ring, the mass and angular momentum fluxes. The first two are set with the prescriptions (2) and (4), which determine the values of $\Sigma^{\rm 1D}$ and $v_\theta^{\rm 1D}$ . The third variable, $v_r^{\rm 1D}$ is naturally set by Eq. (3)/Eq. (5). It appears here that it is impossible to obtain the correct flux of angular momentum in the 1D grid ghosts.

A possible solution is that, during the advection phase in the 1D part of the algorithm, one imposes that the angular momentum flowing through the interface corresponds to the flux computed in the 2D grid, ${F_h^{\rm 2D}}_{\rm computed}$ . However, as the flux of angular momentum carried by a non axisymmetric wave cannot be represented with 1D Navier-Stokes equations, F_h' would be deposited abruptly in the first ring of the 1D grid. This would lead to the formation of a spurious gap at the interface, and possibly lead to a numerical instability.

It has been shown that pressure-supported waves travel through the disk and deposit their angular momentum smoothly, as they get damped viscously (Takeuchi et al. 1996) or through shocks (Goodman & Rafikov 2001). Thus, a better solution is to model this wave propagation through the 1D grid. We first perform advection as usual in each grid and evaluate at the interface $F_h' = F_h^{\rm 2D}-F_h^{\rm 1D}$ (or, better, $\delta F_h \equiv {F_h^{\rm 2D}}_{\rm computed}-{F_h^{\rm 1D}}_{\rm computed} \approx F_h'$ , which might be slightly different to the former according to the numerical scheme). Then, over a time-step $\delta t$ , we spread in the 1D grid the amount of angular momentum $\delta F_h~\delta t$ . A prescription for the deposition of the flux can be found in Goodman & Rafikov (2001) as a function of the distance from the planet. This could be used when there is only one planet, but not if there are several planets, because it is difficult to know what fraction of $\delta F_h$ is due to each planet. Thus, we adopt an exponential function, because it is scale free so that, once $\delta F_h$ is known at the interface, its deposition does not depend on the position of the planet(s) and is a function of the distance from the interface only. In practice, we assume that, over a time-step $\delta t$ , the angular momentum deposited by the wave in a ring of width ${\rm d}r$ located at a distance d from the interface between the 1D and 2D grids is:

$\begin{displaymath}\delta h=\delta F_h~\delta t ~\lambda \exp\left(-\frac{d}{\lambda}\right)~{\rm d}r , \end{displaymath}$

(6)

We assume that $\lambda=0.5$ in natural units, which is the value obtained by fitting with an exponential the decay of the wave-carried flux in Appendix C of Crida et al. (2006). For simplicity, we assume that $\lambda$ does not depend on the disk's viscosity and scale-height (an assumption partially justified in Goodman & Rafikov 2001).

The deposit of the quantity $\delta h$ of angular momentum in a ring is simulated by applying a suitable torque, namely by adding to $v_\theta^{\rm 1D}$ the quantity $\delta v_\theta^{\rm 1D}= \delta h/r$ , where r is the radial distance of the ring from the star.

Notice that the integral of (6) from d=0 to infinity is equal to the total angular momentum $\delta F_h~\delta t$ carried by the wave at the grid interface. However, the 1D disk is not infinite in radial extent. Thus, a fraction of the angular momentum carried by the wave will not be deposited in the disk, but will outflow from the system. This outflowing momentum, as well as the angular momentum and mass advected through the inner and outer radius of the 1D grid, are recorded, in analogy with what was done in Sect. 2 for the 2D grid alone.

3.3.3 Mean viscous torque

The coupling of the 1D and the 2D grid described in the two subsections above ensures a smooth, conservative evolution of the disk. The gas is free to accrete and spread from the star to its physical outer edge, through the interfaces between the grids. However, at the interfaces, the shear gives an azimuthal viscous stress. This appears as a torque exerted on a grid by its ghosts (and reciprocally). The torque exerted by the part of the disk inside a given radius r₀ on the part of the disk outside r₀ reads, from the Navier-Stokes equations:

For angular momentum conservation, the torque felt by the 1D grid from its ghosts must equal the one exerted by the 2D grid on its ghosts. This requires $T_{r\theta}^{\rm 1D}(r_{\rm interface}) = \frac{1}{N_{\rm s}}\sum_j T_{r\theta}[j](r_{\rm interface})$ , where the index j goes through the $N_{\rm s}$ cells of the 2D ring next to the interface. However, given that the expression of $T_{r\theta}$ implies the product of velocity gradients by the density, its average is not equal to the product of the averages, and the required equality is not necessarily true. Thus, one has to replace the value of $T_{r\theta}^{\rm 1D}$ at the interface by the average of $T_{r\theta}$ on the interface in the 2D grid.

3.4 Results and discussion

We consider the final state of the simulation of Sect. 2: a Jupiter mass planet initially at $r_{\rm p}=1$ evolves in a gas disk represented by a 2D grid extending from r=0.25to 3 for 16 000 time units ( $\approx$ 2500 orbits). The final density profile of the gas disk is shown in Fig. 5 as a dashed line. One can see that the planet has opened a gap and migrated inward (the planet is located in the middle of the gap). Moreover, the surface density of the gas over the considered range is strongly reduced relative to its initial value (dot-dashed line in Fig. 5). We then compute another simulation with the same planetary system and the same 2D grid, but introducing a 1D grid extending from 0.117 to 20 length units (approximately from 0.5 to 100 AU, assuming the unit length equal to 5 AU). The final profile of the gas distribution obtained in this new simulation is shown as a solid line in Fig. 5. Its bold part corresponds to the 2D portion of the disk.

We see two important aspects in Fig. 5. First, the radial profile of our new solution is clearly smooth, which indicates that there are no artifacts in the passage of information from the 2D grid to the 1D grid and vice versa. The small kink visible at the inner boundary of the 1D grid ( $r\sim 0.15$ ) is due to the implementation of the open boundary condition (see Appendix A); this artifact also appears in Figs. 4, 7 and 8. The change of sign of the second order derivative of the radial profile near the outer interface ( $r\sim 3$ ) is a real feature at the considered time-step, as we will discuss later (see Fig. 9 and the related discussion). Second, despite that the disk in the new simulation has also significantly evolved relative to the initial profile, the new surface density of the disk is very different from that obtained in the simple 2D simulation (compare the solid and dashed curves). In particular, in the new simulation the outer part of the disk has not significantly evolved, which is due to the long viscous time there, as the latter scales with $r^2/\nu$ . Conversely, in the old simulation, because of the open boundary condition, we observed a significant disk erosion. This shows the importance of the boundary conditions for the global evolution of the system.

$\begin{figure} \par\includegraphics[width=8.5cm,angle=270]{5870_f05.ps} \end{figure}$	Figure 5: Gas surface density profile after 16 000 time units $\sim$ 2500 orbits of the Jupiter mass planet initially placed at $r_{\rm p}=1$ . The solid line corresponds to a 2D (bold) and a 1D (thin) grid coupled, while the dashed line stands for a 2D grid alone. The dot-dashed line shows the initial profile.
Open with DEXTER

The global conservation in the system (namely the evolution of the mass and angular momentum of the planetary system, plus those of the gas in the two grids, plus those advected through the inner and outer radius of the 1D grid and recorded as outflow) is presented in Fig. 6, and compared to the conservation obtained in the code using the 2D grid only. One can see that the use of the extended 1D disk does not change the excellent result obtained in Sect. 2. Thus, the coupling between the two grids that we described above is correct in terms of conservation properties. Actually, the error in angular momentum is slightly larger in the case with the 1D grid because the density at the inner interface of the 2D grid does not vanish, as can be seen in Fig. 5; so, the major source of error (the non axisymmetric outflow from the 2D grid, discussed at the end of Sect. 2.4) does not disappear, while it does for the simulation using the 2D grid alone.

$\begin{figure} \par\includegraphics[width=8.5cm,angle=270]{5870_f06.ps} \end{figure}$	Figure 6: Relative variation of mass and angular momentum in the case of a 2D grid with open boundaries and a 2D grid coupled with a 1D grid.
Open with DEXTER

4 Performance of the code

In this section, we check the accuracy of the results, as a function of the size of the 2D grid, and we discuss the CPU cost of this hybrid scheme. We show that the coupling of a 1D grid increases by far the realism of the results of the hydro-code, for a low additional computation cost.

4.1 Position of the interfaces with respect to the planet

As already mentioned it in Sect. 3.1, the interface between the 2D and the 1D grids has to be sufficiently far from each planet so that the flux of angular momentum carried by the wave at the grid interface and the azimuthal dependence of the disk surface density are both small enough. An obvious test of our code is to check the dependence of the results on the positions of the interfaces between the two grids.

We ran simulations like the one described in Sect. 3.4, with different inner radii for the 2D grid, but keeping the grid resolutions constant. Figure 7 shows the density profiles of the inner disk obtained at the end of the computations (16 000 time units $\sim$ 2500 orbits). It appears that the differences among the results are small, even in the cases where the planet is not very far from the interface ( $r_{\rm p}\approx 0.75$ ). Only the cases with a transition radius greater than 0.5 show not negligible differences with the three other runs, which are remarkably similar to each other, in particular in the density profile of the inner edge of the gap. The case with an interface at 0.5833 is clearly not accurate, but note that it is not unrealistic, while the 1D grid starts only at 4 Hill radii from the planet.

$\begin{figure} \par\includegraphics[width=8.5cm,angle=270]{5870_f07.ps} \end{figure}$	Figure 7: Inner disk profiles at the end of simulations with different locations of the inner interface between the 1D and the 2D grid. The bold part of the profiles corresponds to the 2D grid.
Open with DEXTER

$\begin{figure} \par\includegraphics[width=8.5cm,angle=270]{5870_f08.ps} \end{figure}$	Figure 8: Inner disk profiles at the end of simulations similar to those of Fig. 7, but without the implementation of the wave damping algorithm (6). A strong sensitivity on $R_{\rm interface}$ appears.
Open with DEXTER

All cases, however, give a quite consistent representation of the surface density profile of the inner disk, which is - on the contrary - very different from that obtained by using the 2D grid alone (compare with the dashed curve in Fig. 5).

This would not be the case if we had not implemented in the 1D disk the exponential damping of the angular momentum carried by the density wave launched by the planet. For instance, Fig. 8 shows 3 of the simulations, recomputed by switching off the calculation of (6). One sees a kind of discontinuity at the interface, where the angular momentum deposition abruptly stops in the disk. This change implies a modification of the local equilibrium and of the shape of the density profile. Consequently, the results are strongly dependent on $R_{\rm interface}$ . This also highlights the importance for the gap structure of the flux of angular momentum carried away by the pressure supported wake. This flux is equivalent to the pressure torque studied in Crida et al. (2006).

The outer interface has a much smaller influence on the disk's profile than the inner one, because it is further from the planet in the studied case. In Fig. 5, it appears that the outer interface corresponds to a change of sign in the second order derivative of the density profile. Figure 9 shows that this is a real feature, specific to the chosen output time. Indeed, this figure shows the density profile at different times in the same simulation; its evolution convincingly demonstrates that the second derivative of the density at the interface ( $R_{\rm interface}=2.9167$ , right vertical dashed line) varies with time and can be non-zero. We also repeated the same simulation, moving the outer interface to $R_{\rm interface}=2.3167$ (left vertical dashed line). The obtained surface density profiles overlap exactly with those of Fig. 9 at the corresponding time, so they are not plotted. This test also demonstrates that the position of the outer interface has a negligible influence on the global disk evolution.

$\begin{figure} \par\includegraphics[width=8.5cm,angle=270]{5870_f09.ps} \end{figure}$	Figure 9: Outer disk profile close to the outer interface ( $R_{\rm interface}=2.9167$ , marked by a vertical line) at 3 different times. The vertical line at r=2.3167 marks the outer interface used in a second simulation, whose results are indistinguishable from those plotted in this figure.
Open with DEXTER

The evolution of the disk has a strong influence on the type II migration of the planet. Figure 10 shows the different migration rates in the simulations of Fig. 7, and in the simulation obtained with the 2D grid alone. Once again, the two first cases ( $R_{\rm interface}=0.3333$ and 0.4333) are remarkably similar, while the other ones show a slightly slower inward motion of the planet. However, the simulation done with only a 2D grid gives a migration that is sensibly faster, as expected due to the disappearance of the inner disk.

$\begin{figure} \par\includegraphics[angle=-90,width=8.8cm,clip]{5870_f10.ps} \end{figure}$	Figure 10: Type II migration of the Jupiter mass planet in the disk. The semi-major axis is plotted as a function of time for all the cases already discussed: various positions of the inner interface between the two grids, and no 1D grid.
Open with DEXTER

The results presented in this section show that our method for coupling the 2D and 1D calculations not only ensures the conservation of the angular momentum, but also allows a robust (i.e. weakly dependent on the grid interface position) modeling of the disk's structure and evolution.

4.2 Computational cost

In this subsection, we show how the use of a 1D grid allows the study of the whole disk for a negligible extra cost while the use of an extended 2D grid would be prohibitive.

From a theoretical complexity point of view, the number of elementary operations for the computation of a time-step is proportional to the number of cells of the considered grid. Thus, the computation of the disk evolution in the 1D grid is negligible with respect to the computation in the 2D grid. In addition, one has to consider the three operations required to couple the grids, which also need to be done at every time-step. They are the filling of the ghost rings, the computation of $\delta F_h$ and its spreading, and the change of $T_{r\theta}$ in the 1D grid. They require computation on $N_{\rm G}$ 2D rings or less. As the size of the ghost area is usually negligible with respect to the size of the 2D grid, the computational cost of the coupling is also negligible with respect to the one of a time-step in the 2D grid.

In the previous paragraph, we studied the number of elementary operations to be computed during a time-step. In most codes, the length of the time-step is adapted to the conditions imposed by the required resolution and the existing perturbation. In general, it is determined by the Courant-Friedrichs-Lewy (CFL) condition (Stone & Norman 1993, Sect. 4.6). The limit most likely comes from the cells with smallest radii, where the angular velocity is the largest. To avoid numerical instabilities, the time-step length $\delta t$ is set so that information does not propagate more than one cell over one time-step; thus one has $\delta t~\Omega < 2\pi/N_{\rm s}$ , where $\Omega$ is the angular velocity, which is about Keplerian in a gas disk. Consequently, denoting by $R_{\rm min}$ the inner boundary radius of the 2D grid, one has $\delta t \propto R_{\rm min}^{\;3/2}$ . This shows that extending the 2D grid toward the star shortens the time-step and slows down the simulation significantly. The use of the FARGO algorithm (Masset 2000a,b) enables one remove the mean angular velocity and to consider only the perturbed motion and the shear, but this leads to a similar conclusion, although the scaling of $\delta t$ with $R_{\rm min}$ is generally different.

The CFL condition in a 1D grid is much less constraining. Consequently, the addition of a 1D grid inside the inner edge of the 2D grid has no influence on the time step length, still determined by what happens in the 2D grid. It is thus important to increase $R_{\rm interface}$ as much as possible, to speed up the computation. Some benchmarking confirmed this reasoning (see Appendix B for more detailed results). The accuracy loss of the computation with the increase of $R_{\rm interface}$ has been discussed in Sect. 4.1, so that an acceptable trade-off can be found.

5 Astrophysical applications

The viscous evolution of the protoplanetary disk is thought to govern type II migration. Our code, which has been designed in order to correctly reproduce this evolution, while resolving the planet-disk interactions, is therefore very useful for studying any problem related to type II migration and the feedback exerted by the presence of planets on the evolution of the disk.

5.1 Outward migrating planets

A demonstration of how type II migration depends on the evolution of the disk has been provided by Veras & Armitage (2004). A disk evolving under its own viscosity accretes onto the central star, while spreading outward under the constraint of angular momentum conservation. Thus, at any instant in time there is a boundary in the disk within which the radial motion of the gas is negative and beyond which it is positive - see the right panel of Fig. 11 and Lynden-Bell & Pringle (1974). If a giant planet is located beyond this boundary it should move outward, as its migration has to follow the local evolution of the disk (see also Lin & Papaloizou 1986). Veras & Armitage (2004) showed this using a 1D model, where the disk evolved under its own viscosity, the torque exerted in the planet-disk interaction had the form (1), and the effect of waves carrying angular momentum flux was neglected.

Using our code, we perform more precise simulations of this process. The interaction between the planet and the disk is simulated in 2D, and the effect of waves is taken into account. This implies that the gap opened by the planet in the disk is less wide and deep than in the Veras & Armitage (2004) representation (Crida et al. 2006). Moreover, the planet also feels a corotation torque, which is otherwise neglected in the analytic expression of the planetary torque given by Eq. (1).

We find that the situation is not as simple as illustrated by Veras & Armitage (2004). For instance, the left panel of Fig. 11 shows the result of a simulation in which the planet is initially put at $r_{\rm p}=1.4$ , which is deeply in the outward spreading zone as shown by the right panel. Nevertheless, the planet migrates inward, in what seems to be a runaway type III migration (Masset & Papaloizou 2003). This is because the spreading disk forces enough gas to pass through the coorbital region of the planet so that the planet feels a strong corotation torque (Masset 2001) and decouples from the disk evolution.

We do obtain outward planet migrations, but only in specific cases. For instance when we place the planet initially outside of the spreading disk (which is unrealistic), or if we hold it for a sufficiently large number of orbits, so that it can open a deep gap that effectively truncates the disk at its inner edge (which, in practice, again places the planet outside of the spreading disk).

A detailed description of this mechanism and a comprehensive exploration of the parameter space are beyond the scope of this paper. However, we think this example shows that the evolution of a Jovian mass planet is not an ideal type II migration. Because the gap is not extremely clean, the planet is not fully locked in the disk's evolution. The planet feels the global motion of the disk (accretion or spreading) but at the same time it also feels a non-negligible corotation torque. Only a code like the one that we have developed in this paper can simulate the two effects correctly and hence allow a quantitative study of the planet's evolution.

$\begin{figure} \par\includegraphics[width=8.5cm,angle=270]{5870_f11.ps}\end{figure}$	Figure 11: Left panel: migration of a Jupiter mass planet initially put on a circular orbit at $r_{\rm p}=1.4$ . Right panel: radial velocity in the gas disk when the planet is introduced.
Open with DEXTER

5.2 Cavity opening

If the migration of a giant planet is strongly influenced by the evolution of the disk, the presence of the planet also influences the evolution of the disk, with a non negligible feedback on its own migration.

Perhaps the best example is that of the formation of an inner cavity in the disk. Once a planet has opened a wide and deep gap, the accretion time of the inner disk onto the star is smaller than the time scale of planet's migration, because of the negative torque exerted by the planet (Varnière et al. 2006). This leads to the disappearance of the inner disk. This in turn enhances the imbalance between the inner and outer torques felt by the planets, and accelerates the planet's migration. The presence of cavities are deduced in some protoplanetary disks, because the spectral energy distribution (SED) presents a lack of emissivity at the wavelengths corresponding to hot gas in the vicinity of the star (see for instance Calvet et al. 2005). Rice et al. (2003) suggest the SED of GM Aur could be explained by a cavity maintained by a two Jupiter mass planet.

Our algorithm clearly enables a more quantitative analysis of the cavity opening process than ever done before, for a low computational cost. As we have shown in Fig. 5, when the 1D grid extends down to the real inner edge of the disk, the evolution of the surface density of the inner disk is found to be slower than in a classical code with a truncated 2D grid, which makes the appearance of a deep cavity more difficult. Therefore, we think that the estimates of the planet mass required to achieve cavities of given depth, and the estimates of the lifetimes of these cavities, need to be revised. A detailed study of the cavity opening mechanism as a function of planetary mass, disk viscosity and aspect ratio is in progress and will be addressed in a forthcoming article.

6 Conclusion and perspectives

It is well known that the migration of giant planets (i.e. planets massive enough to open significant gaps in the disk's density distribution) is governed - or at least strongly affected - by the global evolution of the disk under its own viscosity.

The usual simulation algorithms solve the hydrodynamical equations over a 2 dimensional polar grid that is truncated at an inner and outer radius to keep the computing time within reasonable limits. Thus, they cannot correctly reproduce the evolution of the disk, nor the planet's migration.

In this paper, we have shown how the use of a 1D grid surrounding a classical 2D grid allows us to simulate the global evolution of the disk, while resolving the local planet-disk interactions with the accuracy of the usual algorithms. Coupling the two grids via a system of ghost rings, with special attention paid to the conservation of angular momentum leads to a smooth evolution of the disk from its inner radius (the truncation radius or the surface of the star) to its outer edge, several hundreds of AU away. This increases by far the accuracy of the simulation with essentially no additional CPU cost.

Consequently, this algorithm is the tool of choice to properly simulate type II migration and related problems. It also enables the study of the effect of the presence of planets on the disk's global evolution, which in turn affects the migration of the planets in a feedback effect.

In our simulations, we observe that the disks slowly disappear through the open boundaries of the 1D grid. However, disks are believed to disappear rapidly after only a few million years, under the action of photo-evaporation. This phenomenon could easily be introduced in our code via a simplified prescription consisting of removing a fraction of the gas in each cell, depending on time and location. This will open the possibility of simulating the evolution of planets in a globally evolving gas disk, over the disk's lifetime.

Appendix A: Open boundary condition

The "open'' boundary condition that allows outflow from the grid but not inflow is implemented as follows at the inner edge of the 1D grid.

In the innermost ring (number 0), the radial velocity is always 0 and the density is that of the neighboring ring. In the next (ring number 1), the radial velocity is set equal to that of the following ring (number 2) if and only if it corresponds to outflow from the grid (negative radial velocity), and to 0 otherwise. This gives:

Appendix B: Benchmarking

We ran 3 simulations of a Jupiter mass planet initially placed at $r_{\rm p}=1$ , evolving in a gas disk for 1000 time units on an Intel Xenon 2.66 GHz processor. The results, summed up in Table B.1 confirm our reasoning of Sect. 4.2.

Also, the computing time is not directly proportional to the number of cells, and is strongly increased by the higher shear in the vicinity of the star, otherwise the third simulation would not have been more that twice as long as the first one.

Third, the second simulation was slightly shorter than the first one because with the 1D grid, the density at the boundaries of the 2D grid does not tend to zero, the profile is smoother and consequently the CFL condition is less constraining there.

Simulating planet migration in globally evolving disks

1 Introduction