# 3D MHD Simulations of Planet Migration in Turbulent Stratified Disks

## Abstract

We performed 3D MHD simulations of planet migration in stratified disks using the Godunov code PLUTO, where the disk is turbulent due to the magnetorotational instability. We study the migration for planets with different planet-star mass ratios . In agreement with previous studies, for the low-mass planet cases ( and ), migration is dominated by random fluctuations in the torque. For a Jupiter-mass planet for , we find a reduction of the magnetic stress inside the orbit of the planet and around the gap region. After an initial stage where the torque on the planet is positive, it reverses and we recover migration rates similar to those found in disks where the turbulent viscosity is modelled by an viscosity. For the intermediate-mass planets ( and ) we find a new and so far unexpected behavior. In some cases they experience sustained and systematic outwards migration for the entire duration of the simulation. For this case, the horseshoe region is resolved and torques coming from the corotation region can remain unsaturated due to the stresses in the disk. These stresses are generated directly by the magnetic field. The magnitude of the horseshoe drag can overcome the negative Lindblad contribution when the local surface density profile is flat or increasing outwards, which we see in certain locations in our simulations due to the presence of a zonal flow. The intermediate-mass planet is migrating radially outwards in locations where there is a positive gradient of a pressure bump (zonal flow).

## 1 Introduction

Understanding why and how fast planets migrate is fundamental to explaining the observed distribution of exoplanets and constraining planet formation timescales and efficiencies (Alibert et al., 2005). The basic principle behind the migration of planets in protoplanetary disks is the transfer of angular momentum between the planet and its disk. This transport process occurs at Lindblad resonances and, in a locally isothermal disk, typically leads to fast inwards migration. (Goldreich and Tremaine, 1980; Papaloizou and Lin, 1984; Ward, 1986; Tanaka et al., 2002). This is the standard type I migration scenario which applies to low- to intermediate-mass planets where the specific torque is a linear function of the planet mass (Ward, 1997).

If the planet is massive enough (the mass depending on the viscosity and the pressure scale height), the tidal forces on the disk can eventually overcome the pressure gradient and the viscous transport, causing gap opening around the planet orbit. This migration regime, referred to as type II, in which the planet-disk interaction can no longer be described by linear perturbation theory, is then conceptually very different from the type I regime. Due to the gap opening, it is possible for the torques on the planet to cancel in such a way that the evolution of the planet’s position is determined by the viscous transport of gas in the disk, making the planet move with the disk on viscous timescales (Ward, 1997; Bate et al., 2003).

Current numerical and analytical calculations estimate migration timescales to be a small fraction of the expected disk lifetime, which creates a problem for the survival of planetary cores. Gas planet cores need to reach a critical mass before the onset of runaway gas accretion (Ida and Lin, 2008). It is well established that planet population synthesis models together with giant planet formation models require a much less efficient type I migration to reproduce the observed distribution of exoplanets (Benz et al., 2008; Mordasini et al., 2009; Alibert et al., 2005; Ida and Lin, 2004, 2007; Trilling et al., 2002). Nevertheless, core survival mechanisms have also been proposed to solve the timescale problem without resorting to an artificially reduced Type I migration rate (Terquem, 2003; Fromang et al., 2005; Thommes and Murray, 2006).

Deviations from linear theory have been found in a number of three dimensional calculations. Masset (2002) and D’Angelo et al. (2003) and later Masset et al. (2006a) found that for intermediate-mass planets (around ), the torques on the planet can be significantly lower and even reverse sign when the local surface density profile of the disk is flatter ( with ) than in the usually assumed Minimum Mass Solar Nebula (MMSN) model . This is found for a certain range of the disk viscosity. In this case, the torques from the corotation region can become important. The fluid elements that are librating (moving in horseshoe orbits in the corotation region) orbit on a U-turn trajectory at trailing and leading sides of the planet. These fluid elements exert a torque on the planet at each U-turn, which is a symmetric effect on both sides of the planet in an inviscid disk; therefore there is no net torque coming from this region after a few librating periods. This is referred to as the ”saturation” of the corotation torque. In the presence of viscosity, if the viscous crossing timescale across the horseshoe region of accreting disk material is smaller than the libration timescale, the torques exerted by fluid elements around the U-turn are not symmetric at each side of the planet, creating a net positive torque that can be sustained. This is refereed to as the ”unsaturated” corotation torque, and it can depend on the surface density and on the width of the horseshoe region that in turn depends on the planet mass (Ward, 1992).

