A&A 459, L17-L20 (2006)
DOI: 10.1051/0004-6361:20066304

LETTER TO THE EDITOR

Halting type I planet migration in non-isothermal disks

S.-J. Paardekooper1 - G. Mellema2,1


1 - Leiden Observatory, Postbus 9513, 2300 RA Leiden, The Netherlands
2 - Stockholm Observatory, AlbaNova University Center, Stockholm University, 106 91 Stockholm, Sweden

Received 28 August 2006 / Accepted 20 September 2006

Abstract
Aims. We investigate the effect of including a proper energy balance on the interaction of a low-mass planet with a protoplanetary disk.
Methods. We use a three-dimensional version of the RODEO method to perform hydrodynamical simulations including the energy equation. Radiation is included in the flux-limited diffusion approach.
Results. The sign of the torque is sensitive to the ability of the disk to radiate away the energy generated in the immediate surroundings of the planet. In the case of high opacity, corresponding to the dense inner regions of protoplanetary disks, migration is directed outward, instead of the usual inward migration that was found in locally isothermal disks. For low values of the opacity we recover inward migration and show that torques originating in the coorbital region are responsible for the change in migration direction.

Key words: radiative transfer - hydrodynamics - methods: numerical - planets and satellites: formation

1 Introduction

Migration has become a standard ingredient in planet formation theory, providing an explanation for the extrasolar giant planets that were found orbiting very close to their central star, the so-called Hot Jupiters. It is believed that these planets were formed at a distance of several Astronomical Units (AU), after which they moved in due to tidal interaction with the surrounding protoplanetary disk.

Depending on the mass of the planet, three different migration modes can be distinguished:

1.
Type I migration (Goldreich & Tremaine 1979; Ward 1997), valid for low-mass planets. The disk response is linear in this case, and the resulting torques can be calculated in a semi-analytic way (Tanaka et al. 2002). Planets that are less massive than a few Earth masses are subject to this mode of migration;
2.
Type II migration, for which the planet orbits in a deep annular gap (Lin & Papaloizou 1986). In this case, the planet is tidally locked inside the gap, and the time scale for inward migration equals the viscous time scale, which approximately equals 106 yr. This mode of migration holds for planets comparable to or more massive than Jupiter. See, however, Rafikov (2002) on gap formation for low-mass planets;
3.
In between these two regimes, interesting things happen. The onset of non-linear behavior can influence the torque to a large extent (Masset et al. 2006a). However, when the disk is massive enough, the planet enters a new migration mode that can be very fast and that relies on dynamical corotation torques (Artymowicz 2004; Masset & Papaloizou 2003). This is called type III migration, and it holds for planets that open up a partial gap.
In this Letter, we focus on deeply embedded low-mass planets and, therefore, on the regime of type I migration.

The time scale for type I migration is inversely proportional to the mass of the disk and the planet (Tanaka et al. 2002), and it can be much shorter than the disk lifetime of approximately 107 yr. A planet of 1 Earth mass (1  ${{M}_\oplus }$), for example, embedded in the minimum mass solar nebula (MMSN) at 5 AU, would migrate into the central star within 106 yr. Within the widely accepted core-accretion model of giant planet formation (Pollack et al. 1996), the gas giants also pass through a stage where the core is a few ${{M}_\oplus }$, which thus should disappear into the central star even faster. Therefore, some sort of stopping mechanism is required for planets to survive type I migration.

