Two-Dimensional Study of the Propagation of Planetary Wake and the Indication to Gap Opening in an Inviscid Protoplanetary Disk
aa r X i v : . [ a s t r o - ph . E P ] S e p ApJ, Accepted
Preprint typeset using L A TEX style emulateapj v. 11/10/09
TWO-DIMENSIONAL STUDY OF THE PROPAGATION OF PLANETARY WAKE AND THE INDICATION TOGAP OPENING IN AN INVISCID PROTOPLANETARY DISK
Takayuki Muto Department of Earth and Planetary Sciences, Tokyo Institute of Technology,2-12-1 Oh-okayama, Meguro-ku, Tokyo, 152-8551, Japan
Takeru K. Suzuki and Shu-ichiro Inutsuka
Department of Physics, Nagoya University,Furo-cho, Chikusa-ku, Nagoya, 464-8602, Japan
ApJ, Accepted
ABSTRACTWe analyze the physical processes of gap formation in an inviscid protoplanetary disk with anembedded protoplanet using two-dimensional local shearing-sheet model. Spiral density wave launchedby the planet shocks and the angular momentum carried by the wave is transferred to the backgroundflow. The exchange of the angular momentum can affect the mass flux in the vicinity of the planetto form an underdense region, or gap, around the planetary orbit. We first perform weakly non-linear analyses to show that the specific vorticity formed by shock dissipation of density wave canbe a source of mass flux in the vicinity of the planet, and that the gap can be opened even forlow-mass planets unless the migration of the planet is substantial. We then perform high resolutionnumerical simulations to check analytic consideration. By comparing the gap opening timescale andtype I migration timescale, we propose a criterion for the formation of underdense region around theplanetary orbit that is qualitatively different from previous studies. The minimum mass required forthe planet to form a dip is twice as small as previous studies if we incorporate the standard values oftype I migration timescale, but it can be much smaller if there is a location in the disk where type Imigration is halted.
Subject headings: planet and satellites: formation — protoplanetary disks — planet-disk interactions [email protected] JSPS Research Fellow
Muto, Suzuki and Inutsuka INTRODUCTION
Disk-planet interaction is one of the important topics in the planet formation theory. A low mass planet embeddedin a disk excites the density wave (Goldreich and Tremaine 1979), and the backreaction from the density wave causesthe planet to migrate in the disk (e.g., Ward 1986, Tanaka et al. 2002). The excitation of density wave at Lindbladresonances can be understood by linear analyses, although it has recently been pointed out that non-linear effects arealso important at corotation resonances (Paardekooper and Papaloizou 2009). For a high mass planet, the interactionbetween the planet and the disk becomes nonlinear, and the gap opens around the planetary orbit (Lin and Papaloizou1986ab, Ward and Hourigan 1989, Rafikov 2002b, Crida et al. 2006).Gap formation around the planet is important both theoretically and observationally. From the theoretical pointof view, gap formation determines the regime of planetary migration. If there is no gap opening, planetary migrationis in the regime called “type I”, where the interaction between the planet and the spiral density wave is important(Goldreich and Tremaine 1979, Ward 1986, 1997, Tanaka et al. 2002). Type I planetary migration timescale isconsidered to be faster than disk dispersal timescale, which poses a serious problem in the theory of planet formation.If the gap opens around the planet, the migration is in the regime called “type II”, where the planet migrates as thedisk accretion onto the central star occurs (Lin and Papaloizou 1986). The timescale of type II migration, which is ofthe order of viscous timescale, is generally longer than type I migration, and there may be a possibility for the planetsto survive in the disk. From observational point of view, the gap around the planet may be able to be observed by directimaging. Recent progress of disk observation by direct imaging has reached the stage that it is possible to compare thenumerical simulations and observation directly (Mayama et al. 2009). Moreover, the dynamical interaction betweenthe circumstellar dust or gas and the planet can be used to estimate the mass of a low-mass object embedded in thedisk (Kalas et al. 2008). If the gap in the disk can extend to the disk scale, the gap structure can be a very goodindicator of the existence of the planet.Conventionally, gap opening is understood as a balance between the torque exerted by the planet and the viscoustorque (Lin and Papaloizou 1986b). In addition to the viscous torque, Crida et al. (2006) has shown that “pressuretorque” also acts to balance the planetary torque. Gap opening processes are mainly investigated using one-dimensionalmodel. Rafikov (2002b) calculate the evolution of disk surface density in the vicinity of the planet using a classicalone-dimensional model (Lynden-Bell and Pringle 1974, Pringle 1981). He has shown that in the vicinity of the low-massplanet, the torque exerted by the planet is carried away by density wave , and the gap is not opened until the waveshocks to deposit the angular momentum to the mean flow. Considering the pressure effects as well as viscous effects,Crida et al. (2006) have derived the gap opening criterion generalized for arbitrary values of kinematic viscosity usingtwo-dimensional numerical simulation (equation (15) of their paper). Using the criterion of Crida et al. (2006), thegap-opening mass for an inviscid disk reads M p M ∗ ∼ > . × − (cid:18) H/r . (cid:19) , (1)where M p is the mass of the planet, M ∗ is the mass of the central star, H is the scale height of the disk, and r is theorbital semi-major axis of the planet. Therefore, planets more massive than Saturn can open up the gap in the disk.However, recently, Li et al. (2009) have performed high resolution numerical simulations and showed that a partialgap is formed in the vicinity of the planet even if the mass of the plant is smaller than the mass limit given by equation(1) in case the disk viscosity is very low. This motivates us to investigate the physical processes of gap opening in aninviscid disk in detail.In this paper, we show that low mass planets can potentially open up a partial gap. We investigate the processesby means of analytical models taking into account weak non-linearity. We also perform numerical simulations tolook at to what extent numerical calculations and analytic studies agree, and then finally we suggest the criterionfor the gap opening in an inviscid disk. For analytic studies, we study the propagation of the spiral density waveand subsequent shock formation. We then investigate the mass flux in the vicinity of the planet by means of second-order perturbation theory. We note that the second-order perturbation is necessary to study the mass flux since it isessentially a second-order quantity.The plan of this paper is as follows. In Section 2, we describe the basic equations. In Section 3, we investigate thegap opening processes using a second-order perturbation theory. We show that the shock dissipation of density waveand the subsequent formation of specific vorticity can lead to the mass flux in the vicinity of the planet, resultingin the gap formation. We then show the results of numerical simulations in Section 4. In Section 5, we discuss thecondition for gap opening in an inviscid disk. We also discuss that commonly used one-dimensional models of diskevolution may overlook gap opening processes in the vicinity of the planet. Section 6 is for summary. BASIC EQUATIONS
In this paper, we focus on the two-dimensional local shearing-sheet analysis for simplicity. We set up a local Cartesiancoordinate system corotating with a planet. We take the origin of the coordinate system at the planet’s location, andthe x - and y -axes are the radial and the azimuthal direction, respectively. We use isothermal ideal hydrodynamicequations ∂ Σ ∂t + ∇ · (Σ v ) = 0 (2) Note that Crida et al. (2006) explain this in terms of the balance between the pressure torque and the planetary torque. ap Opening in Inviscid Disk 3 ∂ v ∂t + v · ∇ v = − c Σ ∇ Σ − p e z × v + 3Ω p x − ∇ ψ p (3)where we have assumed that the gas is rotating at the Kepler velocity. Notations are as follows: Σ is surface density, v is velocity field, c is sound speed, Ω p is the angular velocity of the planet, and ψ p is the gravitational potential ofthe planet. For ψ p , we assume the form ψ p = GM p ( x + y + ǫ ) / , (4)where G , M p , and ǫ are the gravitational constant, the mass of the planet, and the softening parameter, respectively.Local shearing-sheet approximation is only an approximation of global model, and it may not be appropriate forinvestigating the global evolution of the disk structure. However, local shearing-sheet approximation and full globalmodel share many essential physics in common. Excitation and the propagation of density wave can be understoodusing local approximation. We also show later that one-dimensional disk evolution model constructed from globalmodel and local model are very similar. Local approximation also has the advantages that analyses are greatlysimplified and high-resolution calculations become possible.Here, we comment on the two-dimensional approximation. It is known that three-dimensional processes are im-portant in disk-planet interaction especially when considering the immediate vicinity of the planet. For example,numerical simulations by Paardekooper and Mellema (2006) clearly show that disk structure around the planet isthree-dimensional. However, for the study of density wave, which is the structure away from the location of the planet,the density perturbation by the planet is nearly two-dimensional. In this paper, we shall investigate the physical prop-erties of density wave excited by the planet in detail, and we discuss the gap formation processes that can be derivedfrom the study of density wave. Therefore, we expect that essential physics can be captured by local two-dimensionalmodel. ANALYTIC STUDY OF PLANETARY WAKE
In this section, we investigate the disk-planet interaction and subsequent gap opening processes using analyticmethods. We consider a weakly non-linear stage, where the mass of the planet is GM p Hc ∼ < , (5)which is approximately smaller than the Saturn mass in case of Minimum Mass Solar Nebula (Hayashi et al. 1985)with H/r ∼ . . We note that non-linearity of disk-planet interaction and gap opening is strongly related. Fromequation (5), the onset of the non-linearity is given by M p M ∗ ∼ > . × − (cid:18) H/r p . (cid:19) , (6)where H = c/ Ω p is the scale height of the disk. This non-linear criterion is analogous to equation (1), which is thegap opening criterion given by Crida et al. (2006)The overall picture of gap formation we suggest is as follows (see also Figure 1).1. Density wave excited by the planet will shock at some location away from the planet.2. The shock formation leads to the formation of specific vorticity.3. The change of specific vorticity results in net radial mass flux. This mass flux exists in the place closer to theplanet, even at the place where the change of specific vorticity is not significant.4. Gap opens in the vicinity of the planet.We investigate the shock formation process and the mass flux separately. Shock formation processes are investigated byGoodman and Rafikov (2001) and Rafikov (2002a) using Burgers equation model, and the mass flux is investigated byLubow (1990), using second-order perturbation theory. We show how these two theories can be combined to investigatethe gap opening processes.Physically, it is natural that the gap opens when spiral density wave damps, since the angular momentum flux carriedby the density wave should be transferred to the background flow. However, we shall show that the decay of the spiraldensity wave away from the planet can affect the mass flux in the immediate vicinity of the planet, in contrast to themodel suggested by Rafikov (2002b), in which there is no gap formation in the vicinity of the planet in inviscid cases. Equation (5) is equivalent to the condition where Hill’s radius becomes comparable to the disk scale height. Departure from linearcalculations can be observed at this mass range, see e.g., Miyoshi et al. (1999)
Muto, Suzuki and Inutsuka
Linear Analysis
In this section, as a preparation for the non-linear analyses performed in the subsequent sections, we briefly summarizethe results of linear study of the spiral density wave (e.g., Goldreich and Tremaine 1979, 1980). We denote backgroundstate by subscript “0”, and the perturbation by δ : Σ = Σ + δ Σ , (7) v = v + δ v = −
32 Ω p x e y + δ v . (8)The planet potential is regarded as a perturbation. The linearized equations are (cid:18) ∂∂t −
32 Ω p x ∂∂y (cid:19) δ ΣΣ + ∂∂x δv x + ∂∂y δv y = 0 (9) (cid:18) ∂∂t −
32 Ω p x ∂∂y (cid:19) δv x = − c ∂∂x δ ΣΣ + 2Ω p δv y − ∂∂x ψ p (10) (cid:18) ∂∂t −
32 Ω p x ∂∂y (cid:19) δv y = − c ∂∂y δ ΣΣ −
12 Ω p δv x − ∂∂y ψ p (11)We now assume the stationary state in the frame corotating with the planet, ∂/∂t = 0, and Fourier transform in the y -direction. The perturbed values are given by δf ( x, y ) = X n y ∈ Z δf ( x ) e ik y y , (12)where δf denotes δ Σ, δv x , or δv y , and k y is the wave number in the y -direction. Assuming the periodicity in the y -direction, k y = 2 πn y /L y , where n y is an integer and L y is the box size in the y -direction. The summation in theabove equation is taken over the relative integers denoted by n y . It is possible to derive a single second-order ordinarydifferential equation (Artymowicz 1993). d dx δv y +
94 Ω k y c x − k y − Ω c ! δv y = 32 Ω p k y c xψ p − c Ω p dψ p dx , (13)and other perturbed quantities are given by δ ΣΣ = 1 D (cid:20) Ω p ddx δv y + 32 Ω p k y xδv y − k y ψ p (cid:21) (14)and δv x = 1 D (cid:20) − c ik y ddx δv y + 34 Ω ik y xδv y − ik y Ω p ψ p (cid:21) , (15)where D = Ω c k y . (16)The boundary condition for Equation (13) is that wave should propagate away from the planet. At the locationaway from the effective Lindblad resonance given by94 Ω k y c x − k y − Ω c = 0 , (17)the approximate solution can be written analytically using WKB approximation . It takes the form, for | x | → ∞ , δv y ∼ C ( k y ) p | x | exp (cid:20) ± i
34 Ω p k y c x (cid:21) , (18)where upper sign is for x > x <
0, and C ( k y ) denotes the amplitude (but not necessarily real)of the wave with mode k y . We note that since we consider the linear perturbation excited by the planet’s gravitationalpotential at the origin of our coordinate system, the amplitude C ( k y ) is proportional to the planet mass M p . Fromhere on, using the symmetry of the shearing-sheet, we only consider x > δ Σ and δv x at | x | → ∞ , δ ΣΣ = C ( k y ) ik y Ω p √ xD (cid:20)
34 Ω p c − i k y (cid:21) exp (cid:20) i
34 Ω p k y c x (cid:21) (19) The exact solution of equation (13) is given by parabolic cylinder function, see Artymowicz (1993) for detail. ap Opening in Inviscid Disk 5 δv x = c δ ΣΣ , (20)where we have retained only the leading terms in x in equations (14) and (15). Note that, in the real space, solutionsare given by δf = X k y F ( k y ) x q exp (cid:20) ik y (cid:18) y + 34 x H (cid:19)(cid:21) , (21)where H is the scale height H = c/ Ω p , q = 1 / δ Σ and δv x , q = − / δv y , and F ( k y ) is a function of k y .Therefore, along the line y = − x H + y , (22)where y is constant, perturbed quantities take the same values except for the dependence x q . Therefore, we can writethe form of the WKB solution in real space, δ Σ( x, y ) = M p x / f ( y + (3 / x /H ) , (23) δv x ( x, y ) = M p x / g ( y + (3 / x /H ) , (24) δv y ( x, y ) = M p x − / h ( y + (3 / x /H ) , (25)where we write the dependence on the planet mass explicitly, and f , g , and h are the functions that determinesthe form of the perturbation. We note that from equation (20), f and g are proportional in the place where WKBapproximation is valid.The exact linear solution can be obtained by solving equations (9)-(11) numerically with non-reflecting boundaryconditions. The profiles of density and v x obtained in such a way are shown in Figure 2. In this figure, we calculatenon-axisymmetric modes ( k y = 0) and assumed GM p /Hc = 1. The amplitude of perturbation is proportional to theplanet mass. We note that the results obtained in equations (23)-(25) are based on the WKB approximation, and theyare applicable only in the region | x/H | ∼ > Shock Formation
Goodman and Rafikov (2001) performed a non-linear analysis of the propagation of the spiral density wave, andtheir local approach was later extended to the global model by Rafikov (2002a). Using several approximations basedon the linear theory, they derived the Burgers’ equation which describes the propagation of spiral density wave andconcluded that the spiral density wave eventually shocks as it propagates in the radial direction. They derived thatthe location of shock formation is proportional to M − / . In this section, we derive this relationship using a slightlydifferent consideration.Shock is formed when two characteristics cross. For two-dimensional supersonic steady flow, the gradient of charac-teristic curves is given by (Landau and Lifshitz 1959) (cid:18) dydx (cid:19) ± = v x v y ± c √ v − c v x − c . (26)In the background state of the shearing-sheet, this gives (cid:18) dydx (cid:19) ± = ∓ c (cid:18)
94 Ω x − c (cid:19) . (27)For x >
0, ( dy/dx ) + is the perturbation propagating away from the origin. At x ≫ (2 / H , the characteristic curvesare given by y ∼ − x H + y , (28)where y is a constant. It is to be noted that the outgoing characteristics coincides the curve of the same phase of theperturbation of the density wave (see equation (22)). This is because the density wave is essentially the sound wavepropagating on the disk.For the flow distorted by the perturbation of the planet, the characteristic curve is also distorted. Using the resultsof linear perturbation, the gradient of the characteristics is given by (cid:18) dydx (cid:19) ± ∼ c " ∓ c (cid:18)
94 Ω x − c (cid:19) + 32 Ω p x ( δv x ± cδv y (cid:0) (9 / x − c (cid:1) / ) , (29) Muto, Suzuki and Inutsukaupto the lowest order of perturbation. Since the amplitude of δv y decreases with x − / , δv y in the second term of theright hand side can be neglected for | x/H | ≫
1. Assuming that the perturbed characteristic curve is given by the form y ( x ) = − x H + y + δy ( x ) , (30)where δy ( x ) is given by ddx δy = 32 Ω p c xδv x ( x, y ( x )) . (31)Approximating y ( x ) ∼ − (3 / x /H + y in the argument of δv x , we have δy ( x ) ∝ M p x / g ( y ) , (32)where g ( y ) is the function that appears in the linear solution of δv x given by equation (24). From Figure 2, g ( y ) has azero point. In the vicinity of g ( y ) = 0, g ( y ) is negative for y > y and g ( y ) is positive for y < y . This indicates thatthe separation between two characteristic curves shrinks compared to the unperturbed case. Assuming that when δy exceeds a critical value, the shock forms, we obtain, for the location of the shock formation, x ∝ M − / , (33)which is the same condition as derived by Goodman and Rafikov (2001). Since the flow is supersonic only at x > (2 / H , condition (cid:12)(cid:12)(cid:12)(cid:12) x − H (cid:12)(cid:12)(cid:12)(cid:12) ∝ M − / (34)may be more appropriate. Second-Order Perturbation Analysis and Mass Flux
In this section, we consider second-order linear perturbation in order to derive the mass flux in the vicinity ofthe planet. Although angular momentum flux can be calculated using the results of linear perturbation only, it isnecessary to perform second-order analysis to derive the mass flux, since axisymmetric ( k y = 0) mode of the second-order perturbation contributes to the mass flux. Lubow (1990) performed the time-dependent analysis using Fourierand Laplace transformation. In this paper, we calculate second order perturbation in real space.The equations for second-order perturbation are given by (cid:18) ∂∂t −
32 Ω p ∂∂y (cid:19) δ Σ (2) Σ + ∂∂x δv (2) x + ∂∂y δv (2) y = − ∂∂x (cid:18) δ Σ (1) Σ δv (1) x (cid:19) − ∂∂y (cid:18) δ Σ (1) Σ δv (1) y (cid:19) (35) (cid:18) ∂∂t −
32 Ω p ∂∂y (cid:19) δv (2) x + c ∂∂x δ Σ (2) Σ − p δv (2) y = − δv (1) x ∂∂x δv (1) x − δv (1) y ∂∂y δv (1) x + c δ Σ (1) Σ ∂∂x δ Σ (1) Σ (36) (cid:18) ∂∂t −
32 Ω p ∂∂y (cid:19) δv (2) y + c ∂∂y δ Σ (2) Σ + 12 Ω p δv (2) x = − δv (1) x ∂∂x δv (1) y − δv (1) y ∂∂y δv (1) y + c δ Σ (1) Σ ∂∂x δ Σ (1) Σ . (37)The superscripts “(1)” and “(2)” denote the first- and second-order perturbation, respectively. We assume that thefirst-order results are already known.The mass flux is given by F M ( t, x ) = Σ δv (2) x + δ Σ (1) δv (1) x , (38)where bar denotes the integral over y . Assuming the periodic boundary condition in the y -direction, it is possible toderive the equation for F M which reads ∂ F M ∂t − c ∂ F M ∂x + Ω F M = S ( t, x ) , (39)ap Opening in Inviscid Disk 7where S is the source term consisting of two parts, S ( t, x ) = S v ( t, x ) + ∂∂t S t ( t, x ) , (40)where S v ( t, x ) = 2Ω p (cid:20) Ω p δ Σ (1) δv (1) x − Σ δv (1) x ∂ x δv (1) y (cid:21) (41)and S t ( t, x ) = ∂∂t (cid:16) δ Σ (1) δv (1) x (cid:17) − (cid:20) δv (1) x ∂ x δv (1) x + δv (1) y ∂ y δv (1) x − c Σ δ Σ (1) ∂ x δ Σ (1) (cid:21) (42)The term S v is related to the formation of specific vorticity. In two-dimensional ideal flow in a rotating frame, thespecific vorticity conserves along the streamline ddt ( ∇ × v ) z + 2Ω p Σ = 0 , (43)where d/dt ≡ ∂ t + v ·∇ is the Lagrangian derivative. In the linear perturbation analyses, this reduces to ∂ t − (3 / p ∂ y ,see also equations (35)-(37). The background value of the specific vorticity is Ω p / , and the linear perturbation is (cid:18) ∂∂t −
32 Ω p ∂∂y (cid:19) (cid:20) ∂∂y δv (1) x − ∂∂x δv (1) y + 12 Ω p δ Σ (1) Σ (cid:21) = 0 . (44)If there is no formation of specific vorticity, ∂∂y δv (1) x − ∂∂x δv (1) y + 12 Ω p δ Σ (1) Σ = 0 . (45)From equation (41), S v ( t, x ) becomes a total derivative with respect to y in this case, and therefore S v = 0. However,if there is a formation of vorticity by, for example, shock damping of the spiral density wave, this term can not beneglected. We also note that if stationary state is assumed a priori, and if there is no formation of specific vorticity,the source term S ( t, x ) is zero, leading to the zero mass flux (Lubow 1990, Muto and Inutsuka 2009a).However, if time-evolution effects are taken into account, we have non-zero mass flux. The solution for equation (39)is given by F M ( t, x ) = 1 c Z Z dt dx G ( t, t ; x, x ) S ( t , x ) , (46)where G ( t, t ; x, x ) is the Green’s function G ( t, t ; x, x ) = J (cid:18) H p c ( t − t ) − ( x − x ) (cid:19) | x − x | < c ( t − t )0 otherwise , (47)where J is the Bessel function of zeroth order. From linear perturbation, it is possible to predict that the mass fluxscales with M , since the source term is the second order of perturbation. Later in this paper, we compare this resultwith numerical calculations to understand how much mass flux is excited by the planet.We now consider the model in which the source term is given by S ( t, x ) = t < S [ δ D ( x − x s ) − δ D ( x + x s )] t > , (48)where S > x s > δ D ( x ) denotes the Dirac’s delta function. We later see that thesource term is positive in the region x >
0, and negative for x <
0. This form of the source term is the simplest casewhere we can obtain an analytic solution for the mass flux.Changing the integration variable from ( t , x ) to ( r, θ ) via x − x = r cos θ (49)1 H p c ( t − t ) − ( x − x ) = rH sin θ, (50)equation (46) can be rewritten to F M ( t, x ) = H c Z ct dr Z π dθ sin θJ (cid:16) rH sin θ (cid:17) S h t − rc , x − r cos θ i . (51) Muto, Suzuki and InutsukaIn order to see the mass flux only in the vicinity of the planet, we assume 0 < x < x s . Substituting equation (48), andintegrate over r , we obtain F M ( t, x ) = − HS c ( Z π cos − [( x − x s ) /ct ] dθ (cid:20) tan θJ (cid:18) x − x s H tan θ (cid:19)(cid:21) + Z cos − [( x + x s ) /ct ]0 dθ (cid:20) tan θJ (cid:18) x + x s H tan θ (cid:19)(cid:21) ) (52)We take the limit t → ∞ . Then, we can approximate cos − [( x ± x s ) /ct ] ∼ π/
2. Using the formula (Abramowitz andStegun 1970) Z π/ dθ tan θJ ( a tan θ ) = Z ∞ du u u J ( au ) = K ( a ) , (53)where a >
0, the integration over θ can be performed to obtain F M ( t, x ) ∼ HS c (cid:20) K (cid:18) x s − xH (cid:19) − K (cid:18) x s + xH (cid:19)(cid:21) , (54)where K is the modified Bessel function of the zeroth order. In the vicinity of the planet, x ≪ x s , the mass fluxchanges with F M ∝ x , which indicates that the gap opens up. We note that divergence at x = ± x s is the artefact ofour simplification where we have used delta function as a source term. Instability of a Disk with Surface Density Variation
We have seen that if there is a source of specific vorticity, mass flux appears in the vicinity of the planet, and itleads to the change of surface density to open a gap. We now briefly discuss how this gap opening process ends.An inviscid disk with a rapid surface density variation is prone to a linear instability, which is referred to as Rossbywave instability (Li et al. 2000, de Val-Borro et al. 2007). In Appendix B, we show a brief outline of the linear stabilityanalyses for a disk with a gap, and derive the necessary conditions for the instability. Gap induced by a planet in thedisk naturally excites the variation of specific vorticity, and gap edges are likely places where such instability occurs.Rossby wave instability may stop the gap-opening processes described in previous subsections. We do not expect thatthe instability leads to the complete closing of the gap, since the planet always try to repel the fluid element by theexcitation of the density wave. There may be at least “underdense region” around the planetary orbit, which may becalled a “gap”.
Summary of Analytic Study
We have discussed the non-linear shock formation of density wave in Section 3.2. The formation of shock leads tothe formation of specific vorticity in general, and the perturbation of specific vorticity can lead to the radial mass fluxas discussed in Section 3.3. The gap formation may end by the onset of the linear instability of the disk with a gap.In order to investigate this qualitative picture of gap formation, we perform two-dimensional numerical calculations inthe subsequent section.We note that in this picture of gap formation, the decay of angular momentum flux carried by the spiral densitywave is not directly connected to the mass flux. For example, when we model the source of mass flux by δ -functionin equation (48), we implicitly assume that the spiral density wave does not dissipate at | x | < x s , but we still end upwith non-zero mass flux in this region, see equation (54). One may consider that this result violates the conservationof angular momentum. However, it is actually satisfied, when time evolution effects are taken into account. We shalldiscuss further on the angular momentum conservation in Section 5.1. NUMERICAL STUDY OF PLANETARY WAKE
Numerical Setup
We solve Euler equations (2) and (3) with isothermal equation of state using second-order Godunov scheme (e.g.,Colella and Woodward 1984). In order to investigate the density wave propagation while resolving all the lengthscales (Hill radius, Bondi radius, and radial wavelength of the density wave), we need to use rather large box size withrelatively high resolution. We choose the box size ( L x , L y ) = (16 H, H ) with mesh number ( N x , N y ) = (512 , x = ∆ y = H/
32. We have also used the box with different L y to investigate the boxsize effect.From equation (21), the radial wavelength of the wave with mode k y may be approximated by λ ∼ H π k y H Hx . (55)The most important modes are k y H ∼ O (1). For mode with k y H = 10 at x/H ∼
4, the wavelength of the densitywave is given by λ ∼ H/
4, and therefore, the radial wavelength is resolved by eight meshes if we use ∆ x = H/ | x | < H in the analyses of the numerical results.If we normalize the length by H = c/ Ω p and the time by Ω − , the only dimensionless parameter in this calculationis the planet mass, or Bondi radius r B /H = GM p /Hc and the softening parameter ǫ/H . We use the softening lengthused in previous local shearing-sheet calculations ǫ = r H /
4, where r H is the Hill’s radius r H H = (cid:18) M p M ∗ (cid:19) / r p H = (cid:18) GM p Hc (cid:19) / (56)(Miyoshi et al. 1999). We have performed calculations with five different planetary mass M p shown in Table 1. Notethat Bondi radius, Hill radius, and the softening parameter are resolved (at least marginally) for the smallest massmodel. We have increased the planet mass linearly from t Ω p = 0 to 12. The result does not depend significantly onthis timescale if the timescale is longer than this.In the x -direction, we use non-reflecting boundary used by FARGO (Baruteau 2008), modified for the shearing-sheet.For the y -direction, we adopt the periodic boundary condition. Shearing-sheet calculations performed by Miyoshi etal. (1999) or Tanigawa and Watanabe (2002) used a different boundary condition in the y -direction, which is thecombination of Keplerian inflow and supersonic outflow. This boundary condition may be appropriate to study theflow structure only in the vicinity of the planet, but for the study of gap formation, this boundary condition isinappropriate because this forces the unperturbed gas flowing into the computational domain. We also note thatperiodic boundary condition in the y -direction is useful for the purpose of comparison with the linear analyses, whichassume the periodicity in the y -direction. We have checked that our code reproduces the results of Miyoshi et al.(1999) well when the same parameter and the boundary conditions are used. Results
Disk Structure and Evolution
Figures 3 and 4 show the δ Σ and δv x at t Ω p = 200 obtained by numerical simulations. For the numerical calculationswith GM p /Hc = 0 .
4, we can see that a gap is already formed. Figure 5 shows the evolution of the azimuthally-averaged density profile for GM p /Hc = 0 . .
4. It is possible to see that even for low-mass calculations, densitygap is being gradually formed. In subsequent sections, we provide the physical interpretations of disk structure andevolution using analytical theories provided in Section 3.
Spiral Shock Wave
We first investigate the shock formation. Figure 6 shows the density profiles at various x at t Ω p = 200 obtained bynumerical calculations. It is possible to see the shock-like structure for the calculation with GM p /Hc = 0 .
4, whilefor the run with GM p /Hc = 0 .
1, the structure is not very obvious.It is possible to see whether dissipation acts or not by looking at the perturbation of specific vorticity. In theabsence of dissipation, the specific vorticity conserves along the streamline, and since the background specific vorticityis constant in the calculation box, we expect that it is also constant in the presence of the planet. Note that anypotential force, whether it is time-dependent or not, does not produce specific vorticity. Therefore, specific vorticityarises only if dissipative mechanisms come into play. In our numerical calculation, the dissipation is implemented inthe shock-capturing scheme. If the formation of specific vorticity is driven by the formation of shock, relation givenby equation (34) should be observed at least for the weak shock cases.Figure 7 shows the evolution of the perturbation of azimuthally averaged specific vorticity for calculations with GM p /Hc = 0 . .
4. It is possible to see that the specific vorticity is formed around x/H ∼ .
5. We later seethat this specific vorticity becomes a dominant source of mass flux, which leads to the formation of a gap.For the numerical calculation with GM p /Hc = 0 .
1, we also see the formation of specific vorticity around x ∼ GM p /Hc = 0 .
4, similarthing happens, but the vortices formed by shock dissipation are stronger. In subsequent sections, we show that vorticesjust in the vicinity of the planet location does not produce a significant amount of mass flux.Figure 8 shows the azimuthally-integrated specific vorticity perturbation as a function of ( x − (2 / H ) M / . If shockdissipation occurs as equation (34), the peak of the specific vorticity perturbation comes at the same location of thehorizontal axis. It is possible to observe that the for calculations with GM p /Hc < .
2, equation (34) is marginallysatisfied. For calculations with larger planet mass, it is not the case. This is because the shock formation occursimmediately after the excitation of the wave at x ∼ (2 / H , and therefore in the unit of ( x − (2 / H ) M / , the shockoccurs relatively further away. Mass Flux
We now look at the radial mass flux Σ v x excited by the planet. In figure 9, we show the azimuthally averaged massflux obtained by numerical calculation. It is possible to see that there is a spike of the mass flux just in the vicinity of0 Muto, Suzuki and Inutsukathe planet, and then the mass flux depends linearly with the distance from the planet within | x/H | ∼ <
2. The spike inthe vicinity of the planet is due to the gas falling onto the planet. This causes the spike of the surface density profilein the vicinity of the planet as shown in figure 5. We expect that this feature depends on the treatment of the planet,and such mass flux should be investigated carefully when one wants to look at the processes such as gas accretion ontothe planet. However, as we will see later, the feature in this spike region does not affect the mass flux outside thisregion. In the discussion below, we focus on the structure outside this spike region.We first compare the results obtained by numerical simulation and analytic calculation. Figure 10 compares themass flux derived by numerical calculation and by equation (46). Since it is not possible to predict the amount ofspecific vorticity only from linear perturbation theory, we have used the results of numerical calculation in obtainingthe source term S ( t, x ). It is possible to observe that the second-order perturbation theory is valid for calculationswith low-mass planet, while for calculations with GM p /Hc = 0 .
4, the perturbation theory fails to explain the amountof mass flux. We find that the second-order perturbation theory can be used for GM p /Hc ∼ < . GM p /Hc = 0 .
1. We see that the term with ∂ t S ( t, x ) is not significant after the planet mass isfully introduced at t Ω p = 12. The term S v ( t, x ) depends on time strongly in the vicinity of the planet, while nearlytime-independent contribution is observed from the region | x/H | > .
5. The source arising in the vicinity of the planetcomes from the vortex making a horseshoe orbit, while the source at | x/H | > . | x/H | >
1. Figure 12 compares the mass flux obtained by doing so and the mass fluxobtained by using the full source term. We have used equation (46) with the source terms obtained by numericalcalculation to derive the mass flux. It is clear that the source in the vicinity of the planet does not contribute to themass flux very much. Therefore, we conclude that the mass flux arises from the specific vorticity generated by theshock dissipation of density wave.We now discuss the dependence of mass flux on the box size in the y -direction, L y . Figure 13 shows the azimuthallyintegrated mass flux R dy Σ v x with GM p /Hc = 0 . L y . These mass fluxes are directlycalculated from the numerical simulations. The lines show a perfect match for different box sizes. This can beexplained if we notice that the angular momentum flux carried by the density wave does not depend on the box size.As the spiral density wave damps, the angular momentum carried by the wake is deposited to the background disk todrive the mass flux. The total amount of the angular momentum deposited over the y -direction does not depend onthe box size, since the amount of angular momentum flux carried by the wave does not depend on the box size.We now discuss the timescale for the gap-opening in the vicinity of the planet. As stated at the end of Section 3.3,the mass flux in the vicinity of the planet is expected to be proportional to x , and if second-order analysis is valid,we expect that the mass flux is proportional to the square of the planet mass. We therefore fit the mass flux within | x/H | < L y Z dy Σ v x = K Σ c xH (cid:18) GM p Hc (cid:19) (cid:18) L y /H (cid:19) − , (57)where K is a constant that is derived by fitting the results of the numerical simulations. We found K = 4 . × − can fit the numerical results upto GM p /Hc ∼ < . | x/H | ≤
2. Figure 9 compares the mass fluxobtained by numerical simulation and the fitting function (57).Since the gap opening timescale may be determined by τ − = (1 /L y ) ∂ x R dy Σ v x (1 /L y ) R dy Σ , (58)and the denominator is approximated by Σ , the gap-opening timescale is proportional to L y . Using the result ofequation (57), we have τ − = 4 . × − Ω p (cid:18) GM p Hc (cid:19) (cid:18) L y /H (cid:19) − . (59)Note that this timescale applies to the initial phase of gap opening, and it does not state about how deep the gap willbe. The box size in the y -direction may be interpreted as a circumference of the disk. The above expression can berewritten as τ − = 1 . × − Ω p (cid:18) GM p Hc (cid:19) (cid:18) H/r p . (cid:19) (cid:18) πr p L y (cid:19) . (60) Instability of Gap Profile
In Section 3.4, we have seen that the disk with a density bump can be prone to linear instability. Figure 14 showsthe evolution of azimuthally averaged density perturbation and mass flux in the calculation with GM p /Hc = 0 . − H forthe parameter GM p /Hc ∼ . DISCUSSION: GAP FORMATION IN A PROTOPLANETARY DISK
In this section, using the results of numerical calculations, we discuss the validity of conventional one-dimensionalmodel for gap formation and the condition of gap opening.
One-Dimensional Model
The gap formation processes are modeled using one-dimensional disk evolution model. Thommes et al. (2008)uses a model in which the torque exerted at the Lindblad resonances are directly deposited to the disk to changethe semi-major axis of the fluid particle. Rafikov (2002b) uses a one-dimensional model based on Lynden-Bell andPringle (1974) type surface density evolution equation to model gap formation processes. We investigate whether suchtreatment can reproduce the evolution of the density profile obtained in two-dimensional calculations.Although we have performed shearing-sheet analyses, we first note that global calculation and local shearing-sheetcalculation share very similar property if we consider the azimuthally-averaged one-dimensional model. This can beinvestigated by comparing the one-dimensional model derived using global model and local model (see Appendix A).Exact one-dimensional evolution model is given by equations (A12) and (A15) in local approximation. Whetherthe evolution model (A17) can be used depends on whether the approximation made in deriving this equation isappropriate. This can be seen whether mass flux, angular momentum, and the torque exerted by the planet arerelated as equation (A16). In Figure 15, we compare each term of equation (A16). We observe that the relationshipof this equation is not satisfied, indicating the time evolution of Σ δv y occurs. We have checked that the exact timeevolution equation (A15) is satisfied.We therefore argue that time evolution of density and rotation profile is given by equation (A12) and (A15), and inconstructing one-dimensional model, it is necessary to specify the model of the mass flux, as well as the torque andangular momentum flux. It is especially important in considering the gap formation in the vicinity of the planet orbit.For example, Rafikov (2002b) predicts nothing happens in the vicinity of the planet in an inviscid model unless thespiral density wave shocks. In two-dimensional simulation, in contrast, we also observe gap opening even in low-massplanet calculations. As the gap opens, rotation velocity of the gas is also modified in order to balance the pressuregradient, as Crida et al. (2006) pointed out. Gap Opening Condition in an Inviscid Disk
We have seen that low-mass planets can potentially open a gap, or at least a ring with low surface density, in aninviscid disk. We now discuss the minimum mass of the planet that can open the gap. Low mass planets are proneto type I migration and therefore, if the planet migrates before it opens a gap, gap formation is impossible. Sincethe width of the gap is the order of scale height H , we compare the timescale of the planet migration over length H and the gap formation timescale. The planet migration timescale in an isothermal disk is estimated by Tanaka et al.(2001) as τ − ∼ p M p M ∗ Σ r M ∗ (cid:16) r p H (cid:17) . (61)Note that the power of H/r p is 3 because we consider the migration over the radial length H . For gap formationtimescale, we use the results of numerical calculations given in Section 4.2. The gap opening timescale τ − is givenby equation (59). If τ − < τ − , we expect the gap will open in the vicinity of the planet. Comparing τ mig and τ gap ,we obtain M p M ∗ ∼ > . × − (cid:18) H/r p . (cid:19) (cid:18) L y /H (cid:19) Σ2 × gcm − (cid:16) r p (cid:17) (cid:18) M ∗ M ⊙ (cid:19) − . (62)We have used typical values of protoplanetary disk at 1AU to estimate the number. We note that the gap-opening massscales with the box size L y . If we assume L y = 2 πr p and H/r p = 0 . L y /H = 125 .
66, which gives M p /M ∗ ∼ > × − ,2 Muto, Suzuki and Inutsukawhich is below the criterion given in previous studies, e.g., equation (1). Equation (62) rewritten by using 2 πr p reads(see also equation(60)) M p M ∗ ∼ > . × − (cid:18) H/r p . (cid:19) (cid:18) L y πr p (cid:19) Σ2 × gcm − (cid:16) r p (cid:17) (cid:18) M ∗ M ⊙ (cid:19) − . (63)It is interesting to compare equation (62) with the previous criterion (1). Equation (62) is derived by comparingthe gap opening timescale and type I migration timescale, and this is qualitatively different from the way the previouscriterion is derived. For example, the gap opening criterion derived in this study, equation (62), depends on the diskmass, since type I migration timescale scales with the disk mass. We also note that the type I migration timescalemay be much longer if the corotation torque and non-barotropic effects are considered (Paardekooper et al. 2010ab).In this case, our results indicate that the gap may be opened for lower mass planets, although more thorough analysesof non-barotropic effects on spiral density wave and mass flux are necessary.If the disk is in turbulent state, the effective viscosity can also act to fill the gap. If we use standard α -prescriptionfor turbulent viscosity ν = αcH , the timescale of viscous diffusion over length scale ∼ H is τ − ∼ α Ω p . (64)If τ vis is shorter than τ gap , we expect that the gap is filled. This leads to another condition, α ∼ < . × − (cid:18) H/r p . (cid:19) − (cid:18) M p /M ∗ × − (cid:19) (cid:18) L y /H (cid:19) − , (65)or if we use 2 πr p instead of L y , α ∼ < . × − (cid:18) H/r p . (cid:19) − (cid:18) M p /M ∗ × − (cid:19) (cid:18) πr p L y (cid:19) . (66)The disk must be quiet in order for a low-mass planets to open up the gap, but we note that the condition is verysensitive to the disk aspect ratio. We note that the previous results for gap-opening criterion derived by Crida etal. (2006) is somewhat different from equation (65). Since we derive equation (65) using a crude order-of-magnitudeestimate for viscous diffusion, it is necessary to perform more systematic parameter study in order to investigate howmuch viscosity is necessary to halt gap-opening by a planet.The gap-opening criteria, equations (62) and (65), are more exactly the conditions for the “formation of under-denseregion around the planetary orbit”. We compare the timescale for the decrease of surface density, which is derivedfrom the mass flux induced by the planet, and the type I migration timescale or viscous diffusion timescale. Therefore,these conditions concern the initial phase of the gap-opening, and the timescale argument does not concern how deepthe gap will be as a result of long-term evolution. The condition by Crida et al. (2006) is, on the other hand, thecondition for the formation of “deep gap”, or more quantitatively the formation of the gap where 90% of the originalamount of gas is depleted. Therefore, the conditions are not necessarily in contradiction with each other. For example,in the case of the numerical calculation with GM p /Hc = 0 .
6, 3 H/ r H ∼ .
3, which is slightly above the gap-openingcriterion by Crida et al. (2006). In figure 14, we observe a gap where approximately 50 −
60% of the original gasmass is depleted, which is just below the depth of the gap which is predicted by Crida et al. (2006) This may suggestthat the prediction by Crida et al. (2006) and our calculation is in qualitative agreement, although the setting of ourcalculation is different from theirs. Our condition, however, still has a new aspect since it states that the gap-openingprocess may be different in case of inviscid case. Non-linear evolution of density wave and the feedback onto the diskcan be important when we consider an inviscid disk. SUMMARY AND FUTURE PROSPECTS
In this paper, we have investigated the disk-planet interaction and the evolution of surface density profile of the diskusing analytic methods and high-resolution numerical calculations within the framework of shearing-sheet approxima-tion. We have shown some analytic framework to understand the non-linearity of disk-planet interaction. Formationof specific vorticity by shock dissipation of density wave can be a source of disk mass flux in the vicinity of the planet,which results in the formation of gap around the embedded planet. We estimate the conditions of the formation of theunderdense region, or gap formation, around the planetary orbit, and it is indicated that a planet with 20-30 Earthmass at 1AU will open a gap in a standard Minimum Mass Solar Nebula with
H/r ∼ .
05. We note that our conditionapplies to the initial phase of the gap opening, and the condition depends on the type I migration timescale, whichcan be much longer than that given by Tanaka et al. (2002) formula. If type I migration is halted, much lower massplanets can open a gap. In an inviscid case, which we have explored in this paper, the width of the gap is of the orderof the disk scale height. We note that we have not discussed the final gap depth in detail in this paper, although it isindicated that the gap opening process ends by the onset of hydrodynamic instability when approximately half of themass is depleted. One may call such a shallow underdense region a “dip” or “partial gap”, rather than gap. We havealso discussed that classical one-dimensional model of disk evolution fails to describe gap opening process, and morerefined one-dimensional is necessary to explain gap opening.We have focused on the two-dimensional local shearing-sheet calculations. We have discussed that the final width ofthe gap may be larger than we have obtained in this calculation. The gap depth and width we have obtained may beap Opening in Inviscid Disk 13interpreted as minimum values. Global calculations with high resolution are necessary to construct more quantitativemodel for gap formation. We also note that the discussion using the specific vorticity we have mentioned a number oftimes in this paper may be appropriate only in the two-dimensional calculations. Although we expect that the two-dimensional model is adequate for the investigation of spiral density wave qualitatively, three-dimensional calculationsmight result in the quantitatively different condition for gap opening.Authors thank the referee for useful comments that improved the paper. This work was supported by the Grant-in-Aid for the Global COE Program “The Next Generation of Physics, Spun from Universality and Emergence” fromthe Ministry of Education, Culture, Sports, Science and Technology (MEXT) of Japan. The numerical calculationswere in part carried out on Altix3700 BX2 at YITP in Kyoto University. Data analyses were in part carried out onthe data analysis server at Center for Computational Astrophysics, CfCA, of National Astronomical Observatory ofJapan. The page charge of this paper is supported by CfCA. T. M. is supported by Grants-in-Aid (22 · APPENDIX
DERIVATION OF ONE-DIMENSIONAL MODEL
In this section, we compare one-dimensional dynamical evolution model derived from global and local model, andshow that they share very similar properties.
Global Model
We use a cylindrical coordinate system ( r, φ ). In an inertial frame, a set of fluid equations in two-dimensional globalmodel is given by ∂ Σ ∂t + 1 r ∂∂r ( r Σ v r ) + 1 r ∂∂φ (Σ v φ ) = 0 , (A1) ∂v r ∂t + v r ∂v r ∂r + v φ r ∂v r ∂φ − v φ r = − ∂p∂r − d Ψ ∗ dr − ∂ψ p ∂r , (A2) ∂v φ ∂t + v r ∂v φ ∂r + v φ r ∂v r ∂φ − v r v φ r = − r ∂p∂φ − r ∂ψ p ∂φ , (A3)where Ψ ∗ is the gravitational potential of the central star which is assumed to be axisymmetric. The rotation profileof the background disk is given by r Ω ( r ) = d Ψ ∗ dr . (A4)We now derive a one-dimensional model for disk evolution by taking the azimuthal average of equation of continuity(A1) and conservation of angular momentum, which is essentially equation (A3). We decompose the azimuthal velocityinto background part and perturbation, v φ = r Ω( r ) + δv φ (A5)but we do not use linear approximation.Azimuthally averaged (denoted by bar) equation of continuity is ∂ Σ ∂t + 1 r ∂∂r (cid:0) r Σ v r (cid:1) = 0 , (A6)and the azimuthally averaged angular momentum conservation is ∂∂t (cid:0) r Σ v φ (cid:1) + 1 r ∂∂r (cid:0) r ΩΣ v r + r Σ v r δv φ (cid:1) = T ( r, t ) , (A7)where T ( r, t ) = − Σ ∂ φ ψ p (A8)is the torque exerted by the planet.Decomposing v φ as equation (A5), equation (A7) can be further rewritten, with the aid of equation (A6), ∂∂t (cid:0) r Σ δv φ (cid:1) + 1 r ∂∂r (cid:0) r Σ v r δv φ (cid:1) = − ( r Ω) ′ T ( r, t ) , (A9)where ′ denotes the derivative with respect to r .Upto equation (A9), our treatment is exact. In the standard one-dimensional model, approximation r Ω ≫ δv r isused to neglect the time-derivative of equation (A9) to obtain the relationship between the mass flux and angular4 Muto, Suzuki and Inutsukamomentum flux (Balbus and Papaloizou 1999). This approximation means that one neglects the evolution of therotation profile. We then obtain the relation between mass flux, angular momentum flux, and torque, (cid:0) r Ω (cid:1) ′ Σ v r = − r ∂∂r (Σ v r δv φ ) + T ( r, t ) , (A10)and the evolution of surface density is obtained from equation (A6) ∂ Σ ∂t − r ∂∂r (cid:20) r Ω) ′ ∂∂r (cid:0) r Σ v r δv φ (cid:1)(cid:21) = − r ∂∂r (cid:20) rT ( r, t )( r Ω) ′ (cid:21) . (A11)If α -prescription is used for the angular momentum flux and the planet is absent, this is the model by Lynden-Belland Pringle (1974). Rafikov (2002b) used this equation to study gap formation. Local Model
It is possible to derive equations analogous to global model in averaging over the y -direction in the shearing-sheetapproximation. Taking the y -average of equation of continuity, we obtain ∂ Σ ∂t + ∂∂x Σ v x = 0 (A12)and from the equation of motion in the y -direction in a conservation form, ∂∂t Σ v y −
23 Ω p x ∂∂x Σ v x + ∂∂x Σ v x δv y = − p Σ v x + T loc ( t, x ) , (A13)where T loc = − Σ ∂ y ψ p . (A14)Decomposing the first term in equation (A13) into the background and perturbation, we obtain ∂∂t Σ δv y + ∂∂x Σ v x δv y = −
12 Ω p Σ v x + T loc ( x ) . (A15)This corresponds to equation (A9) in the global model.If we approximate that | (3 / p x | ≫ δv y , mass flux and angular momentum flux are related by12 Ω p Σ v x = − ∂∂x Σ v x δv y + T loc ( x ) . (A16)we obtain ∂ Σ ∂t + ∂∂x (cid:20) − p ∂∂x Σ v x δv y (cid:21) = − p T loc ( x ) , (A17)which corresponds to equation (A11) in global model.We also note that the equation for mechanical energy loss is also analogous in global and local model. LINEAR STABILITY ANALYSIS OF A DISK WITH A GAP
In this section, we show the outline of the linear stability analysis of a disk when the surface density structure is notuniform.We consider a disk without a planet for simplicity. The equations we consider are equations (2) and (3) without ψ p .We assume that the background disk is axisymmetric with density profile Σ ( x ). We denote the background valueswith subscript “0”. In the background state, the pressure gradient must be balanced by Coriolis force and therefore, v y, ( x ) = −
32 Ω p x + c p d Σ dx ≡ U ( x ) , (B1)and v x, = 0 . (B2)We now consider linear perturbation. Perturbed values are denoted by δ and consider the solution proportional toexp[ − iωt + ik y y ]. Linear perturbation is then − i ˜ ω δ ΣΣ + 1Σ d Σ dx δv x + ddx δv x + ik y δv y = 0 , (B3) − i ˜ ωδv x + c ddx δ ΣΣ − p δv y = 0 , (B4)ap Opening in Inviscid Disk 15 − i ˜ ωδv y + c ik y δ ΣΣ + (cid:18) p + dUdx (cid:19) δv x = 0 , (B5)where ˜ ω ( x ) ≡ ω − k y U ( x ). From these equations, we can derive a single second-order ordinary differential equationfor Φ ≡ Σ δv x , ddx (cid:20) (˜ ω − c k y ) d Φ dx (cid:21) + " k y p ˜ ω (cid:18) κ Σ (˜ ω − c k y ) (cid:19) ′ + ˜ ω − c k y − κ c Σ (˜ ω − c k y ) Φ = 0 , (B6)where κ ( x ) = 2Ω p (2Ω p + U ′ ( x )) (B7)and ′ denotes the derivative with respect to x . Equation (B6) with an appropriate boundary condition constructs aneigenvalue problem for ω . Since full analyses of equation (B6) is not the scope of this paper, we briefly outline thequalitative results about the instability of gap profile that can be derived from equation (B6).For an axisymmetric mode ( k y = 0), equation (B6) becomes ddx (cid:20) d Φ dx (cid:21) + ω − κ c Σ Φ = 0 . (B8)We multiply Φ ∗ to this equation and integrate over x . Assuming that the perturbation vanishes at | x | → ∞ andintegrating by part, we obtain Z dx (cid:12)(cid:12)(cid:12)(cid:12) d Φ dx (cid:12)(cid:12)(cid:12)(cid:12) + Z dx κ c Σ | Φ | = Z dx ω c Σ | Φ | . (B9)We first note that from this equation, the eigenvalue ω must be real. For instability, ω must be negative andtherefore, κ must be negative at some point in x . This necessary condition for instability is the well-known Rayleighcriterion.We can further proceed by using the independent variable g ( x ) = 1 √ Σ Φ . (B10)Equation (B8) now becomes d gdx + " ω c − Ω c + 1 √ Σ d dx p Σ ! g = 0 (B11)Defining E = ω c (B12)and V ( x ) = Ω c + 1 √ Σ d dx p Σ , (B13)equation (B11) can be looked as a Scr¨odinger equation with energy E and potential V ( x ), d gdx + [ E − V ( x )] g ( x ) = 0 . (B14)If there is a “bound state” with energy eigenvalue E <
0, the system is unstable. Necessary condition for instabilityis therefore that there exists a point where V ( x ) <
0, or1Σ / d dx Σ / < − Ω c . (B15)Therefore, if there exists a sharp pressure bump, the system is unstable against axisymmetric perturbation.For a non-axisymmetric mode ( k y = 0), the analysis becomes more complex. However, we can make a progress whenwe assume that there is a mode confined close to corotation, ˜ ω = 0. Assuming ˜ ω ∼ ω − c k y ∼ − c k y , equation(B6) becomes ddx (cid:20) d Φ dx (cid:21) − " c k y + κ c Σ − ω k y p (cid:18) κ Σ (cid:19) ′ Φ = 0 (B16)6 Muto, Suzuki and InutsukaMultiplying Φ ∗ to this equation and integrate over x , we obtain Z dx " (cid:12)(cid:12)(cid:12)(cid:12) d Φ dx (cid:12)(cid:12)(cid:12)(cid:12) + c k y + κ c Σ | Φ | = Z dx ω k y p (cid:18) κ Σ (cid:19) ′ | Φ | (B17)Taking the imaginary part of this equation, we haveIm( ω ) Z dx | ˜ ω | k y p (cid:18) κ Σ (cid:19) ′ | Φ | = 0 . (B18)Therefore, should imaginary part of ω exist, ( κ / Σ ) ′ must change the sign somewhere. This is Rossby wave instabilitypreviously investigated by a number of authors (Lovelace and Hohlfeld 1978, Lovelace et al. 1999, Li et al. 2000). Ouranalysis presented here for non-axisymmetric disk is the shearing-sheet version of that given by Lovelace and Hohfeldt(1978), although we have used a different variables from their analysis. We note that κ / Σ is proportional to thebackground vortensity, which is (2Ω p + U ′ ) / Σ in the shearing-sheet approximation. REFERENCESAbramowitz, M. & Stegun, I. A. 1970, Handbook OfMathematical Functions (New York: Dover)Artymowicz, P. 1993, ApJ, 419, 155Balbus, S. A., & Papaloizou, J. C. B. 1999, ApJ, 521, 650Baruteau, C. 2008, PhD. ThesisChandrasekhar, S. 1981, Hydrodynamic and HydromagneticStability, (New York: Dover)Colella, P., & Woodward, P. R. 1984, Journal Comp. Phys., 54,174Crida, A., Morbidelli, A., & Masset, F. 2006, Icarus, 181, 587Goldreich, P., & Tremaine, S. 1979, ApJ, 233, 857Goldreich, P., & Tremaine, S. 1980, ApJ, 241, 425Gooodman, J., & Rafikov, R. R. 2001, ApJ, 552, 793Hayashi, C., Nakazawa, K., & Nakagawa, Y. 1985, in Protostarsand Planets II, ed. Black & Matthews (Tucson: Univ. ArizonaPress)Kalas, P. et al. 2008, Science, 322, 1345Landau, L. D., & Lifshitz, E., M. 1959, Fluid Mechanics. Courseof Theoretical Physics, (Oxford: Pergamon)Lovelace, R. V. E., & Hohlfeld, R. G. 1978, ApJ, 221, 51Lovelace, R. V. E. Li, H., Colgate, S. A., & Nelson, A. F. 1999,ApJ, 513, 805Li, H. Finn, J. M., Lovelace, R. V. E., & Colgate, A. 2000, ApJ,533, 1023Li, H. Lubow, S., & Lin, D. N. C. 2009, ApJ, 690, L52Lin, D. N. C., & Papaloizou, J. C. P. 1986a, ApJ, 307, 395Lin, D. N. C., & Papaloizou, J. C. P. 1986b, ApJ, 309, 846Lubow, S. H. 1990, ApJ, 362, 395 Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603Mayama, S. et al. Science, in pressMiyoshi, K., Takeuchi, T., Tanaka, H., & Ida, S. 1999, ApJ, 516,451Muto, T., & Inutsuka, S.-i. 2009a, ApJ, 695, 1132Muto, T., & Inutsuka, S.-i. 2009b, ApJ, 701, 18Paardekooper, S.-J., & Mellema, G. 2006, A&A, 459, L17Paardekooper, S.-J., & Papaloizou, J. C. P. 2009, MNRAS, 394,2283Paardekooper, S.-J., Baruteau, C., Crida, A., & Kley, W. 2010a,MNRAS, 401, 1950Paardekooper, S.-J., Baruteau, C., & Kley, W. 2010b, MNRAS,Accepted: arXiv:1007.4964Pringle, J. E. 1981, ARA&A, 19, 137Rafikov, R. R. 2002a, ApJ, 569, 997Rafikov, R. R. 2002b, ApJ, 572, 566Tanaka, H., Takeuchi, T., & Ward, W.R. 2002, ApJ, 565, 1257Tanigawa, T., & Watanabe, S.-i. 2002, ApJ, 580, 506Thommes, E. W., Matsumura, S., & Rasio, F. A. 2008, Science,321, 814de Val-Borro, M., Artymowicz, P., D’Angelo, G., & Peplinski, A.2007, A&A, 471, 1043Ward, W. R. 1986, Icarus, 67, 164Ward, W. R. 1997, Icarus, 126, 261Ward, W. R. & Hourigan, K. 1989, ApJ, 347, 490 ap Opening in Inviscid Disk 17 radialazimuthal1. density wave shocks2. specific vorticity formation3. mass fluxPlanetdensity wave
Fig. 1.—
Overall picture of the physical processes involved in gap opening. The density wave launched by the planet eventually shocksas it propagates in the radial direction. Dissipation of the density wave then leads to the formation of the specific vorticity. The specificvorticity leads to the mass flux around the planet, which results in the gap formation.
Fig. 2.—
Linear density perturbation δ Σ / Σ (left) and δv x (right). We use GM p /Hc = 1, but the perturbation scales with the mass ofthe planet. ap Opening in Inviscid Disk 19 Fig. 3.—
Simulation results of δ Σ / Σ (left) and δv x (right) at t Ω p = 200. We use GM p /Hc = 0 . Fig. 4.—
Simulation results of δ Σ / Σ and δv x at t Ω p = 200. We use GM p /Hc = 0 . ap Opening in Inviscid Disk 21 a z i m u t ha ll y a v e r aged den s i t y x/H t Ω p =50t Ω p =100t Ω p =150 0.75 0.8 0.85 0.9 0.95 1 1.05 1.1-4 -3 -2 -1 0 1 2 3 4 a z i m u t ha ll y a v e r aged den s i t y x/H t Ω p =50t Ω p =100t Ω p =150 Fig. 5.—
The evolution of azimuthally averaged density profile for GM p /Hc = 0 . . den s i t y y/H 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 -6 -4 -2 0 2 den s i t y y/H Fig. 6.—
Density profile at t Ω p = 200. Profiles at x/H = 1, 1 .
5, 2, 2 .
5, and 3 are shown from right to left. Left panel is for GM p /Hc = 0 . GM p /Hc = 0 . ap Opening in Inviscid Disk 23 Fig. 7.—
Time evolution of azimuthally averaged specific vorticity perturbation for GM p /Hc = 0 . . -20 0 20 40 60 80 100 120 140 160 180 200 0 0.2 0.4 0.6 0.8 1 s pe c i f i c v o r t i c i t y pe r t u r ba t i on ( a r b i t r a r y un i t ) (x/H-2/3)*(GM p /Hc ) GM p /Hc =0.05GM p /Hc =0.1GM p /Hc =0.2GM p /Hc =0.4GM p /Hc =0.6 Fig. 8.—
Perturbation of specific vorticity (in arbitrary unit) as a function of ( x − (2 H/ × M / . ap Opening in Inviscid Disk 25 -0.02-0.015-0.01-0.005 0 0.005 0.01 0.015 0.02-4 -3 -2 -1 0 1 2 3 4 no r m a li z ed m a ss f l u x x/H numerical result, GM p /Hc =0.1fitting function -0.02-0.015-0.01-0.005 0 0.005 0.01 0.015 0.02-4 -3 -2 -1 0 1 2 3 4 no r m a li z ed m a ss f l u x x/H numerical result, GM p /Hc =0.4fitting function Fig. 9.—
Comparison between the mass flux obtained from numerical calculation (dots) and the fitting function (57) with K = 4 . × − .We plot the value ( GM p /Hc ) − (1 /L y ) R dy Σ v x . Left panel shows the results with GM p /Hc = 0 . GM p /Hc = 0 .
4. The horizontal axis shows x/H . -0.01-0.005 0 0.005 0.01-4 -3 -2 -1 0 1 2 3 4 a v e r aged m a ss f l u x x/H numerical calculation2nd-order perturbation -0.1-0.05 0 0.05 0.1-4 -3 -2 -1 0 1 2 3 4 a v e r aged m a ss f l u x x/H numerical calculation2nd-order perturbation Fig. 10.—
Comparison of the mass flux obtained by numerical calculation and second-order perturbation theory. Left panel shows theresult with GM p /Hc = 0 .
1, and the right panel shows the results with GM p /Hc = 0 . ap Opening in Inviscid Disk 27 Fig. 11.—
The evolution of the source terms of mass flux obtained by the calculations with GM p /Hc = 0 .
1. Left panel shows S v ( t, x )given by equation (41). Right panel shows ∂ t S t ( t, x ) where S t ( t, x ) is given by equation (42). -0.004-0.002 0 0.002 0.004 -4 -3 -2 -1 0 1 2 3 4 a v e r aged m a ss f l u x x/H source term only at |x/H|<1full source term Fig. 12.—
Mass flux obtained by second order perturbation theory. Solid line shows the case where the source term is restricted to | x/H | <
1, while the dashed line shows the result with full source term. ap Opening in Inviscid Disk 29 -0.03-0.02-0.01 0 0.01 0.02 0.03 -4 -3 -2 -1 0 1 2 3 4 i n t eg r a t ed m a ss f l u x x/H L y =16HL y =32HL y =64HL y =128H Fig. 13.—
Integrated mass flux R dy Σ v x for box sizes with L y /H = 16 , , , GM p /Hc = 0 . t Ω p = 100 are used. -0.6-0.5-0.4-0.3-0.2-0.1 0 0.1 0.2 0.3 -4 -2 0 2 4 a z i m t ha ll y a v e r aged den s i t y x/H t Ω p =150t Ω p =250t Ω p =350t Ω p =400 -0.2-0.15-0.1-0.05 0 0.05 0.1 0.15 0.2 -4 -2 0 2 4 a z u m u t ha ll y a v e r aged m a ss f l u x x/H t Ω p =150t Ω p =250t Ω p =350t Ω p =400 Fig. 14.—
Evolution of density profile (left) and the mass flux (right) for the calculation with GM p /Hc = 0 . ap Opening in Inviscid Disk 31 -0.001-0.0005 0 0.0005 0.001 -4 -2 0 2 4x/H (mass flux)/2derivative of angular momentum fluxtorquetotal Fig. 15.—
Comparison of each term in equation (A16). Data with GM p /Hc = 0 . t Ω p = 100 is used. The line “(mass flux)/2” isthe left hand side of equation (A16), “derivative of angular momentum flux” is the first term of the right hand side (sign inverted), and“torque” is the second term of the right hand side. The line “total” is the value when all the terms are taken to the left hand side, andshould be zero if equation (A16) is strictly satisfied. TABLE 1Mass parameters used in numerical calculations
Model Number GM p /Hc Planet Mass a r H /H . M ⊕ . M ⊕ . M ⊕ M ⊕ . M ⊕ a Assuming 1AU and
H/r p = 0 ..