So far most numerical studies of migration and/or gap formation have concentrated on the quasi-laminar disk case, where Navier-Stokes shear viscosity is included in order to model the viscous stresses resulting from including turbulence in the disk (e.g. Nelson et al. (2000); Bate et al. (2003); Bryden et al. (1999); de Val-Borro et al. (2006); Papaloizou and Larwood (2000); Crida et al. (2006) and many more). There has been strong interest in simulating planet-disk interactions in turbulent disks, where the turbulence is magnetically generated by the magneto-rotational instability (MRI) (Balbus and Hawley, 1991, 1998). Only ideal MHD has been considered so far in global simulations. The disk is assumed to be fully ionized and the magnetic diffusivity is negligible. Winters et al. (2003) looked at gap formation by intermediate- and large-mass planets in turbulent unstratified disks and the local internal stresses around the planet. In the MHD case, they found the gap to be shallower and wider compared to the laminar HD case; the Maxwell stresses in the disk dropped in the vicinity of the planet’s orbit. Papaloizou and Nelson (2003) performed a comprehensive study of protoplanets embedded in MHD-turbulent unstratified disks. They found that for low mass planets, Type I migration is no longer effective due to large fluctuations in the torque. No convergence was reached due to fluctuations of the torque on timescales longer than the orbital period and short simulation timescales. However, the torques for planets more massive than were found to converge to the standard Type I migration torques after long-time averaging (Nelson and Papaloizou, 2003; Papaloizou et al., 2004; Nelson and Papaloizou, 2004). For low-mass protoplanets, Nelson (2005) studied the long-term evolution of the orbital elements and particularly the excitation of eccentricity by turbulent fluctuations. The evolution of the orbital elements of particles in MHD turbulence has also been studied using shearing unstratified boxes (Yang et al., 2009) and stratified boxes including a dead zone (Oishi et al., 2007). To avoid the expensive MHD simulations, other approaches have been taken, such as modeling the turbulence as a time and space varying forcing in a laminar disk model (Laughlin et al., 2004). In this case, depending on the amplitude of the forcing, type I migration can be overcome by the random fluctuations in the torque, and random walk motion will be superimposed on the smooth inward migration. Baruteau and Lin (2010) used a similar turbulent forcing model and studied the unsaturation of the corotation torque due to turbulence. Depending on the amplitude of the turbulence, the corotation torque is found to be unsaturated to a certain level, making the total torque increase accordingly (become less negative), slowing down inwards migration. Other approaches include the analytical description of stochastic migration of low-mass planets using a diffusion-advection equation (Johnson et al., 2006; Adams and Bloch, 2009) and coupling N-body simulations with a random forcing to study the accretion and formation of low-mass planets (Ogihara et al., 2007). Recently, Nelson and Gressel (2010) examined the velocity dispersion of 1 m to 10 km planetesimals embedded in a turbulent disk, using 3D MHD simulations and neglecting stratification, and characterized the stochastic gravitational perturbations felt by planetesimals due to MHD turbulence.

In this paper, we study planet migration in stratified 3D MHD-turbulent disks for planet masses in the Type I and II mass range. In Section 2 we describe the numerical setup of our simulations, and the initial conditions for the disk and the magnetic field, before the addition of a planet. In section 3 we present the results of our simulations and finally in section 4 we discuss our results.

## 2 Simulation Setup

Simulations where performed using the finite volume fluid dynamics code PLUTO (Mignone et al., 2007). In the code, time stepping is done using a second order Runge Kutta scheme, while the spatial integration is performed using linear interpolation through the second order TVD scheme. The Riemann fluxes are computed using the HLLC and HLLD solvers for the HD and MHD cases, respectively. The code uses the Constrained Transport method for preserving a divergence-free magnetic field (Gardiner and Stone, 2005). The numerical setup for the MHD case follows the setup presented in (Flock et al., 2010). The MHD equations in the isothermal approximation (no energy equation) are given by

(1) | |||

(2) | |||

(3) |

The potential includes contributions from the star and the planet. We work in spherical coordinates , where the computational domain is given by , and . The grid resolution is and it is centered in the center of mass of the planet-star system. The boundary conditions for the velocities and magnetic field are periodic in the vertical ( boundary) and azimuthal directions and reflective in the radial direction, except for the transverse magnetic field component, which reverses its sign at the radial boundary. Buffer zones are defined at the radial boundaries to avoid boundary effects, where for the magnetic resistivity is given by and for the resistivity is .

### 2.1 Disk and Planet Models

As an initial condition we take a gas disk in sub-Keplerian rotation around a solar mass star. The azimuthal velocity is given by