Proposed stopping mechanisms include magnetic turbulence (Nelson & Papaloizou 2004), a toroidal magnetic field (Terquem 2003), and an inner hole in the disk (Masset et al. 2006b; Kuchner & Lecar 2002). However, most analytical and numerical work on planet migration has focused on disks with a fixed radial temperature profile. For these disks, hydrodynamical simulations can reproduce type I migration in both two and three dimensions (Bate et al. 2003; D'Angelo et al. 2003,2002).

The isothermal assumption is only valid if the disk can radiate away all excess energy efficiently. When the disk is optically thick, the radiative energy flux through the surface of a sphere of radius Haround the planet is given by $F_{\rm R}=\sigma T^4/\tau$, where $\sigma$ is the Stefan-Boltzmann constant and $\tau$ is the optical depth over a distance H. Therefore, when the opacity increases, less energy can be radiated away. This makes the cooling time scale much longer than the dynamical time scale in the dense inner regions of protoplanetary disks, and therefore the isothermal assumption is likely to be invalid there.

One may wonder what effects releasing the isothermal assumption may have on the torque. The effects of a more realistic temperature profile on the Lindblad torques were investigated by Jang-Condell & Sasselov (2005) and Menou & Goodman (2004), who found that the migration rate is very sensitive to the temperature and opacity structure and that in some cases migration may be very slow compared to the isothermal type I case. In this Letter, we present the first results of including a detailed energy balance on the migration behavior of low-mass planets in a radiation-hydrodynamical context. In Sect. 2 we discuss the disk model, in Sect. 3 we present the results, and we give a discussion of the results in Sect. 4. We then conclude in Sect. 5.

   
2 Disk model

We work in spherical coordinates $(r,\theta,\phi)$ with the central star in the origin, and the unit of distance is the planet's orbital radius. In the plots we also use a Cartesian coordinate frame with the central star in the origin and the planet located at (x,y,z)=(-1,0,0). The planet stays on a fixed circular orbit throughout the simulation. The disk is three-dimensional, and the computational domain is bounded by r=0.4, r=2.5, $\theta=\pi/2$, and $\theta=\pi/2-5h/2$, where h is the relative scale height of the disk (h=H/r, with H the pressure scale height). The disk spans the full $2\pi$ in azimuth. For the radial and meridional boundaries, we use non-reflective boundary conditions (Paardekooper & Mellema 2006), except for the boundary at $\theta=\pi/2$, which is fully reflective. The azimuthal boundary conditions are periodic. This domain is covered by a grid with 256 radial cells, 768 azimuthal cells, and 16 meridional cells. This configuration leads to cubic cells near the planet. On top of this main grid, we put 5 levels of adaptive mesh refinement (AMR) near the planet. Each level increases the local resolution by a factor of 2, which leads to a local resolution of $\Delta r \approx 0.00025$. We focus on a planet of 5  ${{M}_\oplus }$ located at 5 AU from the central star, for which the Hill radius is then covered by 70 cells. While this resolution is not needed for accurate torque calculations, it is needed for correctly modeling the planet's envelope and for the process of accretion, which we will consider in a future publication. In the torque calculations, we neglect the material orbiting inside the Hill sphere. For deeply embedded objects, material could be bound in the smaller Bondi sphere (Masset et al. 2006b), so that we may exclude a region that is too large. However, we checked that the results are not sensitive to this.

The disk is taken to be inviscid, so as to filter out effects of viscous heating. The density follows a power law initially with index -3/2. We take h=0.05, a typical value used in simulations of disk-planet interaction. This sets the temperature at approximately 50 K at 5 AU. The disk is in Keplerian rotation initially with a slight correction for the radial pressure gradient. The potential contains terms due to the star, the planet, and the acceleration of the coordinate frame due to the planet and the disk. The last two contributions are of negligible importance for a planet of 5  ${{M}_\oplus }$. We smooth the potential of the planet over 2 grid cells.

We solve the flow equations using a three-dimensional version of the RODEO method (Paardekooper & Mellema 2006), with an energy equation included. Radiative transfer is treated in the flux-limited diffusion approach (Levermore & Pomraning 1981), together with the flux limiter of Kley (1989). The three-dimensional hydrodynamics solver (including AMR) was tested through isothermal disk-planet interaction, which showed that we could reproduce type I and type II migration in the relevant mass regime. The radiative transfer module was tested against multi-dimensional diffusion problems. We used opacity data from Bell & Lin (1994). In the temperature range of interest, the opacity is proportional to the density and temperature squared. We vary the initial midplane density to study different cooling regimes while keeping the initial temperature fixed.

The code is parallelized using MPI, and the simulations were run on 126 1.3 GHz processors at the Dutch National Supercomputer. Each radiation-hydrodynamical model took approximately 4 days of computing time.

   
3 Results

In Fig. 1 we show the evolution of the total torque on our 5  ${{M}_\oplus }$ planet for three different midplane densities, where $\rho_0=10^{-11}$  ${\rm g~cm^{-3}}$ corresponds to the MMSN at 5 AU. For this high density, the torque is positive, indicating outward migration. The absolute value of the torque is approximately a factor 2 lower than the type I torque. For lower densities, and therefore for lower opacities, we recover inward migration; and for the lowest value of the density, we approach inward type I migration of isothermal disks.


  \begin{figure} \par\includegraphics[width=8.4cm,clip]{6304fig1.eps}\end{figure} Figure 1: Total torque on a 5  ${{M}_\oplus }$ planet as a function of time for three different midplane densities, together with the isothermal result. The torques are normalized to the analytical value found by Tanaka et al. (2002), which is reproduced by the isothermal simulation. For high densities (and thereby for high opacities) the torque becomes positive,indicating outward migration.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=8.4cm,clip]{6304fig2.eps} \end{figure} Figure 2: Radial torque distribution for a 5 ${{M}_\oplus }$ planet for two different densities, together with the result for a locally isothermal equation of state. Both the inner wave and the outer wave are reduced in strength for higher densities, but the major difference comes from the corotation region.
