3 D Simulations of Type-I Migration in Nearly Laminar Disks

We report results from 3D hydrodynamic disk simulations with an embedded planet. We examine the planet’s migration dependence on the disk viscosity. We find that the Type I migration can be slowed down significantly in low viscosity disks. The density feedback effect that was studied in 2D simulations still exists in 3D disks and causes the slowing down of Type I. Future work will address further the quantitative differences between 3D and 2D in slowing down the migration.


Introduction
The standard Type I migration for low mass planets (Tanaka et al. 2002) presents a challenge for understanding planet formation.Several studies have suggested ways to slow down or even reverse this migration.These include the effects of magnetic fields (e.g., Terquem 2003) and the non-thermal effects (e.g., Paardekooper & Mellema 2006;Barteau & Masset 2008).In addition, Type I migration in low viscosity disks has been investigated and it was found that there exists a feedback effect that results in a slight mass redistribution (e.g., Rafikov 2002).This causes a reduction of the total torque on the planet, slowing down the migration.The effects of viscosity were investigated in Li et al. (2009) and Yu et al. (2010).Extensive 2D hydrodynamic simulations indeed show that the density feedback effect can slow down or even halt the Type I migration.It was found that the planet migration is strongly dependent on disk viscosities.One key finding is that there exists a critical planet mass above which the feedback effect can half the migration.
In this paper, we present new 3D hydrodynamic simulations of low-mass planet migration in low viscosity disks.We will first describe the numerical simulations in §2 and results are shown in §3.Conclusions are briefly summarized in §4.