(4) |

where is the Keplerian velocity and and are the exponents of the radial power law distribution of the density and sound speed . The initial density distribution is given by

(5) |

The disk is described by a locally isothermal equation of state . The ratio of the pressure scale height to the radial coordinate of the disk is taken to be a constant such that .

The gravitational potential of the planet is given by a softened point-mass potential

(6) |

where is the softening parameter, needed to avoid numerical divergence near the position of the planet and

(7) | |||||

is the distance between the planet and a gas particle in the disk. For all the simulations is set to be a fraction of the Hill radius with . Table 1 shows the parameters of our simulations. Distances are given in units of , density is given in units of , and velocity is given in units of Keplerian speed at , . The surface density have been scaled such that the total disk mass is . Magnetic fields are given in units of . In the cases where the planet is not on a fixed orbit (runs see Table 1), the equations of motion are integrated with a simple leap frog integrator. For the calculation of the torque, we include the entire disk in the integration (without Hill sphere tapering, except for simulation R9). The components of the torque vector in cartesian coordinates are given by

(8) |

where is any of the three cartesian indices and , with being the cartesian unit vectors. Of course, for studying migration, we are mostly interested in the component of the torque vector. Specific torques are given in units of .

### 2.2 Magnetic Field Configuration

Before introducing the planet in the simulations, a weak toroidal magnetic field is imposed on the disk given by

(9) |

where is the initial thermal pressure. This gives an initial azimuthal field with constant plasma beta . The field is imposed in a subset of the full computational domain given by and . The simulation is then followed until turbulence generated by the MRI has reached a saturated state. After this stage, we reset the density to the initial condition. This is the initial state in which the potential of the planet is incorporated and where all our runs start. The azimuthally, vertically and time averaged value of the effective parameter is shown in Figure 1. However, the stress is not constant throughout the vertical dimension. The upper layers of the disk are the most active. Figure 2 shows the time evolution of , and for simulation R0 (see Table 1) that does not include a planet. The top figure shows the characteristic butterfly diagram for the azimuthal component of the magnetic field in a turbulent stratified disk. ^{1}

#### Zonal flows and pressure bumps

The time-averaged thermal and magnetic pressure and the perturbed (with respect to Keplerian) azimuthal velocity are plotted in Figure 3 for run R0. We plot these quantities in the mid-plane of the disk, and one scale height above the mid-plane, for a simulation without a planet (or, equivalently, a massless planet). The radial gradient of the pressure has been removed and the pressure is averaged in the azimuthal direction. As expected of zonal flows, we see pressure bumps that correlate with bumps in perturbed azimuthal velocity, only phase shifted by one quarter of a period (Johansen et al., 2009). Bumps in thermal pressure correlate with drops in magnetic pressure, a behavior that is seen more clearly above the mid-plane, since the MRI is not resolved in the mid-plane of the disk. Notice also that in the velocity peaks, the azimuthal velocity exceeds the Keplerian value at some radial locations. These structure in the pressure and the velocity lived for the entire duration of our simulations, around 1000 inner orbits. These ”zonal flows” result from an inverse cascade of kinetic energy, e.g. a transport of energy from the MRI unstable medium scales, to the largest scales, which is very typical in accretion disks simulations (see for instance Dzyurkevich et al. (2010) and Lyra et al. (2008)).

## 3 Results

### 3.1 Disk torques and migration

Table 1 summarizes the computational time in local planet orbits for each of the simulations. The torque was calculated by taking into account the entire disk and its value was saved at every time step. We calculated the cumulative average specific torque as

(10) |

where is the cumulative average torque up to timestep and is the total time until timestep .

#### Low Mass Planets ( and )

