Local Linear Analysis of Interaction between a Planet and Viscous Disk and an Implication on Type I Planetary Migration
aa r X i v : . [ a s t r o - ph . E P ] J un ApJ Accepted
Local Linear Analysis of Interaction between a Planet andViscous Disk and an Implication on Type I Planetary Migration
Takayuki Muto Department of Physics, Kyoto University,Kitashirakawa-oiwake-cho, Sakyo-ku, Kyoto, 606-8502, Japan andShu-ichiro Inutsuka
Department of Physics, Nagoya University,Furo-cho, Chikusa-ku, Nagoya, 464-8602, Japan [email protected]
ABSTRACT
We investigate the effects of viscosity on disk-planet interaction and discusshow type I migration of planets is modified. We have performed a linear calcula-tion using shearing-sheet approximation and obtained the detailed, high resolu-tion density structure around the planet embedded in a viscous disk with a widerange of viscous coefficients. We use a time-dependent formalism that is useful ininvestigating the effects of various physical processes on disk-planet interaction.We find that the density structure in the vicinity of the planet is modified andthe main contribution to the torque comes from this region, in contrast to invis-cid case. Although it is not possible to derive total torque acting on the planetwithin the shearing-sheet approximation, the one-sided torque can be very dif-ferent from the inviscid case, depending on the Reynolds number. This effecthas been neglected so far but our results indicate that the interaction betweena viscous disk and a planet can be qualitatively different from an inviscid caseand the details of the density structure in the vicinity of the planet is criticallyimportant.
Subject headings: planet and satellites: formation — solar system: formation JSPS Research Fellow
1. Introduction
Disk-planet interaction and type I planetary migration is one of the main issues in thetheory of planet formation. A low mass planet or the core of the gas giant planet embeddedin the protoplanetary disk excites spiral density wave in the disk, which carries angularmomentum flux away from the planet. The linear calculation of the spiral density wave isfirst carried out by Goldreich and Tremaine (1979), and there has been a number of extensivework assuming inviscid, isothermal disk such as Ward (1986, 1989), Artymowicz (1993), andTanaka et al. (2002). In the case of an inviscid, isothermal disk, protoplanets migratetowards the central star as a result of disk-planet interaction very rapidly, thereby imposinga serious problem on the theory of planet formation.The decay of the orbital semi-major axis of the protoplanet is a result of the subtledifference between the inner disk (the part of the disk inside the orbit of the protoplanet)and the outer disk (the part of the disk outside the orbit of the protoplanet). The outerdisk exerts torque on the protoplanet in such a way that the semi-major axis decreases,and the inner disk exerts opposite torque. In the lowest order of the local approximation, orshearing-sheet approximation, these two effects cancel each other, and the difference appearsat higher order. Therefore, type I migration can be very sensitive to the physical state of thedisk. Indeed, when non-isothermal effect is taken into account, the direction of migrationcan be reversed (Paardekooper and Mellema 2006, Baruteau and Masset 2008, Paardekooperand Papaloizou 2008). It is also discussed that the global magnetic field can slow down therate of type I migration or even reverse the direction (Terquem 2003, Fromang et al. 2005,Muto et al. 2008).In order to construct a predictable, self-consistent model for planet formation, it isimportant to make clear how different physical processes of the disk affect disk-planet inter-action and type I migration. In this paper, we revisit the effects of viscosity on disk-planetinteraction. There has been a number of work on this topic. Meyer-Vernet and Sicardy(1987) discussed that the viscosity does not affect the torque exerted on the planet as longas its value is small, although the shape of the wake may be modified. Papaloizou and Lin(1984, and a series of their paper) and Takeuchi et al. (1996) calculated that the damp-ing of the density wave and the transfer of the angular momentum from the density waveto the background disk, resulting in the gap formation around a massive planet. Masset(2001) investigated the effects of viscosity on the horseshoe regions and discussed how thecorotation torque may be altered by the effect of viscosity. More recently, Paardekooper andPapaloizou (2009b) investigated the dynamics of corotation region using both linear andnon-linear calculations.Although there has been a number of investigations on the disk-planet interaction in 3 –a viscous disk (for example, Masset 2001, 2002 or Paardekooper and Papaloizou 2009a forcorotation torque, D’Angelo et al. 2002, 2003 for high resolution numerical study), therehas not yet been a study of wide range in the viscous coefficient that requires an analysison the detailed density structure in the vicinity of the planet. In this paper, as a first stepfor the complete investigation, we show results of linear calculation in a local, shearing-sheetanalysis. We use time-dependent methods to calculate the disk response against the planetpotential. Although this method has not been widely used so far, this is useful in investigatingthe effects of various physical processes on disk-planet interaction. We calculate the non-axisymmetric density structure around the planet and investigate how the resulting torqueis altered by the effect of shear viscosity. We have studied wide range of the parametersof viscous coefficient and calculate the density structure with high resolution. We find thatthe density structure in the vicinity of the planet is altered in a viscous disk, with viscouscoefficient of ∼ . α (standard α parameter, see equation (27) for definition),which may be realized as a result of the turbulence induced by magneto-rotational instability(MRI, Balbus and Hawley 1991, Sano et al. 2004).Since it is not possible to calculate the total torque (differential torque between thedisk inside the orbit of the planet and outside) in shearing-sheet approximation, we look atone-sided torque to investigate how type I migration rate can be affected qualitatively byviscosity. One-sided torque is the torque exerted by one side of the disk with respect to theplanet’s orbital radius. We also note that shearing-sheet calculation is a useful first stepbefore differential torque is calculated (Tanaka et al. 2002). We find that as we increasethe values of viscous coefficient α , one-sided torque stays unchanged until α ∼ < .
01, thenincreases until α ∼
1, and finally decreases. The torque can be factor of two larger thaninviscid case when α = 0 .
1, and more importantly, the enhancement of the torque is aresult of the modified density structure in the vicinity of the planet, which has not yet beeninvestigated in detail. Our results indicate that the physical mechanisms of the disk-planetinteraction in a viscous disk may depend on the detailed density structure around the planet.The plan of this paper is as follows. In section 2 we show the basic equations anddescribe the time-dependent methods to obtain a stationary, non-axisymmetric pattern oflinear perturbation. We also discuss how two-dimensional (2D) and three-dimensional (3D)analyses are related in this section. In section 3, we show the results of 2D and 3D calculationsand show how the torque behaves as we vary the amount of viscosity. In section 4, weexplain the qualitative behavior of the torque using a simple analytic model. Section 5 is fordiscussion and summary. 4 –
2. Basic Equations and Numerical Methods2.1. Linear Perturbation Analysis
We consider isothermal Navier-Stokes equations with one planet using shearing-sheetapproximation. ∂ρ∂t + ∇ · ( ρ v ) = 0 (1) ∂ v ∂t + v · ∇ v + 2Ω p e z × v = − c ρ ∇ ρ + 3Ω x e x + ν ∇ v + 13 ν ∇ ( ∇ · v ) − ∇ ψ p (2)where ρ is density, v is velocity, Ω p is the Kepler angular velocity of the planet, ψ p is thegravitational potential of the planet, ν is shear viscosity. The third term of left hand sideof equation (2) is Coriolis force and the second term of right hand side of equation (2) istidal force. We assume the planet is located at the origin and stationary with respect to thiscoordinate system, i.e., ψ p = ψ p ( x, y, z ) , (3)where ψ p does not depend on time. In this paper, we neglect bulk viscosity for simplicity.We also neglect vertical stratification and assume that the background density without aplanet is homogeneous. This gives an uncertainty in the box size in the z -direction, whichwill be discussed in Section 2.2. We also neglect the effect of global gas pressure gradientexerted on the background disk and assume that the background gas is rotating at Keplervelocity.The unperturbed state without a planet is given by ρ = ρ = const (4) v = −
32 Ω p x e y . (5)In the presence of viscosity, if we perform a global analysis, there is a mass accretion ontothe central star in general. However, in the shearing-sheet approximation, where linearbackground shear is assumed, this effect is not taken into account. We expect that thedensity structure only in the vicinity of the planet can be well approximated even if weneglect the global mass accretion.We consider linear perturbation. All the perturbation quantities are denoted with δ ,e.g., δρ for density perturbation. Perturbation equations are given by (cid:18) ∂∂t −
32 Ω p x ∂∂y (cid:19) δρρ + ∇ · δ v = 0 (6) 5 – (cid:18) ∂∂t −
32 Ω p x ∂∂y (cid:19) δ v − p δv y e x + 12 Ω p δv x e y = − c ∇ δρρ + ν ∇ δ v + 13 ν ∇ ( ∇ · δ v ) − ∇ ψ p (7)We solve equations (6) and (7) to obtain a steady state solution and calculate torque exertedon one side (either x > x <
0) of the disk by the planet. The torque exerted on theplanet is obtained as a backreaction of this torque. Since we use shearing-sheet approxi-mation, torque exerted from each side of the planet is the same in magnitude and oppositein sign. Although we do not obtain a net torque, it is still possible to investigate how thedisk structure is affected by the viscosity and qualitatively predict how the disk and planetinteract.In order to obtain a steady state solution, we solve linear perturbation equations (6) and(7) using Fourier transform methods given by Goodman and Rafikov (2001). We transformthe equations into the shearing coordinate ( t ′ , x ′ , y ′ , z ′ ) defined by t ′ = t (8) x ′ = x (9) y ′ = y + 32 Ω p xt (10) z ′ = z. (11)In this coordinate system, temporal and spatial derivatives are given by ∂∂t = ∂∂t ′ + 32 Ω p x ′ ∂∂y ′ , (12) ∂∂x = ∂∂x ′ + 32 Ω p t ′ ∂∂y ′ , (13) ∂∂y = ∂∂y ′ , ∂∂z = ∂∂z ′ . (14)Then, perturbation equations become ∂∂t ′ δρρ + ∇ · δ v = 0 (15) ∂∂t ′ δ v − p δv y e x + 12 Ω p δv x e y = − c ∇ δρρ + ν ∇ δ v + 13 ν ∇ ( ∇ · δ v ) − ∇ ψ p , (16) 6 –where the spatial derivatives in ∇ are given by equations (13) and (14).The coefficients of equation (15) and (16) are now independent of ( x ′ , y ′ , z ′ ). Therefore,if we Fourier transform in spatial directions, we obtain a set of ordinary differential equationsdecoupled for each ( k ′ x , k ′ y , k ′ z ) mode (Goldreich and Lynden-Bell 1965), δf ( t ′ , x ′ , y ′ , z ′ ) = X δf ( t ′ , k ′ x , k ′ y , k ′ z ) exp (cid:2) i ( k ′ x x ′ + k ′ y y ′ + k ′ z z ′ ) (cid:3) . (17)The relationship of wavenumber in ( x, y, z )-coordinate and ( x ′ , y ′ , z ′ )-coordinate is given by k x ( t ) = k ′ x + 32 Ω p k y t. (18) k y = k ′ y , k z = k ′ z (19)Equation (18) indicates that the value of radial wavenumber in ( x, y, z )-coordinate evolveswith time owing to the background shear. The value of radial wavenumber in shearingcoordinate k ′ x gives the initial value of radial wavenumber in ( x, y, z )-plane. Equations ofcontinuity (6) and motion (7) are now ddt δρρ + ik x ( t ) δv x + ik y δv y + ik z δv z = 0 , (20) ddt δv x − p δv y = − c ik x ( t ) δρρ − ν ( k x ( t ) + k y + k z ) δv x − νk x ( t )( k x ( t ) δv x + k y δv y + k z δv z ) − ik x ( t ) ψ p , (21) ddt δv y + 12 Ω p δv x = − c ik y δρρ − ν ( k x ( t ) + k y + k z ) δv y − νk y ( k x ( t ) δv x + k y δv y + k z δv z ) − ik y ψ p , (22)and ddt δv z = − c ik z δρρ − ν ( k x ( t ) + k y + k z ) δv y − νk y ( k x ( t ) δv x + k y δv y + k z δv z ) − ik z ψ p . (23)In equations (20)-(23), we write all the terms in ( t, x, y, z )-coordinate. These equationsdescribe the excitation and evolution of density wave by the source ψ p . The source of thedensity wave becomes zero when wavenumber ( k ′ x , k ′ y , k ′ z ) is large. 7 –As an initial condition, we assume that there is no perturbation: δρ ( t = 0) = δ v ( t =0) = 0, and k ′ x is taken to be sufficiently large in absolute magnitude and negative (positive)in sign for positive (negative) k y . If we use spatial resolution of x -direction ∆ x , we shouldtake k ′ x = ∓ π/ ∆ x , where upper (lower) sign is used for positive (negative) k y . As a result oftime evolution, k x ( t ) increases (decreases) for each positive (negative) k y mode. When k x ( t )reaches ± π/ ∆ x , where upper (lower) sign is for positive (negative) k y modes, we obtain asteady state profile of perturbation quantities in Fourier space for non-axisymmetric ( k y = 0)modes. The Fourier amplitude of each k x mode in steady state is given by the time evolutiondata through equation (18). The profile of physical quantities in Fourier space is then inverseFourier transformed to obtain values in real space. Since we are interested in torque, we donot calculate axisymmetric modes ( k y = 0).The standard procedure to obtain steady state solution is as follows. Firstly, stationarysolution in the frame corotating with the planet is assumed: ∂/∂t = 0. Then, Fouriertransformation in the y - and z -directions is performed to obtain the ordinary differentialequations in the x -direction. Finally, these ordinary differential equations are solved byimposing outgoing boundary condition, or equivalently the boundary condition that admitsonly trailing wave, see e.g., Goldreich and Tremaine (1979), Korycansky and Pollack (1993),or Tanaka et al. (2002).In the presence of viscosity, this method introduces higher order derivative with respectto x , since viscous terms include the second order derivative in the radial direction. Theresulting ordinary differential equations are higher order in the x -derivative than in inviscidcases. The terms with the highest order derivative comes from viscous terms and theircoefficients are viscous coefficient ν . In this case, it is difficult to take natural ν → ν → The reason why we obtain a steady state as a result of time evolution is as follows. Let us assumethat if | k x | > K c , the source ψ p is so small that we can approximate it to be zero. The Fourier amplitudeof the specific k x mode in ( t, x, y, z )-coordinate at time t = t is given by the solution at t = t with k ′ x = k x − (3 / p k y t mode. Therefore, the absolute magnitude of k ′ x that contributes to the Fourieramplitude of fixed k x becomes large as time t increases. If the absolute magnitude of k ′ x is larger than K c ,there is no wave excitation for this mode until | k x ( t ) | becomes smaller than K c . Therefore, when | k ′ x | > K c ,the contribution to the k x mode is always the same regardless of t . x -direction depending on the modesspecified by k y and k z . Details of our methods are given below.We also note that if we assume an initial condition with non-zero perturbation, addi-tional homogeneous wave is introduced. This wave has both leading ( k x k y <
0) and trailing( k x k y >
0) components and the resulting solution is simply the superposition of the specificsolution of Equations (20)-(23) assuming zero initial condition and homogeneous ( ψ p = 0)solution assuming the specified initial condition. The solution with non-zero initial conditionis irrelevant in the present problem because we consider the perturbation that is induced bythe gravitational perturbation by the planet.We write the box size in ( x, y, z )-directions by ( L x , L y , L z ) and the coordinate systemextends − L x / < x < L x / k ′ x = ∓ π/ ∆ x for positive (negative) k y modes.2. Equations (20)-(23) are solved with initial condition δρ ( t = 0) = δ v ( t = 0) = 0.3. Resulting solutions in ( k x , k y , k z ) space are inverse Fourier transformed to real space.In steps 2 and 3, there are some points we need to take care. In calculating the Fouriermodes, we have varied the mesh number in k x directions, which corresponds to the stepsize of time, in such a way that all the oscillations are well resolved for each ( k y , k z ) mode.The mesh number in k x direction becomes as much as 10 for small k y . We perform inverseFourier transform in k x direction using the increased number of mesh in k x direction. Thisprocedure effectively makes the box size of x -direction longer for a specific k y mode, andthe data corresponding to necessary x , which is − L x / < x < L x /
2, is then interpolatedfrom the resulting profile in ( x, k y , k z ) space. This data is then used to perform inverseFourier transform in ( k y , k z ) directions. In this way, it is possible to avoid aliasing effectin x -direction and outgoing boundary condition in the steady state solution is effectivelysatisfied. We do not use any window function, which has been incorporated by Goodmanand Rafikov (2001), in Fourier transform. We have checked that the resulting torque is notaffected by the window function in the absence of viscosity. Also, since window functionintroduces additional artificial dissipation, this affects our results with viscosity.Once we have obtained the profiles of perturbation quantities in real space, we can 9 –calculate the one-sided torque exerted on x > T = − r p Z L x Z L y / − L y / Z L z / − L z / dxdydzδρ ( x, y, z ) ∂ψ p ∂y , (24)where r p is the semi-major axis of the planet, and the torque exerted on the planet isobtained as backreaction. For later convenience, we define “torque distribution” T ( x, y, z )that is simply the torque exerted on the fluid element located at ( x, y, z ), T ( x, y, z ) = − r p δρ ( x, y, z ) ∂ψ p ∂y (25)and “torque density” that is the torque exerted on an annulus of the disk at x , T ( x ) = − r p Z L y / − L y / Z L z / − L z / dzdyδρ ( x, y, z ) ∂ψ p ∂y (26)We normalize equations using c , Ω p , and ρ so that the homogeneous equations ( ψ p =0)contain only one dimensionless variable α that is defined by α = νcH , (27)where H = c/ Ω p is the scale height of the disk. Since we consider the linear perturbationanalysis, the amplitude of the perturbation is proportional to the normalized planet mass: µ = GM p /Hc . (28)We show the results with µ = 1 in subsequent sections.We use the fifth order Runge-Kutta method to solve time evolution and FFT routinegiven by Press et al. (1992) to perform Fourier transform. We use box size of L x = 10 H and L y = 40 H , and the grid number in x - and y -directions ( N x , N y ) = (512 , z -direction, L z , which is discussed in thenext subsection. In this subsection, we discuss how 2D and 3D calculations differ from each other and de-termine the relevant box size in the z -direction with which the effect of vertical stratificationis effectively taken into account. 10 –We compare 2D calculation and 2D mode of 3D calculation. By “2D calculation”, wemean the gravitational potential of the planet is given by ψ p , ( x, y ) = − GM p p x + y + ǫ , (29)where ǫ denotes the softening length of 2D calculation, which is usually incorporated inmost of 2D work. By “2D mode of 3D calculation”, we mean that the analysis is restrictedto k z = 0 mode of 3D calculation. Three-dimensional potential is given by ψ p , ( x, y, z ) = − GM p p x + y + z + ǫ . (30)Therefore, gravitational potential for “2D mode of 3D calculation” is given by ψ p , ( x, y ) = − L z Z L z / − L z / dz GM p p x + y + z + ǫ = − GM p L z log (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x + y + ǫ + L z / L z p L z / x + y + ǫ x + y + ǫ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (31)where we denote softening length in 3D calculation by ǫ . We note that ψ p , and ψ p , coincides if x + y + ǫ ≫ L z , but they differ if x + y + ǫ ≪ L z , since in this case,3D potential behaves as ψ p , ∼ log r while 2D potential behaves as ψ p , ∼ /r , where r = x + y + ǫ . Vertical averaging makes the potential weaker in the vicinity of the planet,thereby producing small perturbation in 2D mode of 3D calculation. Figure 1 shows thetorque obtained by 2D inviscid ( ν = 0) calculation using ψ p , and ψ p , with various L z .We note that in producing Figure 1, one-sided torque is calculated according to equation(24) and the value is normalized byΓ = µ ρ L z r p Hc . (32)It is clear that 2D and 3D results coincide if L z is small, but the 3D result becomes smallwhen L z is large.The reasonable value of L z may be obtained by comparing the results of 2D mode of 3Dcalculation with calculation in which vertical stratification is taken into account. Tanaka etal. (2002) performed such calculations for an inviscid disk. Their method is to decomposevertical modes in terms of Hermite polynomials, and they have found that the verticallyaveraged 2D mode gives most of the contribution to the resulting torque. Their 2D mode ofthe potential, ψ p , TTW , is given by ψ p , TTW = − √ πH Z ∞−∞ dz GM p p x + y + z + ǫ e − z / H = − GM p √ πH exp (cid:2) ( x + y + ǫ ) / (cid:3) K (cid:2)(cid:0) x + y + ǫ (cid:1) / (cid:3) , (33) 11 –where K ( x ) is the modified Bessel function of the zeroth order. Fitting equation (31) withequation (33), we find L z = 2 . H gives a reasonable agreement, see Figure 2. Zeroth orderof modified local approximation incorporated by Tanaka et al. (2002) coincides with theshearing-sheet approximation, and we have exactly the same homogeneous equations withthem for 2D ( k z = 0) mode if we assume ∂/∂t = 0. As is clear from Figure 1, 2D modein 3D calculation gives smaller amount of torque than the calculation using 2D potentialgiven by equation (29). This explains why Tanaka et al. (2002) has obtained slightly smalleramount of torque compared to 2D results. We have also checked that setting L z = 2 . H ,our results of one-sided torque agree with those determined by using Tanaka et al. (2002)methods within an error of 3%. We use L z = 2 . H and the mesh number in z -direction is taken to be 128. In to-tal, we have performed calculation with ( L x , L y , L z ) = (10 H, H, . H ) with mesh number(512 , , x, ∆ y, ∆ z ) = (0 . H, . H, . H ). How-ever, the radial box size is variably extended according to modes in calculation, as discussedin Section 2.1. Effective box size in x -direction is as much as 10 H . We use ǫ = 10 − H for the softening parameter. Our results vary upto 30% for large values of viscosity whensmoothing length is varied upto twice the grid resolution. Therefore, our results give at leasta qualitative view of how disk-planet interaction is altered by the effects of viscosity. Thevariation of one-sided torque as a function of smoothing parameter is further discussed inSection 3.2.
3. Results
In this section, we show our results of density structure and one-sided torque obtainedby linear analysis. In this section, the normalization of the torque is given byΓ = µ ρ r p H c . (34)Note the difference in the normalization from that used in Section 2.2. The results of 2Dmode torque given in Section 3.1 are different by a factor of L z /H = 2 . There is another complication regarding the normalization of the torque when comparing the value ofthe torque obtained by our calculation and that given by Tanaka et al. (2002). They use surface densityto normalize the torque they have obtained. We use normalization given by equation (32). Correspondencebetween the two results is given by reading our ρ L z to σ p of Tanaka et al. (2002), and this is how we haveobtained the agreement with their calculation.
12 –Figure 1. The normalization of torque distribution defined by equation (25) is taken to be T ( x, y, z ) µ ρ r p c H − , (35)and torque density T ( x ) defined by equation (26) is given by T ( x ) µ ρ r p c . (36) We first show the results of calculations restricted to 2D modes ( k z = 0). Although thisis only an approximation, physics involved is made clear. Figure 3 shows the torque obtainedfor various viscosity parameter α . For small viscous coefficients ( α ∼ < − ), torque is notaffected by the viscosity, as discussed by Meyer-Vernet and Sicardy (1987). However, whenviscous coefficient is large, it is shown that one-sided torque increases and peaks around α ∼ However, our resultshows that one-sided torque peaks at α ∼
1. This result originates from two different effectsof viscosity. First, viscosity damps spiral density wave, as discussed by Papaloizou andLin (1984) or Takeuchi et al. (1996). The second effect of viscosity may be observed bylooking at the density structure in the vicinity of the planet. Figure 4 shows the contour ofdensity profile of the xy -plane for α = 10 − and α = 10 − . It is clear that the oval densityprofile in the vicinity of the planet is slightly tilted for large viscous coefficient, while thespiral density wave is damped. This tilt produces asymmetry of density structure in the y -direction, thereby exerting one-sided torque on the planet.Figure 5 shows the torque density for α = 10 − and α = 10 − , respectively. It isclear that the peak of torque density locates slightly closer to the planet for α = 10 − than α = 10 − . The main contribution to the torque comes from the flow structure in the vicinityof the planet, which is evident when we consider the asymmetry of the density perturbation The effect of viscosity on the fluid elements trapped in the horseshoe regions enhances the resultingtorque since viscosity keeps the asymmetry of the potential vorticity. This applies corotation torque, seeMasset (2001) for detail. In shearing-sheet calculations presented here, however, there is no corotation torquesince background values are assumed to be constant.
13 –between y > y < y > y < T ( x, y ) + T ( x, − y ) , (37)as a function of x and y ( > y -directions. In Figure 6, we compare the results for two different values of viscosity, α = 10 − and α = 10 − . In the case of α = 10 − , the forward-back difference of the torque in thevicinity of the planet almost vanishes since the density perturbation around the planet issymmetric in the y -direction. In the case of α = 10 − , however, the most of the contri-bution to the torque comes from the region very close to the planet since the oval densityperturbation structure is tilted in y < In this section, we present the results of 3D calculation. Figure 7 shows the one-sidedtorque as a function of viscosity. The qualitative behavior of torque increasing with viscosityis unchanged, but the enhancement of torque becomes large compared to 2D calculation.We fit the data and obtain the following empirical relation Tρ r p H c = (0 .
94 + 10 α ) e − . α (38)for α < . L z = 2 . H . We have chosen a form of fitting function in such a way thattorque converges to a non-zero value for α ≪
1, peaks at α ∼
1, and decreases to zero for α ≫
1. This fitting formula is in reasonably good agreement with our calculation in thisrange of viscosity coefficient.Just as in the calculations restricted to 2D modes, the density structure mainly in thevicinity of the planet contributes to the one-sided torque if large values of viscosity is as-sumed. Figure 8 shows the torque density profile obtained for different viscosity coefficients.The more viscous is the disk, the closer to the planet is the location of the dominant contri-bution to the torque. Figure 9 compares the density structure in the yz -plane at x = 0 . H for calculations with α = 10 − and α = 10 − . The asymmetry in the y -direction in stronglyperturbed region is present in calculations with large viscosity.In 3D calculations, gravitational potential in the vicinity of the planet is not as stronglysoftened as in 2D modes. In the vicinity of the planet, gas feels the gravitational potentialthat decreases as − /r , where r is the distance from the planet, in full 3D calculations.However, the gravitational potential varies only logarithmically when r ∼ < H (see equation(31)) in calculations restricted to 2D modes. Since the density fluctuation around the planet 14 –is given by δρ/ρ ∼ ψ p /c , the deeper gravitational potential gives the higher value of thetorque exerted on the planet if there is a substantial asymmetry in the density structure.Calculations restricted to two-dimensional modes predict one-sided torque very well ifviscosity parameter is small, since regions close to the planet ( r ∼ < H ) is not very important.However, if large values of viscosity is used, it is necessary to perform full three-dimensionalcalculations since the density structure close to the planet is important.At the end of Section 2.2, we addressed that varying smoothing length changes one-sidedtorque upto 30%. Below, we discuss the effects of smoothing length in our calculation andargue that the values of one-sided torque obtained for large viscosity, α ∼ > .
01, seems to bethe lower limit of the one-sided torque.Figure 10 shows how one-sided torque of three-dimensional calculation varies withsmoothing length ǫ . Results with α = 10 − and α = 10 − are shown. In this calcula-tion, we used the box size with ( L x , L y , L z ) = (10 H, H, . H ) and the grid number with( N x , N y , N z ) = (512 , , y -direction is better by a factorof four compared to the parameters used in Figure 7. Grid resolution is ∼ . H in alldirections in Figure 10.It is shown that one-sided torque converges well for small values of viscosity parameter.If α is as large as 0 .
1, the calculation converges for sub-grid smoothing lengths, and resultsvary approximately 30% if smoothing length with twice the grid scale ( ∼ . H ) is assumed.The qualitative behavior of one-sided torque can be understood if we notice that smallersmoothing length gives a deeper potential in the vicinity of the planet, and if large values ofviscosity is assumed, density structure in the vicinity of the planet is important for one-sidedtorque.Since contribution to one-sided torque in case of large viscosity mainly comes fromthe density structure close to the planet, it may depend on the grid resolution used inthe calculation. Figure 11 shows the one-sided torque obtained for different resolutions.Calculations with ( N x , N y , N z ) = (512 , , , , , ,
32) are shownwhile keeping ( L x , L y , L z ) = (10 H, H, . H ). These values corresponds to grid resolutionswith ∆ x = ∆ y = ∆ z = 0 . H , 0 . H , and 0 . H , respectively. The value of softeningparameter is kept ǫ = 10 − H for all the calculations. It is shown that the value of one-sided torque for α = 10 − is well converged while the value of the torque for α = 10 − variesapproximately 30% when grid resolution is varied by a factor of four. We note that in case ofsmall viscosity, we have obtained well-converged values of one-sided torque since the effectiveLindblad resonances, which are located at | x | ∼ > (2 / H , are all well resolved.We note that the vertical averaging, large softening length, and coarser grid all introduce 15 –more softened potential in the vicinity of the planet and therefore one-sided torque becomessmaller in case of large viscosity. Equations (30) and (31) show that the vertically averagedpotential diverges at the location of the planet more mildly. Gravitational force becomesweaker in the vicinity of the planet if we use larger values of softening length. The largergrid scale cuts off gravitational potential at larger distance.From the behaviors when we vary the softening parameter and the grid scale, we con-clude that the values of the one-sided torque obtained in our calculation for large viscosityis the lower limit. We also argue that the numerical treatment in the vicinity of the planetmay cause a quantitative difference in the one-sided torque since it may change the densitystructure in the vicinity of the planet. If high values of viscosity is assumed, the densitystructure close to the planet is important and therefore, three-dimensional calculation withhigh resolution and small softening parameter is essential. The resolution and the softeningparameter should be determined, in principle, by considering the realistic size of the planet.
4. Analytic Treatment of Density Structure around the Planet
We have seen that the viscosity exerted on the disk can change the density structurein the vicinity of the planet and in a viscous disk, the planet experiences one-sided torquefrom the gas well inside the effective Lindblad resonance. In this section, we show that adissipative force distorts the density structure close to the planet by using a simple analyticmodel.We consider a two-dimensional model with a friction force exerted only in the y -direction.We consider equations (20)-(22) with viscosity terms replaced by friction. We assume k x ( t ) ∼ k z = 0 and consider 2D perturbation for simplicity. The set of equations wesolve is dσdt + ik y δv y = 0 , (39) ddt δv x − p δv y = 0 , (40)and ddt δv y + 12 Ω p δv x = − c ik y (cid:18) ψ p c + σ (cid:19) − γδv y , (41)where σ = δρ/ρ and γ is a drag coefficient.Using equations (39) and (40) to eliminate δv y and integrating once the resulting equa- 16 –tion, we obtain δv x = − p ik y σ + V, (42)where V is a constant of integration. Physically, V is the perturbation of vortensity, whichis actually conserved in this model since we neglect terms with k x ( t ) and we assume frictionforce is exerted only in the y -direction. Therefore, V is zero if we assume there is no vortensityperturbation initially. Since we are interested in the perturbation that arises from the planetpotential, we assume V = 0 hereafter.Using equations (39) and (42) to eliminate δv x and δv y from equation (41), we obtain asingle telegraph equation with a source terms, d σdt + γ dσdt + s σ = − k y ψ p , (43)where s = c k y + Ω . (44)The right hand side of equation (43) is the source of perturbation that is induced bythe planet’s gravity. The source potential ψ ( k x ( t ) , k y ) depends on time through the time-dependence of k x ( t ).Let us now derive the steady state profile of density perturbation in ( t, x, y )-coordinatespace. In order to solve equation (43), we Fourier transform perturbation in t -direction f ( t, k y ) = Z ∞−∞ dsf ( s, k y ) e ist , (45)where f denotes any of perturbation quantities. We note that the Fourier transformation in t -direction is equivalent to the inverse Fourier transformation of k x modes into x -coordinatespace. The “frequency”, s , in equation (45) is the frequency of perturbation experienced bya single mode specified by k y .The real space quantity f ( x, y ) is obtained by f ( x, y ) = Z ∞−∞ dk y Z ∞−∞ dk x f ( k x , k y ) e i ( k x x + k y y ) . (46)In a steady state, the value of k ′ x can be taken arbitrary [see discussion in Section 2.1].Therefore, we can take k ′ x = 0 without loss of generality. Time-dependence of k x is given by k x ( t ) = 32 Ω p k y t, (47)and therefore, we obtain dk x = 32 Ω p k y dt. (48) 17 –Changing the integration variable from k x to t in equation (46), we obtain f ( x, y ) = Z ∞ dk y Z ∞−∞ dt
32 Ω p k y f ( t, k y ) exp (cid:20) k y x Ω p t + k y y (cid:21) − Z −∞ dk y Z ∞−∞ dt
32 Ω p k y f ( t, k y ) exp (cid:20) k y x Ω p t + k y y (cid:21) . (49)The value f ( t, k y ) appeared in this equation is Fourier transformed in the t -direction accord-ing to equation (45). Substituting equation (49) into (45), performing integral in t and s ,we obtain the relationship between f ( s, k y ) and f ( x, y ), f ( x, y ) = 2 π Z ∞ dk y
32 Ω p k y × (cid:2) e ik y y f ( s = − (3 / p k y x, k y ) + e − ik y y f ( s = (3 / p k y x, − k y ) (cid:3) . (50)One-sided torque is calculated by equation (24) (but there is no integral in the z -directionin two-dimensional problem considered in this section). We choose the phase of the potentialsuch that ψ p ( s, k y ) is real. It is possible because the potential of the planet is sphericallysymmetric and the planet is fixed at the origin of the coordinate system. Thus, we obtain T ∝ Z ∞ dk y k y Z ∞ ds Im [ σ ( s, − k y )] ψ p ( s, − k y ) , (51)where constant of proportionality does not depend on k y , γ , and x and “Im” denotes theimaginary part. From this equation, torque is exerted on the planet only when σ ( s, k y ) and ψ p ( s, k y ) are out of phase. In other words, the amount of torque can be estimated by lookingat the imaginary part of σ ( s, k y ), provided that the phase of Fourier transform is taken insuch a way that ψ p ( s, k y ) is real. If density perturbation and the potential are in phase, the y -component of the force exerted on the planet is symmetric in the y -direction and there isno net one-sided torque when integrated over y .Now we consider the solution of equation (43). By Fourier transform in t -direction, weobtain σ ( s, k y ) = − ( s − s ) − isγ ( s − s ) + s γ k y ψ p ( s, k y ) , (52)and as equation (50) indicates, it is sufficient to consider s = (3 / p k y x mode.In the limit of small friction, γ → σ and ψ p are out of phase only at the resonance,which is located at s = 94 Ω k y x = s . (53)From equation (44), we see that this corresponds to the location of effective Lindblad reso-nance (Artymowicz 1993). There is a phase shift of π only in the vicinity of the resonance 18 –and torque is localized at the resonance location. This localization of the torque comes fromour assumption of zero radial pressure term in equation (40). We note that in this case, theresulting torque is independent of the amount of friction if integrated over the resonancewidth (Meyer-Vernet and Sicardy 1987). This can be seen by considering the imaginary partof the density fluctuation σ limits to, for sufficiently small γ ,lim γ → Im[ σ ( s, k y )] = lim γ → sγ ( s − s ) + s γ k y ψ p = πδ D ( s − s ) k y ψ p , (54)where δ D ( x ) is Dirac’s delta function. The last expression is independent of the amount offriction, and the torque exerted by the mode k y is determined by the strength of the source ψ p at the resonance. Using equations (44) and (53), we see that there is no resonance closeto the planet, or region | x | < x c , where x c = 23 c Ω p . (55)In this region, density perturbation and source term are in phase. In terms of real spacecoordinate ( x, y ), this corresponds to ψ p ( x, y ) ∝ σ ( x, y ).In the case of finite friction, the situation is different. From equation (52), we see thatresonance width is given by | s − s | ∼ γ . Therefore, when γ ∼ Ω p , the region in which densityperturbation and source term are out of phase can overlap the location of the planet. Notethat s is a function of k y and x [see equation (53)], so “the width of the resonance” refers tothe width in the x -direction, and this width varies with the mode in the y -direction ( k y ) weconsider. The finite width of the resonance causes the asymmetry in density perturbationin y > y < x ∼
0. Since gravitationalforce exerted by the planet on the disk is large in the vicinity of the planet, the significantamount of one-sided torque can be exerted on the planet.The amplitude of perturbation is suppressed when significant friction is exerted. Fromequation (52), we see that, in the vicinity of the planet | σ | ∼ p s + s γ k y ψ p ( s, k y ) . (56)The suppression is significant only when γ exceeds Ω p .One-sided torque is calculated by equation (24) (but there is no integral in the z -directionin two-dimensional problem considered in this section). Using equations (51) and (52), it ispossible to obtain T ∝ Z ∞ dk y k y Z ∞ ds sγ ( s − s ) + s γ ψ ( s, − k y ) , (57) 19 –where constant of proportionality does not depend on k y , γ , and x . The value of the torqueis determined by the competition between the suppression of the amount of the densityperturbation and the amplification of the amount of the torque as the resonance widthbecomes wider. Therefore, it is expected that the torque peaks at γ ∼ Ω p . This explainsthe peak of the torque for α ∼ ψ p ( s, k y ) does notdepend on s , ψ p ( s, k y ) = Ψ . (58)This corresponds to the case where forcing potential is constant in the x -direction. In thiscase, it is possible to perform s integral analytically, T ∝ Z ∞ dk y k y × p s − γ π + 2 tan − " s − γ γ p s − γ γ < s p γ − s log " γ − s + γ p γ − s γ − s − γ p γ − s γ > s (59)If there is only one k y mode in the potential, this gives the exact amount of torque. Otherwise,this gives the amount of torque exerted by each mode of k y . Limiting values of equation (59)are T → πk y s γ → T ∝ | γ − s | γ ∼ s T → γ → ∞ . (60)Therefore, the amount of torque is independent of the amount of dissipation in the limit of γ →
0, it increases as γ increases and peaks at γ ∼ s , and then it decreases to zero as γ becomes very large. In case of small γ , the integrand of (57) is localized at s ∼ s . When γ is large, however, it is necessary to consider contribution from all the region of s . In general,since the amplitude of forcing term ψ p ( s, k y ) increases as s →
0, it is expected that the valueof torque peaks at γ ∼ s . If the amplitude of forcing term cuts off at Ω ∼ c k y , the peakof the torque is expected at γ ∼ Ω p . 20 –
5. Summary and Discussion5.1. Summary of Our Results
In this paper, we have performed linear analyses of density structure as a result of grav-itational interaction between a planet and a viscous gas disk using shearing-sheet approxi-mation. We have obtained the density structure in the vicinity of the planet and discussedits effect on type I migration by calculating one-sided torque. We have used time-dependentmethods to obtain the density structure and a wide range of viscous parameter α is inves-tigated. We have calculated one-sided torque, which is the torque exerted on the planet bythe gas on one side of the planet’s orbit. One-sided torque converges to a well-defined valueif viscosity is small, but increases as we increase the value of viscosity for α ∼ > .
01, and thendecreases as α exceeds the order of unity. The values of one-sided torque we have obtainedfor large viscosity may be the lower limit of the actual value. We note that total torqueis not actually calculated in this paper since it is not possible within the shearing-sheetformalism. However, we also emphasize that the results of shearing-sheet calculations areuseful in calculating differential torque (Tanaka et al. 2002). We have seen that the torqueexerted on the planet in a viscous disk is primarily determined by the tilted spherical (orspheroidal) density structure in the vicinity of the planet, in contrast to the inviscid case,where main contribution of the torque comes from effective Lindblad resonances and thehorseshoe region. In an inviscid disk, although the value of density perturbation is large inthe vicinity of the resonance, it is symmetric in y -direction and therefore, torque exerted from y > y < y -direction,resulting in increasing torque with increasing viscosity until α ∼
1. By using a simplifiedanalytic calculations in Section 4, we have shown that it is possible to consider these effectsof viscosity as the increase of the width of Lindblad resonances. When the viscous coefficientis larger than α ∼ >
1, however, the torque becomes small with increasing viscous coefficientsince the density perturbation is damped due to large viscosity. The torque can be a factorof two larger than inviscid case. The impact of detailed structure around the planet ontype I migration has not yet been fully analyzed in previous numerical simulations probablybecause relatively large values of viscosity ( α ∼ > − − − ) are required to see this effect.Below, we show some discussions and future prospects. 21 – In this subsection, we discuss whether or not linear approximation well describes thereal density structure around the planet. In the linear approximation, the magnitude ofperturbed quantities such as δρ/ρ or δ v /c must be smaller than the order of unity. Inour calculation, all the non-dimensional perturbed quantity is proportional to the value of GM p /Hc , which is about 10 − − − for typical values of protoplanetary nebula. Therefore,when δρ/ρ becomes the order of ten in the normalized calculations given in this paper, linearapproximation becomes less accurate. Since the contribution of the torque in a viscousdisk mainly comes from such strongly perturbed regions, it may be necessary to performhigh resolution numerical calculation in order to fully obtain the structure of the densityfluctuation. This is one of the future prospects of this study. However, we expect thatthe basic physics and the order of magnitude of one-sided torque is similar to the linearcalculation. The total torque that is exerted on the planet is the differential torque, which is thedifference of the torque exerted by the inner disk ( x <
0) and outer disk ( x > r B = GM p /c . Assumingthat the differential torque T diff mainly comes from the effects of curvature, the order of themagnitude is T diff ∼ T one − side × δxr p , (61)where T one − side is one-sided torque obtained in this paper, and δx is the length scale wheremost contribution to the torque comes from. In the inviscid case, δx ∼ H and therefore T diff , inviscid ∼ T one − side , inviscid × Hr p . (62) 22 –In the viscous case presented in this paper, δx may be as small as Bondi radius r B . Therefore,the lower limit of the differential torque in viscous disk may be estimated as T diff , viscous ∼ T one − side , viscous × r B r p , (63)and the ratio of differential torque in viscous disk and that in inviscid disk may be given by T diff , viscous T diff , inviscid ∼ T one − side , viscous T one − side , inviscid × r B H . (64)Since r B H ∼ M p M c (cid:16) r p H (cid:17) = 8 × − (cid:18) − M p /M c (cid:19) (cid:18) . H/r p (cid:19) − , (65)and T one − side , viscous /T one − side , inviscid ∼ α ∼ In the shearing-sheet formalism, we cannot calculate corotation torque. For two-dimensionalmodes, this is because all the background quantities (and therefore vortensity) are assumedto be constant. For three-dimensional modes, we note that, in Tanaka et al. (2002), thereis a singularity at corotation modes (see equation (A5) of their paper). However, this singu-larity does not account for the torque in itself. Results of shearing-sheet calculations mustbe combined with the higher order solutions of modified local approximation in order to findthe torque exerted at corotation resonance, see equation (57) of Tanaka et al. (2002).If we assume that time derivatives of equations (6) and (7) are zero in order to obtaina stationary solution, we also find a similar singularity at x = 0 for three-dimensional waves( k z = 0), even we have assumed the constant background density in the z -direction. However,within the shearing-sheet formalism, it does not give a torque at the corotation. Also, thissingularity does not seem to play a major role in our time-dependent calculations, since thereis no singularity in our formulation. 23 –If we proceed to the modified local approximation, we expect that it is possible tocalculate the corotation torque in the linear regime. Recently, Paardekooper and Papaloizou(2009a) has shown that large viscosity will push the corotation torque into a linear regime.Therefore, not only from the point of view of differential torque but also from the view of thecorotation torque, extension of our time-dependent methods to modified local approximationis an interesting future work.In summary, in order to predict the precise value of the torque, it is necessary to use amodified local approximation or perform global calculations, which can take into account thecurvature effect. Non-linear analysis may be important since the main contribution comesfrom the region in the vicinity of the planet where density is strongly perturbed. However,we believe that the basic physics that exerts torque onto the planet can be captured by thislinear analysis and therefore, high-resolution study is necessary. In global calculations, wealso note that it is also possible to investigate the effects of gas accretion onto the centralstar on type I migration, which is always present in a viscous disk but is not captured in ourlocal calculations. Recently, Paardekooper and Papaloizou (2009b) has shown an interesting results re-garding the flow structure close to the planet. In this section, we observe streamline of theflow in the vicinity of the planet and compare our results with those previously studied. Wenote that since we have not calculated axisymmetric ( k y = 0) modes, the flow structure wehave obtained is incomplete. However, we can still find some qualitative difference betweeninviscid and viscous calculations.Figure 12 shows the streamlines plotted over the density fluctuations for calculationscorresponding to Figure 4. Results with α = 10 − and α = 10 − are shown. Calculations arerestricted to 2D modes of 3D calculation. Mass ratio between the planet and the central star q = M p /M c and disk aspect ratio h = H/r p are assumed in such a way that q/h = 0 . ∼ . H for α = 10 − , while the widthbecomes slightly narrower in calculations with α = 10 − . We have obtained about a factorof two smaller horseshoe width compared to Paardekooper and Papaloizou (2009b) in cal-culations with small viscosity (see Figure 3 of their paper.) We think that this is probablybecause we have weaker potential in the vicinity of the planet because of the effective soft-ening introduced by vertical averaging, see equation (31). In Paardekooper and Papaloizou 24 –(2009b), they also obtained narrower horseshoe width in calculations with larger softeningparameters. However, since we have neglected axisymmetric modes and non-linear effects,which may be potentially important in determining the streamline, this issue needs furtherinvestigation. Although we have neglected the axisymmetric modes in our calculations, a viscousoverstability in axisymmetric modes has been discussed, especially in the context of Saturn’srings. In this subsection, we discuss whether this overstability affects our results.Local linear analysis of viscous overstability using a simple hydrodynamic model isperformed by Schmidt et al. (2001). Stability criterion for isothermal case depends on thederivative of viscous coefficient ν with respect to surface density perturbation. We notethat the system of equations we have solved, equations (20)-(23), are actually stable againstviscous overstability, since our assumption of viscosity corresponds to β = − In this case, we do not expect viscous overstability to occur and therefore, wecan safely neglect the axisymmetric modes.However, disk may be prone to viscous overstability if different prescription of viscosityis taken into account. If the non-linear consequence of viscous overstability produces non-axisymmetric structure in the vicinity of the planet, this may exert additional torque ontothe planet. However, observation of the particle simulation performed by Schmidt et al.(2001) (their Figure 1 for example) indicates that non-axisymmetric modes are not stronglydriven by the viscous overstability and therefore this may not play an important role in plan-etary migration. Nonetheless, the effects of viscous overstability seems to be an interestingextension of the analysis presented here.
In this subsection, we briefly discuss that our analysis may indicate that the stochastictorque is important in a turbulent disk dominated by magneto-rotational instability (seee.g., Nelson and Papaloizou 2004). The effective turbulent viscosity originated in MRI maybecome as much as α ∼ . Compare our equations (6) and (7) and their equation (8). Note there is no self-gravity and bulkviscosity in our calculations and we assumed isothermal equation of state.
25 –the density structure around the planet is more important than the tidal wave launched atthe effective Lindblad resonances. Since the density structure around the planet is turbulentin the vicinity of the planet, stochastic force may be exerted on the planet embedded in sucha turbulent disk.In order for the α -viscosity prescription [see equation (27)] to be a good approximationof the turbulent flow, the length scale in consideration must be larger than the eddy size. Inthe present problem, the scale in consideration is smaller than disk thickness. In a turbulentflow driven by MRI, eddies of sizes on the order of several tenths of the disk scale lengthsexist. This can be understood if we notice that the wavelength of the most unstable modeis much smaller than the scale height of the disk with weak magnetic field. If weak poloidalmagnetic field, B z , is exerted on the unperturbed disk, there always exists a net magneticflux, h B z i = B z , (66)where h B z i is the horizontally averaged z -component of magnetic field. Therefore, therealways exists eddies with scales of the order of the most unstable mode of magneto-rotationalinstability with net flux B z , λ ∼ v A, Ω , (67)where λ is the eddy scale, v A, = B z / √ πρ is Alfv´en velocity constructed from the averagepoloidal field, and Ω is Keplerian angular frequency. In most cases, Alfv´en velocity is smallerthan sound speed by more than an order of a magnitude, v A, ∼ < − c. (68)This indicates λ ∼ < − H. (69)Thus, we can possibly use the effective turbulent viscosity in order to study the qualitativeeffects of the interaction between a turbulent disk and a planet.However, the effective values of α may be obtained by averaging turbulent stress overseveral scale heights, since it includes all sizes of eddy that is important in turbulent flowdriven by MRI. Thus, the use of α prescription may not be valid in discussing very smallscale structure, and full MHD calculation is necessary. More detailed, quantitative analysesrequire high-resolution numerical study of the magnetized turbulent disk.We have found that the effects of viscosity becomes apparent if α exceeds the value of ∼ . − .
1. Actual values of α in the disk is largely uncertain, but estimated values arearound this value. Numerical simulations by Sano et al. (2004) indicates that the values of α varies from 10 − − − depending on the setup. From the observational constraints of 26 –dwarf novae, which seem to be better studied than protoplanetary disks, it is indicated thatthe values of α can vary from 0 .
01 to 0 . α seemsto be an open issue for both theoretically and observationally.Authors thank Taku Takeuchi and Hidekazu Tanaka for useful discussions. We alsothank the anonymous referee for a number of suggestions that improved the paper. Thiswork was supported by the Grant-in-Aid for the Global COE Program “The Next Generationof Physics, Spun from Universality and Emergence” from the Ministry of Education, Culture,Sports, Science and Technology (MEXT) of Japan. The numerical calculations were carriedout on Altix3700 BX2 at YITP in Kyoto University. T. M. is supported by Grants-in-Aid forJSPS Fellows (19 · REFERENCES
Artymowicz, P., 1993 ApJ, 419, 155Balbus, S. A., & Hawley, J. F., 1991, ApJ, 376, 214Baruteau, C., & Masset, F., 2008, ApJ, 672, 1054Cannizzo, J. K., & Shafter, A. W., & Wheeler, J. C., 1988, ApJ, 333, 227D’Angelo, G., Henning, T., & Kley, W., 2002, A&A, 385, 647D’Angelo, Kley, W., & Henning, T., 2003, ApJ, 586, 540Fromang, S., Terquem, C., & Nelson, R.P., 2005, MNRAS, 363, 943Goldreich, P., & Lynden-Bell, D., 1965, MNRAS, 130, 125Goldreich, P., & Tremaine, S., 1979, ApJ, 233, 857Goodman, J., & Rafikov, R. R., 2001, ApJ, 552, 793Korycansky, D. G., & Pollack, J. B., 1993, Icarus, 102, 150Masset, F. S., 2001, ApJ, 558, 453Masset, F. S., 2002, A&A, 387, 605Meyer-Vernet, N., & Sicardy, B., 1987, Icarus, 69, 157 27 –Muto, T., Machida, M. N., & Inutsuka, S.-i., 2008, ApJ, 679, 813Narayan, R., Goldreich, P., & Goodman, J., 1987, MNRAS, 228, 1Nelson, R. P., & Papaloizou, J. C. B., 2004, MNRAS, 350, 849Paardekooper, S.-J., & Mellema, G., 2006, A&A, 459, L17Paardekooper, S.-J., & Papaloizou, J. C. B., 2008, A&A, 485, 877Paardekooper, S.-J., & Papaloizou, J. C. B., 2009a, MNRAS, 394, 2283Paardekooper, S.-J., & Papaloizou, J. C. B., 2009b, MNRAS, 394, 2297Papaloizou, J. C. B., & Lin, D. N. C., 1984, ApJ, 285, 818Press, W. H., Flannery, B. P., Teukolsky, S. A., & Vetterling, W. T., 1992, NumericalRecipes in Fortran 90: The Art of Scientific Computing , (Cambridge: CambridgeUniv. Press)Sano, T., Inutsuka, S.-i., Turner, N. J., & Stone, J. M., 2004, ApJ, 605, 321Schmidt, J., Salo, H., Spahn, F., & Petzschmann, O., 2001, Icarus, 153, 316Takeuchi, T., Miyama, S. M., & Lin, D. N. C., 1984, ApJ, 460, 832Tanaka, H., Takeuchi, T., & Ward, W. R., 2002, ApJ, 565, 1257Terquem, C. E. J. M. L. J., 2003, MNRAS, 341, 1157Ward, W. R., 1986, Icarus, 67, 164Ward, W. R., 1989, ApJ, 336, 526
This preprint was prepared with the AAS L A TEX macros v5.2.
28 – no r m a li z ed t o r que L z
2D mode of 3D calculation2D calculation
Fig. 1.— Comparison of the torque obtained in 2D calculation and the 2D mode of 3Dcalculation for various L z . In 2D calculations (dashed line), we use the gravitational po-tential given by equation (29), while in 2D mode of 3D calculation (solid line), we use thegravitational potential given by equation (31). The value of the normalized torque, T / Γ is plotted, where Γ is defined by equation (32). 29 – ψ p (x +y + ε ) Tanaka et al. zeroth order2D mode of unstratified disk
Fig. 2.— Comparison of the 2D mode of the potential used by Tanaka et al. (2002) given byequation (33) (solid line) and the k z = 0 mode of the potential given by (31) with L z = 2 . H (dashed line). 30 – -6 -5 -4 -3 -2 -1 no r m a li z ed t o r que α
2D mode torque
Fig. 3.— Variation of torque as a function of viscous coefficient α as a result of 3D calculationrestricted to 2D modes. 31 – α =10 -4 y / H α =10 -1 y / H Fig. 4.— Contour plot of the density fluctuation δρ/ρ around the planet (located at theorigin) for 2D mode of 3D calculation with L z = 2 . H with α = 10 − (left) and α = 10 − (right). Axisymmetric mode is not included. 32 – t o r que den s i t y x/H α =10 -4 α =10 -1 Fig. 5.— Torque density given by equation (26) for α = 10 − (solid line) and α = 10 − (dashed line). Values are normalized according to equation (36). 33 – α =10 -4 y / H α =10 -1 y / H Fig. 6.— Forward-back asymmetry of the torque distribution given by equation (37) for α = 10 − and α = 10 − . 34 – -6 -5 -4 -3 -2 -1 no r m a li z ed t o r que α all n z modesonly 2D modesfitting function Fig. 7.— Variation of the torque as a function of viscosity parameter α obtained for 3Dcalculations (plus). For comparison, we plot the results restricted to 2D modes by cross.Dashed line shows the fitting function given by equation (38). 35 – -2 0 2 4 6 8 10 12 0.01 0.1 1 t o r que den s i t y x/H α =10 -4 α =10 -2 α =10 -1 Fig. 8.— Torque density profile obtained by 3D calculation with various α . Solid line showsthe results of α = 10 − , thick dashed line α = 10 − , and thick dotted line α = 10 − . 36 – α =10 -4 α =10 -1 -0.1 -0.05 0 0.05 0.1y/H-0.1-0.05 0 0.05 0.1 z / H α =10 -4 α =10 -1 -0.1 -0.05 0 0.05 0.1y/H-0.1-0.05 0 0.05 0.1 z / H Fig. 9.— The contour plot of 3D mode density structure in yz -plane at x = 0 . H with α = 10 − (solid line) and α = 10 − (dashed line). The lines show the contours of δρ/ρ = 10. 37 – no r m a li z ed t o r que ε /H α =10 -4 α =10 -1 Fig. 10.— One-sided torque obtained by 3D calculation for different smoothing length.Horizontal axis shows the values of smoothing length ǫ and vertical axis shows the one-sided torque. Calculations with α = 10 − (solid line) and 10 − (dashed line) are shown. 38 – no r m a li z ed t o r que grid size (H) α =10 -4 α =10 -1 Fig. 11.— One-sided torque obtained by 3D calculation for different grid resolutions. Hori-zontal axis shows the values of grid size (the same for all three dimensions) and vertical axisshows the one-sided torque. Softening parameter ǫ = 10 − H is used. Calculations with α = 10 − (solid line) and 10 − (dashed line) are shown. 39 – α =10 -4
0 0.05 0.1 0.15 0.2 0.25x/H-3-2-1 0 1 2 3 y / H α =10 -4
0 0.05 0.1 0.15 0.2 0.25x/H-3-2-1 0 1 2 3 y / H -0.5 0 0.5 1 1.5 2 2.5 3 α =10 -1
0 0.05 0.1 0.15 0.2 0.25x/H-3-2-1 0 1 2 3 y / H α =10 -1
0 0.05 0.1 0.15 0.2 0.25x/H-3-2-1 0 1 2 3 y / H Fig. 12.— Density and streamline for calculations restricted to 2D mode in 3D calculationwith α = 10 − (left) and α = 10 − (right). This corresponds to the calculations presentedin Figure 4. Mass ratio between the central star q = M p /M c and disk aspect ratio h = H/r p are assumed in such a way that q/h = 0 . δρ/ρ divided by q/h3