Open with DEXTER

To find out where this positive torque originates, we show the radial torque profile close to the planet in Fig. 2. An embedded planet excites two spiral waves into the disk. The inner spiral wave is responsible for the positive bump near r=0.95, while the outer spiral wave creates the negative bump near r=1.05. The location and the strength of these bumps varies between the different runs, but the total Lindblad torque remains negative for all simulations. For  $\rho_0=10^{-11}$  ${\rm g~cm^{-3}}$, the magnitude of the Lindblad torque is a factor of 2 smaller compared to the isothermal case. This implies that it is the corotation region, located between r=0.95 and r=1.05, that is responsible for the change in sign of the total torque.


  \begin{figure} \par\includegraphics[width=8.4cm,clip]{6304fig3.eps} \end{figure} Figure 3: Vertical torque distribution for an adiabatic simulation with no radiative cooling compared to the locally isothermal result. The difference starts high above the Hill sphere of the planet, which is approximately at z=0.017.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=7.8cm,clip]{6304f4c.eps} \end{figure} Figure 4: Temperature at z=0.02, slightly above the Hill sphere of the planet. Overplotted is the direction of the velocity field, clearly showing the horseshoe orbits. The point where the two horseshoe legs meet is a region of compression, meaning the temperature is slightly higher at that location.
Open with DEXTER

We can pin down the origin of the positive corotation torque further by looking at the vertical torque distribution. In Fig. 3 we compare an isothermal run to an adiabatic simulation. The latter does not include radiative cooling, so the effect of inefficient cooling is even more pronounced. The important thing to note is that the difference in Fig. 3 starts well above the Hill sphere of the planet, the edge of which is located at approximately z=0.017; therefore, we can filter out the effects of heating deep within the planetary envelope and look at temperature differences above z=0.017.


  \begin{figure} \par\includegraphics[width=7.8cm,clip]{6304f5c.eps} \end{figure} Figure 5: Temperature at x=-1 (or r=1), showing the hot plume behind the planet (y>0, the planet moves towards the left of the figure). This plume extends all the way to z=0, but there the heating inside the planetary envelope dominates the temperature structure. For z>0 we show a simulation without a global pressure gradient (i. e. the gas orbits at the Kepler velocity), while for z<0 we show a simulation with a global radial pressure gradient ( $p\propto \rho T \propto r^{-2}$). The difference in the position of the warm plume is apparent for both simulations. The vertical line at y=0 is there to guide the eye. The color scale has been chosen to focus on the warm plumes.
Open with DEXTER

We show a cut through z=0.02 in Fig. 4, showing the temperature and the direction of the velocities. The horseshoe orbits clearly stand out, and where the two horseshoe legs meet, the gas is compressed and the temperature rises. Due to the slightly sub-Keplerian velocity of the gas, the two horseshoe legs meet behind the planet (see Fig. 5). As the disk tries to maintain pressure equilibrium, the density will be lower in this region. We verified that the temperature rise is enough to account for the positive torque by verifying that