Figure 4 summarizes the density structure of simulations R2, R5 and R9. Runs R1 (), R2 and R3 () shows no significant perturbation of the density by the planet, and no spiral arms are seen. The turbulent perturbations dominate in this case. Figure 5 shows the torque cumulative average torque as a function of local orbits for run R1. Figure 6 shows the torque for simulations R2 and R3. The fluctuations in the torque created by the perturbations in the density can, in both cases, be larger than the mean torque expected for standard Type I migration in a laminar disk. Comparing the torque for the planet at different positions in the disk, we see that the local (in time) evolution depends on the location of the planet. Random variations in the torque can be an order of magnitude larger than torques coming from the Lindblad resonances, in addition to the possibility that the spiral waves excited at Lindblad resonances are partly or totally suppressed by density fluctuations coming from the turbulence, such that the magnitude of this torque can be reduced. For the low-mass planet simulations, we find no convergence of the torque on timescale of the runs. For a run with a massless planet orbitting at , a gaussian fit of the time distribution of the torque gives a standard deviation of . We also calculate the auto-correlation function of the torque and take the correlation time to be given by . This gives local orbits, while the first and second zero crossing of the torque ACF occur at and local orbits. This is in agreement with results by Nelson (2005) and Fromang and Nelson (2009) and with estimates used by Baruteau and Lin (2010) that are based in previous MHD simulations of turbulent low-mass planet migration. We also calculate the power spectrum of the mid-plane density (see Eq. 3 in Baruteau and Lin (2010)). This is shown in Figure 7 and we compare our results from MHD simulations to Figure 1 of Baruteau and Lin (2010), where the spectrum in the result of the forcing model for the turbulence with . This comparison is valuable, since ultimately, a more complete parameter study migration and turbulence will have to be studied in models with forced turbulence. For this value, in our simulations we find that the larger scales carry more power than in the random forcing model, while the two spectrum agree for smaller scales for the case that includes the modes with . The higher power at larger scales can result from the higher compressibility of stratified disks, especially at large scales, as a vertically stratified disk can respond to compression with vertical expansion. However, the overall shape of the spectrum of the MHD simulation agrees better with the HD simulation without the modes included. Therefore with the proper scaling of the amplitude of the turbulence, and a cutoff of these modes, these simulations could reproduce the MHD spectrum. Another possibility is that the power at the small scales in the MHD simulation is lower due to lacking resolution at these scales. Ultimately, there needs to be a physical motivation for the cutoff of the turbulent forcing potential after the first few modes, if this is indeed the model that better matches the global MHD simulations.

#### Intermediate Mass Planets (, and )

Figure 4 shows the log density for run R5. For this simulation (), spiral arms are visible and their amplitude is comparable (or larger closer to the planet) to that of the perturbations generated by the turbulence. Figure 8 shows the cumulative average specific torque for run R4 and Figure 9 shows the torque for simulations R5 and R6. We see that in these three simulations there is an initial stage where the torque is negative followed by a reversal of the direction of the migration where the torque becomes positive and takes a defined value for the rest of the simulation. This happens at different times when we compare two different positions of the planet in the disk (runs R5 and R6). Instead of a random walk variation in semi major axis superimposed on smooth inwards migration, we find that planets of around 30 Earth masses undergo systematic outward migration. This outward migration is sustained for the total duration of the simulation. The simulation times for runs R5 and R6 are around 600 to 1000 orbits at the inner boundary of the disk (1AU). During this time, the density profile in the disk can evolve significantly from the initial state, and although the surface density profile can still be fitted by the initial profile (), there can be changes in the local profile at the position of the planet and accumulation of mass at the disk inside the planet’s orbit due to turbulent stresses. However, for both simulations, the torque is reversed before there is a significant accumulation of mass at the inner boundary and it converges to a constant value for the remaining simulation time. The torque for run R8 is shown in Figure 11. In this run the planet mass () is now able to modify the density profile around its orbit, and opens a partial gap, which affects the convergence of the torque. We don’t find convergence for the simulation time, but there is still a tendency for outwards migration.

Unlike the simulations for the small-mass planets (R1, R2 and R3), for the simulations R4, R5 and R6, the hill radius of the planet and the horseshoe region are resolved (by approximately 4, 4 and 7 grid cells per half width respectively). In this case, the component of the torque originating from the horseshoe region can dominate if there is a mechanism for keeping the corotation torque unsaturated and the local density profile differs from the global profile, possibly increasing outwards, such that the corotation torque can be larger than the Lindblad torque, making the total torque positive. There are special locations in the disk where is possible for the surface density to increase outwards, due to the appearance of zonal flows, as seen in Section 2.2.1.

Comparing the torque values of the simulations with analytical estimates by Paardekooper and Papaloizou (2009) or Masset et al. (2006a) is not straightforward. First, the undergoing evolution of the disk can make the surface density profile at the position of the planet and the effective stress resulting from turbulence vary in time, therefore one torque estimate does not apply at all times. On the other hand, the value of the horseshoe drag is very sensitive to the structure of the horseshoe region, and the estimate used in the analytical calculations is based on a 2D model of the flow around the planet. In our case, the horseshoe region is distorted, making the half-width difficult to define. Additionally, the inclusion of magnetic fields can introduce new magnetic resonances that affect the total torque, as seen by Terquem (2003) and Fromang et al. (2005) for a uniform non-turbulent field. For the sake of the comparison and simplicity, we discard this type of contribution. We attempt a comparison with the analytical estimates for the torque, including a contribution of the horseshoe drag. We take the total torque to be composed of the Lindblad torque in a 3D locally isothermal disk (Tanaka et al., 2002)