Numerical Method and Initial Setup
The protoplanetary disk is treated as a 3D, self-gravitating gas whose motion is described by the locally isothermal Navier-Stokes equations in a spherical coordinate (r, θ, φ) centered on the star.The differential equations for the disk are similar to those given in Kley et al. (2009) with a different gravitational potential that is defined in Nelson et al. (2000).Simulations are carried out using a hydrodynamics code developed at Los Alamos (Li & Li, 2009, 2012) The equations are solved by directional split Godunov method for the inviscid Euler equations plus operator-split method for the EPJ Web of Conferences viscous source terms.We use a sub-cycling technique for the azimuthal sweep to alleviate the time step restriction.We also extend the FARGO scheme of Masset (2000) and modified in Li et al. (2001) to our 3D code to accelerate the transport in the azimuthal direction.Furthermore, we have implemented a reduced 2D (r, θ) and a fully 3D self-gravity solver on our uniform disk grid, which extends our 2D method (Li, Buoni, & Li 2008) to 3D.This solver uses a mode cut-off strategy and combines FFT in the azimuthal direction and direct summation in the radial and meridional direction.An initial axis-symmetric equilibrium disk is generated via iterations between the disk density profile and the 2D (r, θ) disk self-gravity.We do not need any softening in the disk self-gravity calculation as we have used a shifted grid method (Li et al. 2009) to calculate the potential.The motion of the planet is limited on the mid-plane and the equations are the same as given in D'Angelo et al. ( 2005), which we adapted to the polar coordinates with a fourth-order Runge-Kutta solver.The disk gravitational force on the planet is assumed to evolve linearly with time between two hydrodynamics time steps.The planetary potential acting on the disk is calculated accurately with a small softening given by a cubic-spline form (Kley et al. 2009).Since the torque is extremely sensitive to the position of the planet, we adopt the co-rotating frame that allows the planet moving only in radial direction if only one planet is present.This code has been extensively tested on a number of problems.For the earthmass planet with constant aspect ratio h = 0.05, the torque calculated using our code matches quite well with the 3D linear theory results by Tanaka et al. (2002).
The 3D disk is modeled between 0.4 ≤ r ≤ 2.0.In the vertical direction the annulus extends from the disk midplane (at θ = 90 • ) to about 7.5 • (or θ = 82.5 • ) above the midplane.The planet is initially located at r = 1, which corresponds to a physical distance of Jupiter's orbital radius (5.2AU), and orbits about a 1M star.Time is measured by the orbital period at r = 1, which is about 12 yr.A corotating frame that rotates with the initial angular velocity of the planet is used.The central star is fixed at the origin r = 0.The disk is assumed to be isothermal throughout the simulated region, having constant sound speed c s .The dimensionless disk thickness is scaled by the initial orbital radius of the planet h = c s /v φ (r = 1), where v φ is the Keplerian velocity.Without specification, h = 0.035 is used in the simulations.
The initial surface density is normalized to the minimum mass solar nebular model as Σ(r) = 152 f (r/5AU) (−3/2) g cm −2 , where f = 1 is used in our simulations.This surface density corresponds to an approximate 3D density profile ρ(r, θ) = ρ 0 (r) exp(−z 2 /2H 2 ), where the density in the midplane is ρ 0 (r) ∝ r −3 .The initial rotational profile of the disk is calculated so that the disk will be in equilibrium with the disk self-gravity and pressure (without the planet).A dimensionless, constant kinematic viscosity ν between 0 and 10 −6 is used, which corresponds to an effective Shakura and Sunyaev 0 ≤ α ≤ 4 × 10 −4 .(α = ν/h 2 ) Runs are made typically using a grid of (n r × n θ × n φ = 384 × 32 × 1536, though we have used higher resolution to ensure convergence on some runs.We remark that lower resolution (e.g., 192 × 16 × 768) may produce different results that are not be correct.

Results
We have performed one set of 3D runs where a single planet with mass m p = 10M ⊕ being placed in the same disk except with different disk viscosity ν = 10 −6 to 10 −9 , which corresponds to α ≈ 8.2 × 10 −7 to 8.2 × 10 −4 when c s = h = 0.035.We simulate the disk+planet system up to 2000 orbits to study its migration behavior.
Figure 1 shows the evolution of planet's semi-major axis for this set of runs.When the disk α is relatively large, the migration is essentially the same as Type I migration discussed in Tanaka et al. (2002).But when the disk viscosity is reduced, with α ≤ 8.2 × 10 −5 , the migration is slowed down significantly.This result is similar to that obtained in 2D simulations (Li et al. 2009).Careful analyses  We also made another set of runs where we explore the effects of numerical resolution.Figure 2 shows the comparison of a 10M ⊕ planet migrating in a disk with f = 1, c s = 0.035 and ν = 10 8 .The simulations have three different resolutions and they indicate somewhat different rates of migration.For highest resolution, the simulation was only carried out to ∼ 900 orbits.

05003-p.3
Instabilities and Structures in Proto-Planetary Disks EPJ Web of Conferences

Summary
We have performed 3D hydrodynamic simulations of disks with low mass planets where we investigate the density feedback effects originally studied in the 2D simulations.Such feedback effect was shown to be important in low viscosity disks to slow down or even stop the Type I migration in 2D.The 3D simulations qualitatively confirm the conclusion that the feedback effect still exists in 3D low viscosity disks.We also find that it is important to use enough resolutions to study these effects, though it is not clear whether our highest resolution simulations have reached convergence.Future studies will address quantitatively whether and how the feedback effects in 3D differ from that in 2D and whether this is still a mechanism for halting Type I migration in 3D disks.

Figure 1 .
Figure1.Influence of disk viscosity on planet migration.The planet mass is 10M ⊕ in a disk with h = 0.035 and f = 1.The vertical axis is the orbital semi-major axis of the planet in units of the initial orbital radius.The horizontal axis is the time in units of the initial planet orbital period.For dimensionless kinematic viscosity ν = 10 6 , the planet undergoes the typical Type I migration, but its migration slows down significantly when ν is smaller.

Figure 2 .
Figure 2. Migration of a 10M ⊕ planet in a low viscosity disk.Three curves correspond to three different spatial resolutions, indicated by {n r , n θ , n φ }.Even with relatively high resolution, it is not clear that the simulations have reached convergence.