\begin{displaymath}\int_{\rm disk} \frac{{\rm G}M_{\rm p}\partial \rho}{\vert{\v... ...rm {p}} \times ({\vec r}-{\vec r}_{\rm p}) {\rm d}^3{\vec r}, \end{displaymath} (1)

with $\partial \rho=-\rho_0 \partial T/T_0$, is within the same order of magnitude as the total torque on the planet. Equation (1) reproduced the total torque to within 25%, which shows that indeed the temperature change is strong enough to induce a density asymmetry that is capable of producing the observed positive torque. It is this density asymmetry that leads to a large, positive contribution to the torque. Only if the disk can radiate away this heat efficiently does the total torque become negative again.

   
4 Discussion

The corotation region is very sensitive to local conditions (Masset 2002,2001), including the shape of the planetary envelope. This may affect the contribution to the total torque. However, we have found that the positive corotation torque is a robust effect, showing no strong dependence on numerical resolution, for example. Because we take a smoothing length of 2 grid cells, we have also found no dependence on the exact form of the planetary potential.

We estimate that the transition between inward and outward migration happens around 15 AU in the MMSN, based on the simulations with different opacities. Although we parametrized the opacity using the local gas density, it is the opacity that matters. If the opacity is lowered, as when it is due to grain growth (Dullemond & Dominik 2005), the transition radius will be located farther inward.

The computational costs of the radiation-hydrodynamical simulations is too high to run simulations for hundreds of orbits. Therefore, effects that happen on the libration time scale, such as saturation of the corotation torque (Ogilvie & Lubow 2003), are not captured by these simulations. Adiabatic runs showed that indeed saturation occurs on the libration time scale, which considerably reduces the magnitude of the positive torque. This is another confirmation that is indeed the corotation region that is responsible for the positive torque. This does not mean, however, that the positive torque should disappear on a libration time scale. When the radiation diffusion time scale is shorter than the libration time scale, the asymmetry in temperature, and therefore in density, can be maintained. In the MMSN, this is the case for approximately r>1 AU.

For lower-mass planets, the horseshoe region shrinks, so the positive torque should decrease in strength. However, in this case the Lindblad torques also decrease in magnitude. We have verified that even for a planet of 0.5  ${{M}_\oplus }$ the total torque is still positive for $\rho_0=10^{-11}$  ${\rm g~cm^{-3}}$. Further study is required to pin down the exact dependence of the torque on planetary mass. We expect that in the high-mass regime, when the planet starts to open up a gap and when the density around the planet decreases by several orders of magnitude, cooling will always be efficient. Therefore, type II migration should not be very different from the isothermal case (see also Klahr & Kley 2006).

   
5 Conclusions

In this Letter, we have presented the first radiation-hydrodynamical simulations of low-mass planets embedded in a protoplanetary disk. We find that the migration behavior of a 5  ${{M}_\oplus }$ planet depends critically on the ability of the disk toradiate away the heat generated by compression close to the planet. When radiative cooling is not efficient, the torque on the planet is positive, indicating outward migration. Analysis of the torque distribution showed that it is the coorbital region that generates the large positive contribution to the torque, while the Lindblad torque is reduced in strength, but is still negative. Local, adiabatic simulations showed that compression at the turn-over point of the horseshoe orbits creates a warm region behind the planet. As the disk tries to maintain pressure equilibrium, the density behind the planet will become lower, which leads to a positive corotation torque. This happens in the dense inner parts of protoplanetary disks, where the opacity is high. In regions of low opacity, we recover inward type I migration. Therefore, low-mass planets do not migrate inward all the way to the central star, bun instead stop migrating when radiative cooling becomes inefficient. In the MMSN, the transition occurs at approximately 15 AU. As the disk slowly accretes onto the central star, this radius will move inward; as a result, low-mass planets will also eventually continue their inward migration, but on the much longer viscous time scale. This way, low-mass planets are saved from fast inward type I migration into the central star.

Acknowledgements
We thank the anonymous referee for an insightful report. S.P. thanks Yuri Levin, Mordecai-Marc MacLow, and Peter Woitke for useful discussions. We thank Willem Vermin for his assistance at the Dutch National Supercomputer. This work was sponsored by the National Computing Foundation (NCF) for the use of supercomputer facilities, with financial support from the Netherlands Organization for Scientific Research (NWO).

References

 
Copyright ESO 2006