(11) |

plus the fully unsaturated non-linear horseshoe drag^{2}

(12) |

Here is the angular frequency of the planet and is the surface density at the position of the planet. The cumulative average torque at the end of the simulation for run R5 is and we take the half width of the horseshoe region to be , as is measured in our simulations (calculated from the analytical expression, ). Assuming the global surface density profile , with , will always give a negative torque. However, the torque always becomes positive for and matches the simulation value for , which is comparable to the local profile observed in the simulations (see middle plot in Figure 12). We should also note that already a ”close-to-flat” profile can significantly reduce the negative torque or change the sign of the torque. For simulation R6, the cumulative average torque at the end of the simulation is and the measured half width of the horseshoe region is . In this case, a local profile of is enough to obtain a positive torque, while a local profile of matches the value of the torque obtained in the simulation. However for this simulation, we only observe a flatter profile (see Figure 12). This discrepancy can be due to the fact that these values are very sensitive to the value of , since the horseshoe drag scales as , and we stress that the streamlines can be very distorted, therefore making the estimation of difficult. This is a critical parameter, and we find that an increase of to in the simulation value of with respect to the analytical estimation is enough to reproduce the observed positive torques. Therefore, if one assumes that the observed torque is composed of the wave torque plus the corotation torque and neglects any additional effect, we see that the corotation torque is crucial and able to cancel out or overcome the negative Lindblad contribution for standard disk parameters.

To further test the effect of the local density profile, we performed a simulation with , for the planet located at (run R7, Figure 10), initially at the right side of a pressure bump (where pressure and density decrease with radius). In this case, the cumulative average torque does not clearly converge, and we do not see systematic outwards migration, as the cumulative average torque approaches zero. However there is still a significant reduction of the torque as compared to the Type I Lindblad torque, which cannot be explained only in terms of a locally decreasing radial density profile. This result suggest that even in the absence of a pressure bump, inwards migration can be significantly slowed down for this planet mass. We note also that we used the expression for the horseshoe drag valid for an isothermal disk, so that there is an additional contribution due to the locally isothermal profile that we did not take into account.

To see if the transport of mass in the disk is enough to sustain the unsaturated torque, we take the expression for the minimal to mantain the unsaturated corotation torque (Masset et al., 2006a)

(13) |

we obtain for , which is always smaller than what we observe in our simulations^{3}

#### Large Mass Planet ()

The density structure and spiral arms induced by the Jupiter mass planet in run R9 () dominate over the turbulent perturbations and influences the entire disk. Figure 14 shows the cumulative average torque. We have excluded the torques coming from the Hill sphere for the calculation of the torque, to exclude material bound to the planet which is not properly simulated at this resolution or without including other relevant effects. For these simulations, initially we see the same trend as for runs R5 and R6, where the torque becomes positive, however this is not sustained for this planet mass as torques coming from the corotation region are suppressed due to gap opening (see bottom plot in Figure 12 and the right panels in Figure 4). Additionally, the planet modifies the stresses in the disk, and therefore, the accretion behavior, as can be seen in the bottom plot in Figure 13, that shows the evolution of the stresses for run R9. In this case, the stress is progressively suppressed and the Reynolds stress dominates the total stress. The reduction of the Maxwell stress is seen mostly in the part of the disk inside the planet’s orbit and in the gap region (see Figure 15). Is possible that the magnatic stress is suppressed in the gap region due to the modified azimuthal velocity near the planet, which can suppress the MRI locally. Notice that for this run, the density has a more laminar appearance (see Figure 4), consistent with a reduction of the stress due to the presence of the planet. This could also be a numerical effect that appears at this resolution, so further studies at higher resolution are needed.

In Figure 16 we compare the gap opened by the Jupiter mass planet in a magnetized disk (run R9) with an equivalent HD 3D simulation with viscosity where (run R10) and stratification. The time of the snapshots is 100 local orbits. The gap for the hydro case is narrower and slightly deeper than the gap formed in the magnetized turbulent disk. However the gap is not completely cleaned after this time. We observed the same characteristics for lower-mass planets that open only a partial gap in the disk.

