Formation of a disc gap induced by a planet: Effect of the deviation from Keplerian disc rotation

0.97 MB·English

by  Kazuhiro D. Kanagawa

Mon.Not.R.Astron.Soc.000,1–15(2014) Printed23January2015 (MNLATEXstylefilev2.2)
Formation of a disc gap induced by a planet: Effect of the
deviation from Keplerian disc rotation
5
1 K.D. Kanagawa1⋆, H. Tanaka1, T. Muto2, T. Tanigawa3 and T. Takeuchi4
0 1Institute of Low Temperature Science, Hokkaido University,Sapporo 060-0819, Japan
2 2Divisionof Liberal Arts, Kogakuin University,1-24-2, Nishi-Shinjuku, Shinjuku-ku, Tokyo, 163-8677, Japan
n 3School of Medicine, University of Occupational and Environmental Health, Yahatanishi-ku, Kitakyushu, Fukuoka 807-8555, Japan
a 4Department of Earth and Planetary Sciences, TokyoInstitute of Technology, Meguro-ku, Tokyo 152-8551, Japan
J

2
2 23January2015

]
P
ABSTRACT
E
Thegapformationinducedbyagiantplanetisimportantintheevolutionoftheplanet
h. and the protoplanetary disc. We examine the gap formation by a planet with a new
p formulation of one-dimensional viscous discs which takes into account the deviation
- from Keplerian disc rotation due to the steep gradient of the surface density. This
o
formulation enables us to naturally include the Rayleigh stable condition for the disc
r
t rotation. It is found that the derivation from Keplerian disc rotation promotes the
s
radialangularmomentumtransferandmakesthegapshallowerthanintheKeplerian
a
[ case. For deep gaps, this shallowing effect becomes significant due to the Rayleigh
condition.Inourmodel,wealsotakeintoaccountthepropagationofthedensitywaves
1
excited by the planet, which widens the range of the angular momentum deposition
v
to the disc. The effect of the wave propagation makes the gap wider and shallower
2
than the case with instantaneous wave damping. With these shallowing effects, our
2
one-dimensional gap model is consistent with the recent hydrodynamic simulations.
4
5
Key words: accretion, accretion discs, protoplanetary discs, planets and satellites:
0
formation
.
1
0
5
1
1 INTRODUCTION Lubowet al. 1999; Kley 1999), which is a possible mecha-
:
v nism for forming an inner hole in the disc (Zhuet al. 2011;
i A planet in a protoplanetary disc gravitationally interacts
X Dodson-Robinson & Salyk 2011).
with the disc and exerts a torque on it. The torque ex-
r erted by the planet dispels the surrounding gas and forms Becauseoftheirimportance,discgapsinducedbyplan-
a
a disc gap along the orbit of the planet (Lin & Papaloizou ets have been studied by many authors, using simple one-
1979;Goldreich & Tremaine1980).However,agasflowinto dimensional disc models (e.g., Takeuchi et al. 1996; Ward
the gap is also caused by viscous diffusion and hence the 1997; Crida et al. 2006; Lubow& D’Angelo 2006) and nu-
gap depth is determined by the balance between the plan- merical hydrodynamic simulations (Artymowicz & Lubow
etary torque and the viscous diffusion. Accordingly, only a 1994;Kley1999;Varni`ere et al.2004;Duffell & MacFadyen
large planet can createa deepgap (Lin & Papaloizou 1993; 2013; Funget al. 2014). One-dimensional disc models pre-
Takeuchiet al. 1996; Ward 1997; Rafikov2002;Crida et al. dict an exponential dependence of the gap depth. That is,
2006). the minimum surface density at the gap bottom is pro-
The gap formation strongly influences the evolution portional to exp[ A(Mp/M )2], where Mp and M are the
of both the planet and the protoplanetary disc in various masses of the pla−net and th∗e central star, and A ∗is a non-
ways. For example, a deep gap prevents disc gas from ac- dimensionalparameter(seealsoeq.[38]).Ontheotherhand,
creting onto the planet and slows down the planet growth recent high-resolution hydrodynamic simulations done by
(D’Angelo et al. 2002; Bate et al. 2003; Tanigawa & Ikoma Duffell& MacFadyen(2013,hereafterDM13)showthatthe
2007), and also changes the planetary migration from gap is much shallower for a massive planet than the pre-
the type I to the slower type II (Lin & Papaloizou 1986; diction of one-dimensional models. According to their re-
Ward 1997). Furthermore, a sufficiently deep gap in- sults, the minimum surface density at the gap is propor-
hibits gas flow across the gap (Artymowicz & Lubow 1996; tionalto(Mp/M )−2.Varni`ereet al.(2004)andFunget al.
∗
(2014) obtained similar results from their hydrodynamic
simulations.Itsorigin hasnotyetbeenclarified bytheone-
⋆ E-mail:[email protected] dimensional disc model. Funget al. (2014) also estimated
2 K. D. Kanagawa
thegapdepthwitha“zero-dimensional”analyticmodel,by 2 MODEL AND BASIC EQUATIONS
simply assuming that the planetary gravitational torque is
We examine an axisymmetric gap in the disc surface den-
produced only at the gap bottom. Their simple model suc-
sity around a planet by using the one-dimensional model
ceedsin explainingthedependenceoftheminimumsurface
density of (Mp/M )−2. However, the zero-dimensional of viscous accretion discs. Although the Keplerian angular
∝ ∗ velocity is assumed in most previous studies, we take into
model doesnot givetheradial profileof thesurface density
accountadeviationfromKepleriandiscrotationinourone-
(orthewidthofthegap).Itisnotwellunderstoodwhatkind
dimensional model. The deviation cannot be neglected for
ofprofileacceptstheirassumption on theplanetarytorque.
a deep gap, as will shown below. We also assume non-self-
Further development of the one-dimensional gap model is
gravitating and geometrical thin discs. For simplicity, the
required in order to clarify both the gap depth and width.
planet is assumed to be in a circular orbit. We also adopt
Such a model enables us to connect the gaps observed in
simple models for density wave excitation and damping to
protoplanetary discs with theembedded planets.
describe thegap formation.
Oneoftheproblemsof theone-dimensional discmodel
istheassumptionoftheKeplerianrotationalspeed.Thedisc
rotation deviates from the Keplerian speed due to a radial
2.1 Angular velocity of a protoplanetary disc with
pressure gradient (Adachiet al. 1976). When a planet cre-
a gap
atesadeepgap,thesteepsurfacedensitygradientincreases
thedeviationofthediscrotationsignificantly,whichaffects The angular velocity, Ω, of a gaseous disc around a central
theangularmomentumtransferatthegap(seeSections2.1 star with mass M is determined by the balance of radial
∗
and2.2).Furthermore,alargedeviationofthediscrotation forces:
cdaisncsa(lsCohavniodlraatseekthhearR1a9y3l9e)ig.hAsvtiaoblalteiocnonodfitthioenRfaoyrlerigohtactoinng- Ω2R− GRM2∗ − Σ1 ∂∂PR2D =0, (1)
ditionpromotestheangularmomentumtransferandmakes
whereRistheradialdistancefromthecentralstar,Σisthe
the surface density gradient shallower so that the Rayleigh
condition is only marginally satisfied (Tanigawa & Ikoma surface density of the disc and P2D denotes the vertically
averaged pressure. On the left-hand side of equation (1),
2007; Yang& Menou 2010). To examine such feedback on
thefirsttermrepresentsthecentrifugalforceonaunitmass
the surface density gradient, we should naturally include
ofthedisc.Thesecondandthirdtermsarethegravitational
the deviation from the Keplerian disc rotation in the one-
forcebythecentralstarandtheforceoftheradialpressure
dimensional disc model.
Another simplification is in the wave propagation at gradient,respectively.ForP2D,weadoptthesimpleequation
the disc–planet interaction. The density waves excited by of state P2D =c2Σ, where c is the isothermal sound speed.
Using this equation of state, equation (1) can be rewritten
planets radially propagate in the disc and the angular mo-
as
menta of the waves are deposited on the disc by damping.
Thisangularmomentumdepositionisthedirectcauseofthe Ω2 =Ω2 (1 2η), (2)
K −
gapformation.Mostpreviousstudiessimplyassumeinstan-
with
taneousdamping of thedensitywaves after theirexcitation
(e.g.,Ward1997;Crida et al.2006).Ifthewavepropagation h2 ∂lnΣ ∂lnc2
η= + , (3)
is taken into account, the angular momentum is deposited −2R ∂R ∂R
(cid:18) (cid:19)
in a wider region of the disc, which increases the width of
thegap(Takeuchiet al.1996;Rafikov2002).Inawidegap, where ΩK = GM /R3 is the Keplerian angular velocity,
∗
the disc–planet interaction would be weak because the disc and h=c/ΩKp
gas around the planet decreases over a wide region. Hence, For a disc with no gap, the order of magnitude of
wecannotneglecttheeffectofwavepropagationonthegap thenon-dimensionalparameterηisO(h2/R2)(Adachi et al.
formation. 1976), because the term in parentheses in equation (3) is
In thepresent paper, we re-examine thegap formation comparable to 1/R.Ontheotherhand,if aplanet opens
∼
by a planet with the one-dimensional disc model, taking adeepgapwithawidthof h,thesteepgradientofthesur-
∼
into account the deviation from Keplerian rotation and the face density increases η to O(h/R). Hence, it also enhances
effect of wave propagation. To include the deviation from the deviation of the disc rotation from the Kepler rotation
Keplerian disc rotation, we modify the basic equations for inthedeepgap.Weneglecttheterm∂c/∂Rinequation(3)
one-dimensionalaccretiondiscs,detailedinthenextsection. because the temperature gradient would be small. Neglect-
Theeffectofthewavepropagationisincludedusingasimple ing the smaller terms of O(h/R), we approximately obtain
model. In Section 3, we obtain estimates of gap depths for ∂Ω/∂R as
ttwheozseimrop-dleimcaesnessi.oOnanlemesotdimelapterofporosaedwibdyeFgaupngcoertreaslp.o(2n0d1s4t)o. ∂Ω = 3ΩK 1 h2∂2lnΣ . (4)
∂R − 2R − 3 ∂R2
InSections4and5,wepresentnumericalsolutionsofthegap (cid:20) (cid:21)
without and with wave propagation, respectively. We find Note that the second term in the parentheses is of order
thatthegapbecomesshallow duetotheeffectsofthedevi- unity since d2lnΣ/dx2 1/h2 in a deep gap. Therefore, it
∼
ation from Keplerianrotation, theviolation oftheRayleigh is found that ∂Ω/∂R is significantly altered from the Kep-
condition and thewavepropagation. Withtheseshallowing lerianvalueduetothesteepgradientofthesurfacedensity,
effects, our results are consistent with the recent hydrody- thoughthedeviationofΩissmall( h/R).Asshownlater,
∼
namic simulations. In Section 6, we summarize and discuss this deviation promotes radial viscous transfer of angular
our results. momentum and makes a gap shallower.
Formation of a disc gap induced by a planet 3
2.2 Basic equations describing a disc gap around wouldbevalidiftheaccretionrateontotheplanetissmaller
a planet than theradial disc accretion rate, FM.
Under these assumptions, equation (5) shows that FM
Theequationsforconservationofmassandangularmomen-
is constant. Equation (6) yields
tum are given by
∞
∂∂Σt + 2π1R∂∂FRM =SM, (5) FJ =FJ(∞)−ZR ΛddR′, (8)
and where FJ(∞) is the angular momentum flux without the
planet. From equations (7) and (8),we obtain
∂∂t(Σj)+ 2π1R∂∂FRJ =jSM+ 2π1RΛd, (6) jFM 2πR3ΣνdΩ =FJ( ) ∞ΛddR′. (9)
− dR ∞ −
where FM and FJ are theradial fluxesof mass and angular ZR
momentum, and j(= R2Ω) is the specific angular momen- Equation (9) with a constant mass accretion rate describes
tum. In equation (5), the source term SM represents the a steady disc gap around a planet for a given Λd. Since
massaccretionrateontoaunitsurfaceareaofthedisc.The dΩ/dR is given by equation (4), equation (9) is thesecond-
accretion of disc gas onto theplanet can be included in SM order differential-integral equation. Note that equation (9)
as a negative term. In equation (6), Λd(R) represents the is derived from equation (6) and indicates the angular mo-
deposition rate of the angular momentum from the planet mentum conservation.
on thering region with radius R. By differentiating equation (9), we obtain a rather fa-
Todescribethedeposition rateΛd,weconsiderthean- miliar expression for the mass flux:
gular momentum transfer from the planet to the disc. This
transfer process can be divided into two steps. First the FM = ddRj −1 ddR 2πRp3νΣddRΩ +Λd . (10)
planet excites a density wave by the gravitational interac- (cid:18) (cid:19) (cid:20) (cid:18) (cid:19) (cid:21)
tion with the disc (e.g., Goldreich & Tremaine 1980). Sec- Note that this expression is valid only for the steady state.
ond, the density waves are gradually damped due to the In a time-dependent case, equation (10) should include the
disc viscosity or a nonlinear effect (Takeuchiet al. 1996; term 2πRΣ(∂j/∂t)intheparentheses.Asaboundarycon-
−
Goodman & Rafikov 2001). As a result of the wave damp- dition, the disc surface density should approach its unper-
ing,theangularmomentaofthewavesaredepositedonthe turbedvalues at both sides of thegap far from theplanet.
disc. If instantaneous wave damping is assumed, the depo- Here,wealsoconsidertheunperturbedsurfacedensity.
sition rateΛd isdeterminedonlybythewaveexcitation.In In the unperturbed state, Ω can be replaced by ΩK, by ne-
Section 2.4.1, we will describe the deposition rate for the glecting the smaller term O(h2/R2) (see equations [2] and
casewithinstantaneouswavedamping.InSection 2.4.2,we [3]).Furthermore,settingΛd =0inequation(9),weobtain
will give a simple model of Λd for the case of gradual wave theunperturbedsurface density Σ0 as
damping.
TheradialangularmomentumfluxFJ isgivenby(e.g., 3πR2νΩKΣ0(R)=−R2ΩKFM+FJ(∞). (11)
Lynden-Bell& Pringle 1974)
ThusΣ0 is given by
∂Ω
FJ=jFM−2πR3νΣ∂R. (7) Σ0 =−3FπMν (cid:18)1− RF2JΩ(K∞F)M(cid:19). (12)
The first term is the advection transport by the disc radial
Thisagrees with thewell-known solution for steady viscous
mass flow, FM, and the second term represents the viscous accretion discs (e.g., Lynden-Bell& Pringle 1974).
transport.Forthekineticviscosity,weadopttheαprescrip-
tion,i.e.,ν =αch(Shakura& Sunyaev1973).NotethatFJ
does not include the angular momentum transport by the
2.3 Rayleigh condition
density waves in ourformulation.
Equations (5)–(7) describe the time evolution of the For a deep gap around a large planet, the derivative of the
three variables Σ, FM and FJ with the given mass source angular velocity deviates significantly from the Keplerian
term SM, the angular momentum deposition rate from a velocity, as shown in Section 2.1. A sufficiently large de-
planet Λd and the disc angular velocity Ω. Note that Ω de- viation in Ω violates the so-called Rayleigh stable condi-
pendson ∂Σ/∂R, as in equations (2) and (3). tion of dj/dR >0 (see, Chandrasekhar 1961). Such a steep
Nextweconsiderthediscgapinasteadystate(∂/∂t= gap is dynamically unstable, which would cause a strong
0). The time scale for the formation of a steady gap is angular momentum transfer, lessening the steepness of the
approximately equal to the diffusion time within the gap gap.Thiswouldmaketheunstableregionmarginally stable
width, tdiff = h2/ν. For a nominal value of α ( 10−3), (i.e.,dj/dR=0).
the diffusion time is roughly given by 103 Kepler∼ian peri- Using equation (4), we give dj/dR as
ods,whichisshorterthanthegrowthtimeofplanets(105–7 dj 1 d2lnΣ
yprla)n(eKtaorkyubdoisc&s I(d1a06–2700y0r,)2(0H02a)iscohretthael.lif2e00t1im).eHoefncperotthoe- dR = 2RpΩKp(cid:18)1+h2p dR2 (cid:19), (13)
assumption of a steady gap would be valid. In addition, wherethesuffixp indicatesthevalueat R=Rp;thissuffix
we assume SM = 0 for simplicity. Although gas accretion is also used for other quantities. Hence, using the second-
onto the planet occurs for Mp & 10M (Mizuno 1980; derivativeof thesurface density,the marginally stable con-
⊕
Kanagawa & Fujimoto 2013), the assumption of SM = 0 dition dj/dR = 0 can be rewritten as (Tanigawa & Ikoma
4 K. D. Kanagawa
2007). where C = (25/34)[2K0(2/3) +K1(2/3)]2/π 0.798 and
≃
K denote the modified Bessel functions. The sign of equa-
d2lnΣ i
h2p dR2 =−1. (14) tthioenc(lo1s7e)visicpinoistiytivoef tfhoreRpla>neRt,pRor nRegpat.ivehfpo,roRn t6heRopt.hIenr
| − |
Actually,aroundasufficientlylargeplanet,equation(9) hand, the WKB formula is overestimated. Thus, we model
gives h2pd2lnΣ/dR2 < −1 in some radial regions. In such theexcitation torque densityΛex with a simple cutoff as
unstable regions, we have to use equation (14) instead of
equation (9). Λex = ΛWexKB for |R−Rp|>hp∆, (18)
The breakdown of equation (9) indicates that the flux (0 for R Rp 6hp∆.
| − |
mFJenotfuemqudaetpioonsit(e7d)bcyantnhoetptlraannestp.oInrtaarlleaolfstyhseteamn,ghuolawremveor-, The cut-off length hp∆ is determined so that the one-
the instability would enhance the angular momentum flux, sided torque T (= R∞pΛexdR) agrees with the result of the
linear theory for realistic discs (Takeuchi& Miyama 1998;
whichkeepsthegapmarginally stable.Theenhancementof R
Tanakaet al.2002;Muto & Inutsuka2009).Thenweobtain
FJ can beconsidered to be dueto an effectiveviscosity νeff
∆=1.3.
enhancedbytheinstability.Sincesuchaneffectiveviscosity
NotethattheWKBformulaisderivedfordiscswithno
restores equation (9),νeff in theunstableregion is given by
gap. Petrovich & Rafikov (2012) reported that the torque
νeff = −jFM+F4Jπ(∞R2)Σ−Ω R∞ΛddR′, (15) dsietnysibtyeciasuaselteorfedthbeysthhieftstoefepthgeraLdiniednbtlaodf trheesosnuarnfaccees.dFenor-
R
simplicity, however, we ignore this effect in the present pa-
where we use the relation dΩ/dR = 2Ω/R obtained from
− per. Hence, in our model, the excitation torque density
tinhsetemadargoifnνal,lyeqsutaatbiloenco(7n)digtiivoens. tFhuertehnehrmanocreed, baynguuslianrgmνeoff- Λex is simply proportional to the disc surface density at
R, Σ(R), and is independent of the surface density gra-
mentum fluxin theunstable region.
dient even for deep gaps. For a large planet with a mass
The Rossby wave instability may be important for the
gapformation(e.g.,Richard et al.2013;Zhu et al.2013;Lin of Mp/M∗ & (hp/Rp)3; furthermore, the non-linear effect
would not be negligible for wave excitation (Ward 1997;
2014).AswellastheRayleighcondition,theRossbywavein-
Miyoshi et al. 1999).Thisnon-lineareffect is also neglected
stabilityrelatestothediscrotation(Li et al.2000).Because
in our simple model.
it can occur before the Rayleigh condition is violated, how-
ever, the Rossby wave instability may suppress the surface
densitygradientmorethantheRayleighcondition.Forsim-
2.4.2 Case with wave propagation
plicity,weincludeonlytheRayleighconditioninthepresent
study.AfurtherdetailtreatmentincludingtheRossbywave When wave propagation is included, the angular momen-
instability should bedone in future works. tum deposition occurs at a different site from the wave
excitation and equation (18) is not valid. In this case,
the angular momentum deposition is also governed by the
2.4 Angular momentum deposition from a planet damping of the waves. Although the wave damping has
been examined in previous studies (e.g., Takeuchi et al.
Inthedisc–planetinteraction,aplanetexcitesdensitywaves
1996; Korycansky& Papaloizou 1996; Goodman & Rafikov
and the angular momenta of the waves are deposited on
2001),it isnot clear yethowthedensitywavesare damped
the disc through their damping. The angular momentum
in a disc with deep gaps. In the present study, therefore,
depositionrateΛdisdeterminedbythelaterprocess.Firstly, weadopt asimple model of angular momentum deposition,
we will consider the deposition rate Λd in the case with described below.
instantaneous wave damping. In this case, the deposition
Since thewaves are eventually damped in thedisc, the
rate is governed only by the wave excitation. Next, taking
one-sided torque (i.e., the total angular momentum of the
into account the wave propagation before damping, we will
wavesexcited at theouter discin unit time) is equalto the
model the deposition rate in a simple form.
total deposition rate in the steady state. That is,
∞ ∞
T = ΛexdR′= ΛddR′. (19)
2.4.1 Case with instantaneous wave damping ZRp ZRp
Under the assumption of instantaneous wave damping, the Usingtheone-sided torque,theangular momentumdeposi-
angular momentum deposition rate Λd(R) is equal to the tion rate can be expressed by
excitationtorquedensityΛex(R),whichistherateatwhich
a planet adds angular momenta to density waves per unit Λd =±Tf(R), (20)
radial distance at R.That is, wherethedistributionfunctionf(R)satisfies ∞f(R)dR=
Rp
Λd =Λex. (16) 1, andthesign is thesame as in equation (17R).Asasimple
model, we assume a distribution function f(R) given by
Atapositionfarfromtheplanet,theexcitationtorqueden-
1
sity is given by theWKB formula (e.g., Ward 1986) as f(R)= wd for xdhp− w2d <|R−Rp|<xdhp+ w2d,
ΛWexKB =±CπRp2Σ(cid:18)MMp∗(cid:19)2(RpΩKp)2(cid:18)R−RpRp(cid:19)4, (17) 0 otherwise. (21)
Formation of a disc gap induced by a planet 5
Inthissimplemodel,thenon-dimensionalparameterxd de- Rp given by equation (12). Dividing equation (25) by
termines the position of the angular momentum deposition 3πRp2νpΣ0(Rp)ΩKp and using equation (11), we obtain a
andtheparameterwd representstheradialwidthofthede- non-dimensional form:
position site. The waves propagate from the excitation site 1d2lns 1 ∞
twoavtehsepdreoppoasgiatitoenawsitaeyafrroomundthe|xp|l=anextd,.thSeindceeptohseitidoennssiittye (cid:18)1− 3 dx2 (cid:19)s=1− 3Zx λddx′, (27)
is farther from the planet than the excitation site. The pa- where λd is the non-dimensional angular momentum depo-
rameter xd should beconsistent with this condition. sition rate defined by
tionsIn(20t)heancdas(e21)wittohobwtaavine tphreopgaapgasttiorunc,tuwree wuisteh eeqquuaa-- λd = πRp2νpΛΣd0h(Rpp)ΩKp. (28)
tion(9).ItshouldbenotedthatT inequation(20)depends
The marginally stable condition can berewritten as
on the surface density distribution through the definition
of equation (19), because Λex is proportional to Σ. These d2lns
coupledequationsaresolvedasfollows.First,weobtainthe dx2 =−1. (29)
surfacedensitydistributionwithequation(9)foragivenT.
This equation is used instead of equation (27) in the
Next, we determine the corresponding mass of the planet
Rayleigh unstableregion.
from equation (19),using the obtained surface density.
The non-dimensional excitation torque density, λex, is
definedby
2.5 Local approximation and non-dimensional C
equations λex = πRp2νpΛΣe0x(hRpp)ΩKp =0±Kx4s(x) ffoorr |xx|>6∆∆.,
The typical width of a disc gap is comparable to the disc  | |
scale height and much smaller than the orbital radius of (30)

theplanet.Thusitisconvenienttousethelocalcoordinate
where thenon-dimensional parameter K is given by
defined by
x= R−hpRp. (22) K =(cid:18)MMp∗(cid:19)2(cid:18)Rhpp(cid:19)5α−1. (31)
Note that thesuffix p indicates thevalue at R=Rp. In the above, we use νp = αh2pΩKp. In our model, the pa-
rameter K is the only parameter that determines the gap
Weadoptalocalapproximationinwhichtermspropor-
structurefor theinstantaneous damping case.
tionaltohp/Rp andhigherordertermsareneglected.From
Inthecasewithinstantaneouswavedamping,theangu-
equations (2) and (3), the deviation in Ω from ΩK is given
by larmomentumdepositionrateisgivenbyλd=λex(eq.[16]).
In thecase with wave propagation, equation (20) gives
hpΩKpdlnΣ
Ω−ΩK = 2Rp dx , (23) λd(x)=T˜hpf(Rp+hpx), (32)
andis proportional tohp/Rp.Thus,thedisc angularveloc- where thenon-dimensional one-sided torque,T˜, is given by
iutnydΩeritshreeplolcaacledapbpyrotxhiemaantigounl,aranvdeltohceityspoefcitfihceapnlgaunleatrΩmKop- T˜=K ∞ Cs(x)dx. (33)
x4
mentum j also is given by Rp2ΩKp. As for the derivative Z∆
dΩ/dR,wecannotneglect thedeviationfrom theKeplerian Notethatthedepositionrateλdincludestwoparametersxd
value. Equation (4) yields andwd,inadditiontoK.Inthiscase,wesolveequation(27)
foragivenvalueofT˜.Then,wecanobtainKbysubstituting
dΩ = 3ΩKp 1 1d2lnΣ . (24) thesolutions(x)intoequation(33),asmentionedattheend
dR − 2Rp (cid:18) − 3 dx2 (cid:19) of the last subsection. In order to obtain the solution for a
certainK,weneedaniterationoftheaboveprocedurewith
Equation (9) can be rewritten in the local approxima-
trial values of T˜.
tion as
Theboundaryconditionsofequations(27)and(29)are
1d2lnΣ
R2ΩKpFM+3πRp2νpΣΩKp 1− 3 dx2 s=1, at x= . (34)
(cid:18) (cid:19) ±∞
∞ Under the local approximation, the surface density has a
=FJ( ) Λdhpdx′. (25)
∞ − symmetryofs(x)=s( x),sinceboththeabovebasicequa-
Zx −
tions and thedeposition rate are symmetric.
Because of the local approximation, equation (25) cannot
be applied for the wide gap formation. If the half width of
gap is narrower than about 1/3Rp, equation (25) would be
valid. 3 ESTIMATES OF GAP DEPTHS FOR
Hereweintroducethenon-dimensionalsurfacedensity, SIMPLE SITUATIONS
s, definedby
3.1 Case of the Keplerian discs
Σ
s= , (26) Before deriving the gap solution in our model described in
Σ0(Rp) Section 2, we examine the gaps for two simple situations.
where Σ0(Rp) is the unperturbed surface density at R = First,weconsideradiscwithKeplerianrotation,asassumed
6 K. D. Kanagawa
inpreviousstudies.NeglectingthedeviationindΩ/dRfrom Funget al.(2014)1.Inthezero-dimensionalmodel,themin-
theKeplerian(i.e.,thetermofd2lns/dx2)inequation(27), imum surface density is estimated from a balance between
we have the planetary torque and the viscous angular momentum
fluxoutsidethegap.Suchabalanceisalso seen from equa-
1 ∞
s=1 λddx′ (35) tion (40) (and eq. [27]). The first and second terms in the
− 3
Zx right-hand side of equation (40) correspond to the viscous
Herewealsoassumeinstantaneouswavedampingandadopt angular momentum flux outside the gap and the planetary
λd =λex (eq.[30]).Differentiatingequation(35),weobtain torque and the left-hand side is negligibly small for a large
K.
C
dlns K for x >∆, With their hydrodynamic simulations for K . 104,
= ±3x4 | | (36)
dx  DM13 derived a similar result2,
0 for x 6∆,
 | |
29 1
Hence, we obtain thesurface density in the Keplerian discs smin = 29+K = 1+0.034K. (42)
with theinstantaneous wave damping as
It is found that equations (41) and (42) are consistent with
C
exp K for x >∆, each other. Note that these minimum surface densities are
−9x3 | | much larger than that of the Keplerian disc (eq. [38]) for a
s(x)= (cid:18) | | (cid:19) (37)
exp C K for x 6∆, large K because equation (39) is not accepted in the Ke-
−9∆3 | | plerian solution. The wide-limit gaps assume that all the
(cid:18) (cid:19)
Using equation (31), C =0.798 and ∆=1.3, the minimum waves are excited in the bottom region with s ≃ smin, i.e.,
equation (39). In Sections 4 and 5, we will check whether
surface density,smin,is or not this assumption is valid, by comparing it with our
smin=exp 0.040α−1 Rp 5 Mp 2 . (38) one-dimensional solutions.
"− (cid:18)hp(cid:19) (cid:18)M∗(cid:19) #
This solution is almost the same as that in the previous
4 GAP STRUCTURE IN THE CASE WITH
one-dimensionalgapmodel(e.g.,Lubow & D’Angelo2006).
INSTANTANEOUS WAVE DAMPING
For a very large K, the Rayleigh condition is violated
andequations(37)and (38) areinvalid.Tanigawa & Ikoma 4.1 Linear solutions for shallow gaps
(2007) obtained the gap structure in Keplerian discs, in-
Herewepresentthenumericalsolutionofthegapinthecase
cluding the Rayleigh condition. Their solution is described
with instantaneous wave damping(i.e., λd =λex).
inAppendixA.InAppendixB.wealsoderivegapsolutions
First, we consider the case with a small K in equa-
inKepleriandiscs,takingintoaccountthewavepropagation
tion(27),inwhichλd isproportionaltoK.Thiscasecorre-
with thesimple model of equation (20) and (21).
spondstoashallowgaparoundasmallplanet.Since s 1
| − |
is small, it is useful to express thesolution as
3.2 Case of the wide-limit gap s=exp(Ky), (43)
Next,weconsiderasituationimpliedbythezero-dimension or s = 1+Ky. As seen in the next subsection, the former
analysisdonebyFunget al.(2014),whichassumesthatthe expression is better for an intermediate K ( 10). Substi-
∼
waveexcitationoccursonlyathegapbottom.Thisassump- tuting equation (43) into equation (27) with equation (16),
tionwouldbevalidifthegapbottomregioniswideenough. we can expand it into a power series of Ky. The first order
Hence,wecallthissituation’wide-limitgap’case.Sincethe terms give thelinear equation of y:
density waves are excited at thegap bottom with s≃smin, C
theone-sided torqueof equation (33) is simply given by for x >∆,
d2y 3x3 | |
T˜= 3C∆3Ksmin ≃0.121Ksmin. (39) dx2 −3y=∓3∆C| 3| otherwise, (44)
Usingequation (27),we can estimate smin of thewide-limit where the sign in the right-hand side is negative for x > 0
gap. The right-hand side of equation (27) can be rewritten and positive for x6 0. Equation (44) is an inhomogeneous
as 1 T˜/3 at x=0. In the left-hand side of equation (27), linear differential equation, and can be integrated with the
more−over, we can neglect the term d2lns/dx2 when a flat- boundary conditions of equation (34). We do not need to
bottomgapisassumed.Then,therelationbetweensminand take care of the Rayleigh condition in the shallow gaps. A
T˜ is obtained as detailedderivationofthelinearsolution isdescribedinAp-
T˜ pendixC.
smin=1 . (40) Fig.1ashowsy,whichcanbeconvertedintothesurface
− 3
density s by equation (43). In these shallow gaps, the gap
Equations (39) and (40) yield depth is almost the same as for the Keplerian case, though
1 ourmodel gives a smooth surface density distribution.
smin= . (41)
1+0.040K
ForalargeK,smingivenbyequation(41)isproportionalto 1 InthenotationofFungetal.(2014),Kisgivenbyq2/(α[h/r]5)
1/K.Thisresultagreeswiththezero-dimensionalmodelby 2 InthenotationofDM13,K isgivenbyM−1(Msh/Mp)2α−1.
Formation of a disc gap induced by a planet 7
a a
K=50
1
0.00
y s
linear
-0.02
eq.(41)
exact
Linear
-0.04 Kepler
Kepler 0.1
b
0.02 b K=200
1
Linear
0.1
dy s
eq.(41)
linear
0.00
dx Kepler 0.01 exact
0.001
Kepler
-0.02
-4 -3 -2 -1 0 1 2 3 4
-4 -3 -2 -1 0 1 2 3 4
x
−
R R
p
x =
hp Figure 2. SurfacedensitydistributionsforK=50(a)and200
(b). The red line is the exact solution (see text). The chain line
Figure 1. Linear solutionfory (a)and dy/dx (b). Thesurface isthe linearsolutiongiven byequation (43)and the dashedline
densityandangularvelocityaregivenbyyanddy/dxwithequa- is the solution for the Keplerian case (eq. [37]). The dotted line
tions (43) and (45), respectively. The dashed lineis the solution represents the minimum surface densities for the wide-limitgap
fortheKepleriandisc. givenbyequation(41).
Atregions far from theplanet,thesurface density per-
Fig.1bshowsthederivativeofywhichisrelatedto∆Ω,
turbation is rather small and the linear approximation is
as
valid.Thus,weadoptalinearsolutionat x >10.Notethat
∆Ω=Ω ΩK =KhpΩKp dy, (45) thislinearsolutionhasdifferentcoefficien|ts|forthehomoge-
− 2Rp dx neoustermsfromthoseinSection4.1(seeAppendixC).The
using equations (23) and (43). The absolute value of ∆Ω coefficientsofthehomogeneoussolution aregiventosatisfy
the boundary conditions of equation (34). At x 6 10, we
attainsamaximumat x 1.5.Thesecond-orderderivative
| |
| |≃ integrate equation (27) with the fourth-order Runge-Kutta
of y gives theshear, dΩ/dx, as
integrator.IntheRayleighunstableregion,thesurfaceden-
dΩ = dΩK 1 K d2y , (46) sityisgovernedbythemarginallystablecondition(eq.[29]),
dx dx − 3 dx2 instead of equation (27).
(cid:18) (cid:19)
Fig.2showsthesurfacedensitydistributionsoftheex-
as seen from equation (24). At x > 1.5, the shear mo-
| | actsolutionsforK =50(a)and200(b).Ifweassumeadisc
tion is enhanced compared to the Keplerian case, because
d2y/dx2 < 0. Since the shear motion causes viscous angu- withhp/Rp =0.05andα=10−3,thesecasescorrespondto
lar momentum transfer, this enhancement makes the sur- Mp =1/8MJ and1/4MJ respectively,whereMJisthemass
ofJupiter.Forcomparisons,theKepleriansolution(eq.[37])
facedensitygradientlesssteepcomparedwiththeKeplerian
and the linear solution with equation (43) are also plotted.
case, as shown in Fig. 1a.
ForK =50,thelinearsolutionalmostagreeswiththeexact
solution,whileitismuchdeeperthantheexactsolutionfor
K = 200. For K = 200, the Keplerian solution has a much
4.2 Nonlinear solutions for deep gaps
smaller smin than theexact solution.
Next we consider deep gaps around relatively large plan- Fig. 3 illustrates the angular velocities (a) and specific
ets. In this case, we numerically solve the non-linear equa- angularmomenta(b)fortheexactsolutionsforK =50and
tion (27) with theRayleigh condition. Wecall theobtained 200.SimilartothelinearsolutioninFig.1,theshearmotion
non-linear solution the “exact” solution. isenhancedat x &1.4.Thisenhancementoftheshearmo-
| |
8 K. D. Kanagawa
1.5 a
1.5 K=200
1.0 K=200 1
n
50 mi
0.5 s
Ω−Ω 0.5 10
K 0.0 λ =
c /R ex 0 s
p p Kepler
-0.5
-0.5 n
mi
-1.0 s
0
-1 1
=
-1.5
s
b -1.5
2.0
K=200 -4 -3 -2 -1 0 1 2 3 4
50
x
1.0
j −j K e ple r
p Figure 5. Excitation torque density givenbyequation (30)for
0.0
K =200. The two vertical lines indicate the positions with s=
c R
p p 10smin.
-1.0
In the wide-limit gap, it is assumed that the density waves
-2.0 are excited only at the gap bottom with s smin. Fig. 5
≃
shows the excitation torque density given by equation (18)
-4 -3 -2 -1 0 1 2 3 4 fortheexactsolutionwithK =200.Thistorquedensityin-
x dicatesthatthewavesareexcitedmainlyintheregionwith
s> 10smin. Thus the assumption of wave excitation at the
gap bottom is not valid in this case. Since wave excitation
Figure 3. (a) Deviation from Keplerian disc rotation and (b)
with a larger s increases the one-sided torque, this can ex-
specificangularmomentum,forK=50(dashed)and200(solid).
plainwhythegapoftheexactsolutionismuchdeeperthan
Thefilledcirclesindicatetheedgeofthemarginallystableregion
the wide-limit gap in Fig. 2. Note that this result for the
fortheRayleighcondition.
waveexcitationisobtainedinthecaseofinstantaneouswave
damping.Theeffectofthewavepropagationcanchangethe
1.4 gap width and the mode of wave excitation, as seen in the
next section.
1.2
Kepler
K
R 1.0
Ω
d 4.3 Effect of the Rayleigh condition
d 0.8
/
We further examine the effect of the Rayleigh condition on
Ω R 0.6 K=50
the gap structure. Fig. 6 shows the surface densities (a)
d d
0.4 andspecificangularmomenta(b)fortheexactsolutionand
the solution without the Rayleigh condition. The solution
0.2 200
without the Rayleigh condition has unstable regions with
0.0 dj/dx < 0 (i.e., 1.4 < x < 3.1). This comparison between
| |
0 1 2 3 4 these two solutions directly shows how the Rayleigh con-
x dition changes the gap structure. The Rayleigh condition
increasessminbyafactor6forK =200.Thisisbecausethe
marginal condition of d2lns/dx2 > 1 keeps the surface
Figure4. ShearofexactsolutionsforK=50(dashed)and200 −
density gradient less steep and makes the gap shallow.
(solid).
It can beconsidered that themarginally stable state is
maintained by νeff of equation (15). The non-dimensional
tion is also seen in Fig. 4. The enhancement promotes the form of equation (15) is given by
angular momentum transfer and makes the surface density
gradient less steep. For K =200, the Rayleigh condition is νeff = 3− x∞λddx′. (47)
ν 4s
violated.Intheunstableregion,themarginallystablecondi- R
tionfurtherreducesthesurfacedensitygradient.Thismakes Fig. 7 shows νeff in the unstable region for K = 200. The
the gap much shallower than for the Keplerian solution, as effective viscosity is twice as large as the original value at
seen in Fig. 2b. x = 1.8. This enhancement of the effective viscosity causes
Wealsoplottheminimumsurfacedensities,smin,ofthe theshallowing effect in Fig. 6a.
wide-limit gap (eq. [41]) in Fig. 2. Thewide-limit gap gives In Fig. 6a, we also plot the surface density distribu-
a much larger smin than the exact solution for K = 200. tion given by Tanigawa & Ikoma (2007)(hereafter TI07), in
Formation of a disc gap induced by a planet 9
a 1
1 K=200
0.1 D
s eq.(41) 0.1 M13
n
0.01 wT/I0 r7c smi eq.(41)
w/o rc 0.01
T
0.001 I0
K w 7
3.0 b Kepler 0.001 eplerw/o rc/ rc
w/o rc
1 10 100 1000
2.0
w/ rc K
1.0
j −jp 0.0 K e ple r Figure 8. Minimumsurfacedensities, smin, for the exact solu-
c R tion (red line) and the solution without the Rayleigh condition
p p
-1.0 (green line). The dashed line is smin in the Keplerian case. The
chain, dotted and solid lines denote smin given by the model of
TI07, the wide-limit gap (eq. [41]) and the empirical relation of
-2.0
DM13(eq.[42])respectively.
-3.0
-4 -3 -2 -1 0 1 2 3 4
Wealso showthattheKeplerian solution byTI07 does
x not satisfy the angular momentum conservation. The Kep-
lerian solution without the Rayleigh condition (eq. [37]) is
derived just from equation (35) (or eq. [27]), which is orig-
Figure 6. (a) Surface density distribution and (b) specific an-
inated from equation (9). In this solution, thus, the angu-
gular momentum distribution, for K = 200. The red line indi-
lar momentum conservation is satisfied. However, when the
cates the exact solution. The green line is the solution with-
out the Rayleigh condition (see text). The chain line in (a) Rayleighconditionisviolated,themarginalstablecondition
denotes the surface density distribution given by the model of (eq. [29]) is used instead of equation (37). Because of this,
Tanigawa&Ikoma(2007)(eq.[A3],TI07). thesurfacedensityattheflatbottomofTI07’ssolutiondoes
not satisfy equation (35) or theangular momentum conser-
vation, either. This violation is resolved in our formulation
2.2
becauseourexactsolutionalwayssatisfiesequation(27)out-
K=200 side theRayleigh unstableregion.3
2.0
1.8
4.4 Gap depth
ν
eff 1.6 Fig.8showstheminimumsurfacedensities,smin,asafunc-
ν tion of K for the exact solutions. For comparison, we also
p
1.4
plot smin for the solutions without the Rayleigh condition
and the Keplerian solutions. These solutions give deeper
1.2
gaps than the exact solution, similar to the result of Sec-
tion 4.2 . It is found that the shallowing effect due to the
1.0
Rayleigh condition becomes significant with an increase in
K. This is because the Rayleigh condition is violated more
0 1 2 3 4
x strongly for large K.
InFig.8,ontheotherhand,theexactsolutionismuch
deeperthanDM13’s resultsand thewide-limit gap, though
Figure7. Effectiveviscosityνeff oftheexactsolutionwithK= thelattertwocasesagreewellwitheachother.Themodelof
200. TI07 also gives much deeper gaps than DM13. These com-
parisons indicate that in the case with instantaneous wave
damping,ourexactsolutioncannotreproducethehydrody-
which the Rayleigh condition is taken into account (for de- namicsimulationsofDM13.Thisdifferenceinthegapdepth
tails, see Appendix A). Their model gives a shallower gap
than our exact solution. This is because a very steep sur-
facedensitygradientintheKepleriansolutionissuppressed 3 Byintroducingtheeffectiveviscosityofequation(47)andmul-
by the Rayleigh condition to a greater extent than in our tiplying the LHS of equation (27) by νeff/ν0, equation (27) is
model. recoveredintheRayleighunstableregion.
10 K. D. Kanagawa
a
from DM13 is likely to be due to the fact that the assump-
K=200
tion of the wide-limit gap is not satisfied in the case with 1
instantaneous wave damping (see Fig. 5). In the next sec- w =h
tion, we will see that the effect of wave propagation widens d p
the gap and makes the assumption of the wide-limit gap
valid.
x =4
s 0.1
eq.(41) d
3
5 EFFECT OF DENSITY WAVE
PROPAGATION 2
0.01
In this section, we consider the effect of wave propagation.
Wavepropagation changestheradialdistributionofthean- instant
gularmomentumdeposition.Asimplemodelofangularmo- b
x =4
mentum deposition rate altered by wave propagation is de- 4 d
scribed in Section 2.4.2. Using this simple model, we solve
equation(27)withtheRayleighconditioninthesimilarway
2 3
to the previous section. At the region far from the planet
(i.e., x > 10), we use the linear solution to equation (C1) 2
| | λ
with g(x)=0 in this case. ex 0
5.1 Gap structure for K =200 -2
Fig. 9 illustrates the surface densities (a) and the excita-
tion torque densities (b) of the exact solutions in the case -4
instant
withwavepropagation.Theangularmomentaoftheexcited
waves are deposited around x =xd in our model. A large -6 -4 -2 0 2 4 6
| |
xd indicates a long propagation length between the excita-
tion and the damping. The parameter K is set to 200. For x
anincreasingxd,thegapbecomeswiderandshallower.The
gap width is directly governed by the position of the angu-
Figure 9. Surface densities (a) and excitation torque densities
larmomentumdeposition.Forxd=3and4,thegapdepths
(b) in the case with the wave propagation for K = 200. The
areconsistentwiththewide-limitgap(andalsoDM13).For
green,blueandredlinesdenotethesolutionswithxd=2,3and
xd = 4, the density waves are excited mainly at the bot- 4 respectively. The parameter wd is set to hp. The gray dashed
tom region with s smin, as seen in Fig. 9b. Moreover, line is the surface density in the case with instantaneous wave
≃
for xd = 3, a major part of the wave excitation occurs at damping.Thedottedlinein(a)represents theminimumsurface
thebottom.Thatis,theassumptionofthewide-limitgapis densityforthewide-limitgapgivenbyequation(41),i.e.smin=
almostsatisfiedforthesolutionswithxd =3and4.Thisex- 0.109.
plainswhythegapdepthsareconsistentwiththewide-limit
gap for these large xd.
It is also valuable to compare the gap width with hy- casewithxd =3or4,theexcitationat x <xd contributes
| |
drodynamic simulations. DM13 performed a simulation for 55% or 78% of theone-sided torque,respectively.
thecaseofMp =1/4MJ (2Msh intheirnotation),α=10−3 InFig.10,wechecktheeffectofthewidthofthedepo-
and hp/Rp = 0.05. This case corresponds to K = 200. In sition site,wd,forxd=3andK =200. Itisfoundthatthe
thissimulation,theyfoundthatthegapwidthisabout6hp, width wd has only a small influenceon thegap structure.
assuming that these gap edges are located at the position We show that the deviation from the Keplerian rota-
with Σ=(1/3)Σ0(Rp) (i.e., s=1/3). If we adopt the same tion is also important in the case with wave propagation.
definitionofthegap edge,thegapwidthsofourexactsolu- InFig. 11, we plot thesolution with theKeplerian rotation
tions with xd = 3 and 4 are 6.1hp and 7.7hp, respectively. and our exact solution. The Keplerian solution is derived
Hence, if we take into account the wave propagation and from equation (35) with the angular momentum deposition
adopt xd = 3–4, our exact solution can almost reproduce model (eqs. [20] and [21]). When the Rayleigh condition is
both of thegap width and depth of thehydrodynamicsim- violated, the marginal stable condition (eq. [29]) is used. A
ulations by DM13, for K =200. detailderivationofthissolutionisdescribedinAppendixB.
Itshouldbealsonotedthat,forxd=2,thewaveexcita- IntheKeplerian solution of Fig. 11, theRayleigh condition
tionmainlyoccursat x >xd(80%oftheexcitationtorques is violated over thewhole region of the angular momentum
| |
comefrom thisregion). However,thedeposition siteshould deposition. Then the minimum surface density is given by
be farther from the planet than the excitation site because equation (B3), which is much larger than our solution and
thedensitywavespropagateawayfromtheplanet.Thus,the equation(41).Becauseequation(B3)doesnotsatisfyequa-
case with xd = 2 does not represent a realistic wave prop- tion(35),theKepleriansolutiondoesnotsatisfytheangular
agation. From now on, we judge that our simple model for momentum conservation, as pointed out in Section 4.3. On
the wave propagation is valid if more than half of the one- theotherhand,inthezero-dimensionanalysisbyFunget al.
sided torque arises from the excitation at x < xd. In the (2014) (or in eq. [41]), smin is estimated from a balance be-
| |