Winters et al. (2003) studied a similar case of gap opening, but in a unstratified MHD-turbulent disk. In agreement with our results, they found a wider gap when the disk is turbulent, and larger transport of mass from the outer to the inner disk (see our Figure 16). In contrast to our findings, they find a deeper gap in the hydro case. This can be due to a different treatment of the gap opening criteria, since the planet mass in their calculations does not satisfy the viscous criterion. In terms of the reduction of the stresses around the planet, we find agreement with their results.

Nelson and
Papaloizou (2003) studied gap opening by a giant planet in an MHD-turbulent unstratified disk and compared their results with 2D simulations with an viscosity. They found that a run with an equivalent stress to the turbulent run produced a shallower gap. This is different to our own results where we found the same gap depth. However, we should note they only studied how the turbulence affects an already formed gap and did not observe the depletion of the outer disk. The difference to our simulation results could also be due to our choice of in the hydro run. We choose an matching the global average of the turbulent runs, but this is a quantity that varies vertically. In a non-stratified disk simulation, this averaging is not necessary.

## 4 Discussion and Conclusions

For simulations R1, R2 and R3, where , during the simulated time, migration was dominated by random fluctuations in the torque, that can be orders of magnitude larger that what is expected for the value of the Lindblad or corotation torques for this planet mass. This is in agreement with simulations by Nelson (2005) of migration of low-mass protoplanets in cylindrical disk models, where stratification is neglected.
It is unclear if after long term averaging ( orbits), the fluctuations will average out to zero while some component of the systematic torque will remain. Such a calculation is currently too expensive. It will also be difficult to get a steady state without special prescriptions for correcting the density, due to the accretion evolution of the disk, in addition to the decrease in stress for long simulation times due to the limited resolution. Another interesting point for further studies is to investigate this type of migration with enough resolution to resolve the corotation region, to see the impact of the corotation torque in these cases. However, even if this torque is present and well resolved, its magnitude would still be small compared to the amplitude of the fluctuations, since ultimately the torque depends strongly on the width of the corotation region, which approaches zero as the planet mass approaches zero.

As the planet mass is increased by one order of magnitude to , the hill radius is now properly resolved and the systematic torque is now large enough to dominate over the random component of the torque. Outwards migration in a locally isothermal disk can occur due to the viscosity unsaturating the torque coming from the corotation region (where the viscous timescale across the horseshoe region is smaller than the libration timescale), as was found by Masset et al. (2006a) for planets in the intermediate-mass range. Specifically, a planet with mass ratio in a disk with with a flat surface density profile () and an viscosity was found to be the critical mass for which the offset from linear theory was the largest. Additionally, a planet mass of is within the range of masses for which reversal of migration occurs, if one extrapolates their results to a disk with . However, this is the first study that observes the effect of the unsaturated corotation torque due to the accretion of mass in the disk provided directly by turbulence that has been self-generated by the MRI. For runs R4, R5 and R6, the planet is in locations in the disk where the local surface density profile is either close to flat or increasing outwards, due to the pressure bumps seen in Section 2.2.1. This can make the contribution of the corotation torque dominate over the Lindblad torque, which is usually unexpected in an disk, since for realistic density profiles, the Lindblad torque will dominate. This is also consistent with the torque in run R6 being initially negative, since at the beginning of the simulation the local profile is decreasing outwards, but getting shallower as time increases, eventually reaching the point where the torque reverses. For run R5, the slope is almost immediately increasing outwards due to the evolution of the disk, which makes the torque positive from the beginning of the simulations. We can only roughly compare our numerical results with analytical estimates, as was done in the previous section, for the reasons described already there. Also in comparing with previous estimates, we also discarded any possible additional contributions to the torque that might arise because of the turbulent magnetic fields. The detailed structure of the horseshoe region in the presence of turbulence and stratification deserves further study. Our results are summarized on Figure 17, where the torque dependence on planet mass is shown. For each simulation, we plot the last value of the cumulative average torque. Note however that only for part of the simulations the torque converges to a well defined value. It is possible to see a trend of the torque to reverse, corresponding to the addition of the contribution of the fully unsaturated horseshoe drag (). For the plot we assumed values for the width of the horseshoe region that are larger than the analytical estimate given by Paardekooper and Papaloizou (2009) and we use the value of the local surface density profile for the calculation of the torque. We see that the trend breaks down already for , where gap opening starts to become important and there is a transition into the Type II regime. Error bars represent the standard deviation of the time distribution of the torque. Note that since the raw torque is a highly oscillating quantity, the standard deviation does not match directly to the amplitude of the turbulent fluctuations, especially in the high-mass planet cases. For run 10, the standard deviation was found to be only 20% lower that in the turbulent run R9. In the Type II range, we plot the torque corresponding to the viscous timescale of the disk, taking . We find reasonable agreement with our simulation, taking into account the short simulation time, and that the value of the torque is still decreasing in the simulation after 100 orbits. Additionally we use the value of the initial, volume averaged , while the mid-plane value is smaller.

The question remains about the long term behavior of the torque, and whether this is only a transient behavior lasting for the first few hundred orbits (assuming the same local surface density profile), afterwards saturating and returning to standard negative Type I values. This is still a transient behavior in the sense that the planet can migrate out of the part of the disk where the local profile allows for outwards migration and enter a region where migration proceeds inwards again. Additionally it is limited by the lifetime of the pressure bumps, which we weren’t able to determine. We observe a stable pressure bump through the duration of our simulations. If there are other mechanisms such as the ones discussed in Masset et al. (2006b) that produces this type of locally increasing outwards density profile, then, in the presence of turbulence, these density bumps can also act as a protoplanet trap and halt, slow down or reverse inwards migration. Dzyurkevich et al. (2010) performed non-ideal MHD simulations of accretion disks with spatially varying resistivity. They also find zonal flows/pressure bumps not only at the snow line, e.g. a region with a jump in resistivity, but also inside the more active region. They already suggest that small planets should get trapped at those local pressure maxima (see also Kretke and Lin (2007)).

## 5 Summary

We studied the migration of planets under the influence of turbulence that is a result of the magneto-rotational instability. We find that, under the right conditions, planets can undergo systematic outwards migration in a locally isothermal disk. After long term averaging, transient or long term periods of outwards migration can help the survival and influence the mass accretion history of giant planet cores of a certain mass ratio. The contribution of the unsaturated horseshoe drag and the stochastic migration of low-mass planets, which are both consequences of the turbulence, should be incorporated into planet population synthesis models in order to test the influence of this element on the produced populations of planets. On future work we plan on studying low-mass planet migration in detail using similar stratified disk models.

Giant planets significantly decrease the magnetic stresses in the disk (mostly inside its orbit), effectively killing the turbulence, as we observe in our simulations. This is possibly a numerical effect and it will affect the accretion behavior of the disk and possibly the Type II migration rate of the giant planet. This issue deserves further study, with high resolution simulations to determine any possible effects of numerical dissipation of the magnetic fields induced by the presence of the planet. Additionally, in agreement with previous studies, we find that the gap opened by a planet in the presence of turbulence is wider than the gap produced in a quasi-laminar disk with an equivalent viscosity.

Name | Fixed orbit | Run time (local orbits) | |||
---|---|---|---|---|---|

R0 | Massless | 0.3 | 3.3 | Yes | 130 |

R1 | 5.0 | No | 85 | ||

R2 | 5.0 | No | 140 | ||

R3 | 3.3 | Yes | 89 | ||

R4 | 5.0 | No | 90 | ||

R5 | 5.0 | No | 100 | ||

R6 | 3.3 | Yes | 135 | ||

R7 | 4.0 | No | 85 | ||

R8 | 5.0 | No | 95 | ||

R9 | 5.0 | No | 100 | ||

R10 (HD) | 5.0 | No | 100 |

### Footnotes

- A more complete description of the type of model used in this paper and a detailed analysis of the MRI, magnetic fields and turbulent spectra can be found in Flock et al. (2011).
- We take the expression for an isothermal disk, in the zero gravitational softening limit.
- However, this expression for is derived using a 2D model of the HS region, which determines the viscous crossing time across the region, the libration time and the U-turn time.

### References

- Adams, F. C. and Bloch, A. M.: 2009, ApJ 701, 1381
- Alibert, Y., Mordasini, C., Benz, W., and Winisdoerffer, C.: 2005, A&A 434, 343
- Balbus, S. A. and Hawley, J. F.: 1991, ApJ 376, 214
- Balbus, S. A. and Hawley, J. F.: 1998, Reviews of Modern Physics 70, 1
- Baruteau, C. and Lin, D. N. C.: 2010, ApJ 709, 759
- Bate, M. R., Lubow, S. H., Ogilvie, G. I., and Miller, K. A.: 2003, MNRAS 341, 213
- Benz, W., Mordasini, C., Alibert, Y., and Naef, D.: 2008, Physica Scripta Volume T 130(1), 014022
- Bryden, G., Chen, X., Lin, D. N. C., Nelson, R. P., and Papaloizou, J. C. B.: 1999, ApJ 514, 344
- Crida, A., Morbidelli, A., and Masset, F.: 2006, Icarus 181, 587
- D’Angelo, G., Kley, W., and Henning, T.: 2003, ApJ 586, 540
- de Val-Borro, M., et al.: 2006, MNRAS 370, 529
- Dzyurkevich, N., Flock, M., Turner, N. J., Klahr, H., and Henning, T.: 2010, A&A 515, A70+
- Flock, M., Dzyurkevich, N., Klahr, H., and Mignone, A.: 2010, A&A 516, A26+
- Flock, M., Dzyurkevich, N., Klahr, H., Turner, N. J., and Henning, T.: 2011, ArXiv e-prints
- Fromang, S. and Nelson, R. P.: 2009, A&A 496, 597
- Fromang, S., Terquem, C., and Nelson, R. P.: 2005, MNRAS 363, 943
- Gardiner, T. A. and Stone, J. M.: 2005, Journal of Computational Physics 205, 509
- Goldreich, P. and Tremaine, S.: 1980, ApJ 241, 425
- Ida, S. and Lin, D. N. C.: 2004, ApJ 604, 388
- Ida, S. and Lin, D. N. C.: 2007, ArXiv e-prints
- Ida, S. and Lin, D. N. C.: 2008, in Y.-S. Sun, S. Ferraz-Mello, & J.-L. Zhou (ed.), IAU Symposium, Vol. 249 of IAU Symposium, pp 223–232
- Johansen, A., Youdin, A., and Klahr, H.: 2009, ApJ 697, 1269
- Johnson, E. T., Goodman, J., and Menou, K.: 2006, ApJ 647, 1413
- Kretke, K. A. and Lin, D. N. C.: 2007, ApJ 664, L55
- Laughlin, G., Steinacker, A., and Adams, F. C.: 2004, ApJ 608, 489
- Lyra, W., Johansen, A., Klahr, H., and Piskunov, N.: 2008, A&A 479, 883
- Masset, F. S.: 2002, A&A 387, 605
- Masset, F. S., D’Angelo, G., and Kley, W.: 2006a, ApJ 652, 730
- Masset, F. S., Morbidelli, A., Crida, A., and Ferreira, J.: 2006b, ApJ 642, 478
- Mignone, A., Bodo, G., Massaglia, S., Matsakos, T., Tesileanu, O., Zanni, C., and Ferrari, A.: 2007, ApJS 170, 228
- Mordasini, C., Alibert, Y., Benz, W., and Naef, D.: 2009, A&A 501, 1161
- Nelson, R. P.: 2005, A&A 443, 1067
- Nelson, R. P. and Gressel, O.: 2010, MNRAS 409, 639
- Nelson, R. P. and Papaloizou, J. C. B.: 2003, MNRAS 339, 993
- Nelson, R. P. and Papaloizou, J. C. B.: 2004, MNRAS 350, 849
- Nelson, R. P., Papaloizou, J. C. B., Masset, F., and Kley, W.: 2000, MNRAS 318, 18
- Ogihara, M., Ida, S., and Morbidelli, A.: 2007, Icarus 188, 522
- Oishi, J. S., Mac Low, M., and Menou, K.: 2007, ApJ 670, 805
- Paardekooper, S. and Papaloizou, J. C. B.: 2009, MNRAS 394, 2283
- Papaloizou, J. and Lin, D. N. C.: 1984, ApJ 285, 818
- Papaloizou, J. C. B. and Larwood, J. D.: 2000, MNRAS 315, 823
- Papaloizou, J. C. B. and Nelson, R. P.: 2003, MNRAS 339, 983
- Papaloizou, J. C. B., Nelson, R. P., and Snellgrove, M. D.: 2004, MNRAS 350, 829
- Tanaka, H., Takeuchi, T., and Ward, W. R.: 2002, ApJ 565, 1257
- Terquem, C. E. J. M. L. J.: 2003, MNRAS 341, 1157
- Thommes, E. W. and Murray, N.: 2006, ApJ 644, 1214
- Trilling, D. E., Lunine, J. I., and Benz, W.: 2002, A&A 394, 241
- Ward, W. R.: 1986, Icarus 67, 164
- Ward, W. R.: 1992, in S. F. Dermott, J. H. Hunter Jr., & R. E. Wilson (ed.), Astrophysical Disks, Vol. 675 of Annals of the New York Academy of Sciences, pp 314–323
- Ward, W. R.: 1997, Icarus 126, 261
- Winters, W. F., Balbus, S. A., and Hawley, J. F.: 2003, ApJ 589, 543
- Yang, C., Mac Low, M., and Menou, K.: 2009, ApJ 707, 1233