On corotation torques, horseshoe drag and the possibility of sustained stalled or outward protoplanetary migration

0.49 MB·English

by  S. -J. Paardekooper

Mon.Not.R.Astron.Soc.000,1–13(2008) Printed21January2009 (MNLATEXstylefilev2.2)
On corotation torques, horseshoe drag and the possibility
of sustained stalled or outward protoplanetary migration
9 S.-J. Paardekooper⋆ and J.C.B. Papaloizou
0
0 DAMTP, Universityof Cambridge, Wilberforce Road, Cambridge CB3 0WA, United Kingdom
2

n
Draftversion21January2009
a
J

5
ABSTRACT
1
We study the torque on low mass protoplanets on fixed circular orbits, embedded in

] a protoplanetary disc in the isothermal limit. We consider a wide range of surface
P density distributions including cases where the surface density increases smoothly
E outwards.We performbothlineardisc responsecalculationsandnonlinearnumerical
. simulations. We consider a large range of viscosities, including the inviscid limit, as
h
well as a range of protoplanet mass ratios, with special emphasis on the coorbital
p
- region and the corotationtorque acting between disc and protoplanet.
o For low mass protoplanets and large viscosity the corotation torque behaves as
r expected from linear theory. However, when the viscosity becomes small enough to
t
s enable horseshoe turns to occur, the linear corotation torque exists only temporarily
a
after insertion of a planet into the disc, being replaced by the horseshoe drag first
[
discussed by Ward. This happens after a time that is equal to the horseshoe libra-

1 tion period reduced by a factor amounting to about twice the disc aspect ratio. This
v torque scales with the radial gradient of specific vorticity, as does the linear torque,
5 but we find it to be many times larger. If the viscosity is large enough for viscous
6
diffusion across the coorbital region to occur within a libration period, we find that
2
the horseshoe drag may be sustained. If not, the corotation torque saturates leaving
2
only the linear Lindblad torques.As the magnitude of the nonlinear coorbitaltorque
.
1 (horseshoe drag) is always found to be larger than the linear torque, we find that
0 the sign of the total torque may change even for for mildly positive surface density
9
gradients. In combination with a kinematic viscosity large enough to keep the torque
0
fromsaturating,strongsustaineddeviationsfromlineartheoryandoutwardorstalled
:
v migration may occur in such cases.
i
X Key words: planetary systems: formation – planets and satellites: formation.
r
a
1 INTRODUCTION scale (Lin & Papaloizou 1986). For standard disc param-
eters, planets more massive than Jupiter migrate in this
A planet embedded in a gaseous disc is subject to a net
Type II regime (Crida & Morbidelli 2007). Intermediate
torquethat givesrise to orbital evolution. Sincethediscov-
mass planets, comparable to Saturn, embedded in mas-
ery of the first hot Jupiter (Mayor & Queloz 1995), planet
sive discs may undergo runaway or Type III migration
migration due to interaction with the protoplanetary disc
(Masset & Papaloizou 2003),whichmaybedirectedinward
has become a necessary ingredient of planet formation the-
aswell as outward (Peplin´ski et al. 2008).In thispaper,we
ory.Aconsiderableamountofanalyticalandnumericalwork
focus on low mass planets, with masses typical up to a few
has gone into understanding the disc-planet interactions
timesthemass oftheEarth (M⊕).Lowmass planetsexcite
leading to planet migration (see Papaloizou et al. 2007, for
linearwavesinthedisc,theactionofwhichleadstoTypeI
an overview).
migration (Ward1997).Linear,semi-analyticalcalculations
Three modes of migration can be distinguished, in
haveresultedinawidelyusedtorqueformula(Tanaka et al.
most cases leading to migration towards the central star.
2002, hereafter TTW02) for isothermal discs.
High mass planets, for which the Hill sphere is larger
than the disc scale height, open up deep gaps in the Thetotaltorqueonanembeddedplanetcanbedecom-
disc, after which migration proceeds on a viscous time posedintorquesduetowaves,whicharegeneratedatLind-
blad resonances and propagate away from the planet, and
corotation torques, generated near the orbit of the planet,
⋆ E-mail:[email protected] where material on average corotates with the planet (see
2 S.-J. Paardekooper and J.C.B. Papaloizou
Goldreich & Tremaine1979).Forlowmassplanets,thewave The plan of the paper is as follows. We start in Sec-
torqueisusuallythoughttodominatethetotaltorque,and, tion 2 by reviewing the model used throughout the pa-
beingrelativelyinsensitivetobackgroundgradients,leadto per. In Section 3, we review and discuss linear corotation
inward migration. torque estimates as well as torque estimates based on the
Although Type I migration was thought to be math- horseshoedragexperiencedbymaterialexecutinghorseshoe
ematically well-understood, it nevertheless posed a serious turns. These are shown to be distinct phenomena with dif-
problem for planet formation. Applyingthe torque formula ferentdependenciesonthephysicalvariables,thehorseshoe
from TTW02 (see also Korycansky & Pollack 1993, here- dragbeingessentiallynonlineareventhoughtheassociated
after KP93) to planets of a few M⊕ embedded in a typi- torquescalesinthesamewayasthelinearone.Wegoonto
cal disc resulted in inward migration time scales that are perform detailed linear calculations in Section 4, and then
much shorter than the lifetime of the disc. In other words, compare thelinear and non-linear(horseshoe drag) torques
these planets would all migrate to very small orbital radii, inSection5.InSection6,wepresenttheresultsofnumerical
oreven intothecentral star(Ward 1997). Sincegaseous gi- hydrodynamical simulations, torques are also compared to
ant planets are thought to form around a solid core that linearandhorseshoedragestimatesconfirmingtheviewout-
is in this mass range, Type I migration theory essentially lined above.Thehorseshoe dragis foundtobesignificantly
predictsthatthereshouldbenoplanetsat largeradii. This larger than the linear corotation torque and potentially a
hasledtoinvestigationsonhowtoslowdownorstopTypeI very important contributor to the total torque for disc sur-
migration,forexamplethroughtheactionofmagneticfields facedensitiesthat increaseevenmildly outwards.Wegoon
(Terquem2003;Nelson & Papaloizou2004)orsharpsurface toconsiderthelongtermtorqueevolutionandsaturationas
density gradients (Masset et al. 2006). a function of viscosity finding corotation torque saturation
More recently, efforts have been made to relax the inlowviscositycasesandsustainedcorotationtorqueswhen
isothermal assumption that was made in previous works, theviscosityislargeenoughtoresupplyangularmomentum
and to properly account for the energy budget of the disc. tothecoorbitalregion.WediscussourresultsfurtherinSec-
Since radiation is the dominant cooling agent, radiation- tion 7. Finally we present a short summary, together with
hydrodynamical simulations are needed. The outcome of some concluding remarks, in Section 8. In addition we give
thesesimulationsweresurprising:Paardekooper & Mellema an analysis and discussion of the time development of the
(2006a) found that the torque on the planet was positive, linear corotation torque acting on a planet after immersion
leading to outward migration, whenever the opacity of the intoadisc,demonstratingthatthecharacteristictimeisthe
disc at the location of the planet was high enough to make orbital period as was found in the numerical simulations.
cooling inefficient. In Paardekooper & Mellema (2008) it
wasrecognisedthatthispositivetorquewasacorotationef-
fect,relatedtoaradialentropygradientintheunperturbed
2 BASIC EQUATIONS
disc. Baruteau & Masset (2008), through a linear analysis
of the corotation torque, suggested that the linear corota- The evolution of a gaseous disc is governed by the Navier-
tion torque in the presence of an entropy gradient can be Stokes equations. We will work in a cylindrical coordinate
strong enough to overcome the negative wave torque. How- polar frame (r,ϕ), centred on the central star, and we in-
ever, Paardekooper & Papaloizou (2008) showed that this tegrate the equations of motion vertically to obtain a two-
linearcontributionissmall,andthatagenuinelynon-linear dimensional problem. We are then left with the continuity
effectisresponsibleforthechangeofsignofthetotaltorque. equation and two equations of motion:
In this paper, we take a step back and reanalyse the
∂Σ
isothermal case, in a two-dimensional set-up. The simplic- +∇·Σv=0 (1)
∂t
ityofthismodelallowsustoperformnumericalsimulations
∂v ∇p
with high enough resolution investigate in detail non-linear +(v·∇)v=− −∇Φ+fvisc, (2)
∂t Σ
effects on the corotation torque and to run them for long
enough to study its possible saturation. We will show that where Σ is the surface density, v=(v,rΩ)T is the velocity,
linear theory only applies for short times on the order of pistheverticallyintegratedpressure,Φisthegravitational
a few orbits when a protoplanet is inserted into a disc for potentialandfvisc representstheviscousforce.Weusealo-
which the viscosity is not too large and there is a non-zero callyisothermalequationofstate,p=c2Σ,wherethesound
s
corotationtorque.Wealsoshowthatitisanon-lineareffect speedcs maybeaprescribedfunction ofradius.Wetakec2s
associated with horseshoe bends (Ward 1991) that results tobeapowerlawwithindex−β,whichmakes−β basically
in a departure from linear theory. This departure is found thepowerlawindexoftheradialtemperatureprofile.Writ-
to result in torquesthat can havea much larger magnitude ingcs =HΩK,whereH istheverticalpressurescaleheight
thanthelinearones. Asthetorquescales with thegradient and ΩK is the Keplerian angular velocity, we usually adopt
ofspecificvorticity,thismayproducenoticeableeffectsthat either a sound speed that gives rise to a constant aspect
may even lead to stalled or outward migration, when this ratio h=H/r (or, equivalently,β =1), or apurely isother-
gradient is relatively large as for example occurs when the mal disc with constant cs (β = 0). In the latter case, we
surface density increases gently outwards. For these effects usually quotetheaspect ratio at thelocation of theplanet,
to be sustained the viscosity must be large enough to pre- hp = cs/(rpΩp). The initial, or for linear calculations the
venttorquesaturation.Viscositiescomparabletothoseoften background,surfacedensityistakentobeapowerlawwith
assumed in protoplanetary disc modelling are found to be index −α. Thus Σ = Σp(r/rp)−α, where Σp is the surface
large enough to enable these non linear corotation torques density at the location of the planet. The gravitational po-
to besustained for protoplanets in theEarth mass range. tentialΦisthesumofthepotentialduetothecentralstar,
Corotation torques and horseshoe drag 3
the planet’s potential Φp, and an indirect term that arises where Φ′m is the m-th Fourier component of the planet’s
duethetheacceleration of thecoordinate frame centredon potential, m is the azimuthal mode number and ω is the
thestar. ForΦp, we use a softened point mass potential: flow vorticity, being equal to Ω/2 in a Keplerian disc. All
quantities in equation (4) should be evaluated at corota-
Φp =− r2−2rrpcos(GϕM−pϕp)+rp2+b2rp2, (3) ttoionth(erra=dirapl)g,rwadhiiecnhtminakspesectifihce vtootratilcittoyr,qourevporrotepnosrittiyo,nianl
whereGpisthegravitationalconstant,Mp isthemassofthe the unperturbed flow there. Of course unless the softening
planet, and (rp,ϕp) are the coordinates of the planet. We parameterisratherlarge, theperturbingpotentialdoesnot
willalsouseq toindicatethemassratioMp/M∗,whereM∗ vary slowly on the corotation circle and so (4) cannot be
is the mass of the central star. The softening parameter b used as described above. However, it becomes valid if the
should be a sizeable fraction of h, in order to approximate perturbing potential Φ′ is replaced by the generalised po-
m
theresultofappropriateverticalaveragingofthepotential. tential obtained by adding the enthalpy perturbation to it
The exact form of theviscous terms can be found else- (seee.g.Lai & Zhang2006,andtheappendix).Butthenthe
where (e.g. D’Angelo et al. 2002). We take the kinematic linearresponseequationsneedtobesolvedinordertodeter-
viscosityν tobeapowerlawinradius,suchthattheinitial minetheenthalpyperturbationinordertoevaluatetorques
surface density profile is a stationary solution. It is easy to using the modified form of (4), which are nonetheless still
see that the required viscosity law is ν ∝rα−1/2. This way, proportional tothe gradient of specific vorticity.
we can neglect a global radial velocity field in the model,
which greatly simplifies the analysis. Values quoted in the
3.2 The total linear corotation torque in the limit
text refer to ν at the location of the planet, and will be in
unitsofrp2Ωp,whereΩpistheangularvelocityoftheplanet. of zero softening
To obtain the total torque, one needs to sum the contribu-
tionsfrom allvaluesofm(seeWard1989).ForaKeplerian
3 COROTATION TORQUES disc with zero softening, the total torque has been calcu-
lated,byfindingthelinearresponsenumerically,tobegiven
The wave (or Lindblad) torque exerted by a low-mass
by(see TTW02)
planet on a gaseous disc is relatively well-understood
(Goldreich & Tremaine 1979, TTW02). When the Hill 3 q2
sphere of the planet is much smaller than the scale height Γc,lin=1.36 2 −α h2Σprp4Ωp. (5)
„ « p
of the disc, the density waves are excited at Lindblad res-
Fora3Ddiscitisfoundthatthesameexpressionholdsbut
onances are linear and the resulting torque can be calcu-
with thenumerical coefficient beingreplaced by0.632.
lated by performing a linear response calculation and then
An important issue is the time required to develop the
summing the contributions arising from individual Fourier
linear corotation torque. Being linear this should not de-
components (TTW02). This torque usually leads to inward
pend on the mass of the protoplanet but only on intrinsic
migration, andhasbeenthoughttodominatein theTypeI
discparameters.Furthermorewhenaprotoplanetisinserted
regime of low mass planets. There are two approaches that
into a disc, there has to be an initial linear phase and for
can be followed in order to obtain corotation torques. The
sufficiently low protoplanet mass, the full linear response
first, which we will refer to as the linear estimate, is based
should develop.Asit is somewhat involved,we relegate the
onalinearresponsecalculation. Thesecondisafundamen-
discussion of the time development of the linear corotation
tallydifferentapproachbasedonadirecttorquecalculation
torque to the Appendix. From this discussion, the charac-
madeaftermakingsomeassumptionsaboutthetrajectories
teristic development time is expected to be on the order of
of the disc fluid elements, which we refer to as a horseshoe
theorbitalperiod.Wenowgoontoconsiderthesubsequent
drag calculation. The first approach applies at early times
development of thecorotational flow.
after a protoplanet is inserted into a disc as, at least for a
lowmassprotoplanet,theresponsemustfirstofallbelinear.
Thesecondapproachappliesatlatertimesafterthestream- 3.3 The horseshoe drag
line pattern has adjusted to thepresence of the planet.
Ward(1991)derivedaformulaforthecorotationtorqueus-
ing a fundamentally different approach. The argument was
3.1 Formula for the linear corotation torque for based on the expected form of the gas streamlines. These
each azimuthal mode number can be classified into four groups as viewed in a reference
frame corotating with the planet (see Fig. 1). The first
In order to obtain linear estimates for corotation torques,
pass the protoplanet interior to the coorbital region and
the equations are linearised (see Section 4) and Fourier-
the second pass it exterior to the coorbital region. These
decomposed in azimuth. The resulting single second-order
groups extend to large distances from theprotoplanet. The
equation (see Goldreich & Tremaine 1979) governing the
other two groups consist of material in the corotation re-
disc response can be solved to find the corotation torque
gion on horseshoe orbits that execute turns close to the
acting on the planet. Using an approximation scheme that
planet.Thethirdgroupapproachandleavetheprotoplanet
assumes the perturbing potential varies on a scale signifi-
from itsleadingsidewhile thefourthdosofrom itstrailing
cantly longer then the scale height, Goldreich & Tremaine
side. These groups are separated by two separatrices (see
(1979) findthis torquegiven by
e.g.. Masset et al. 2006; Paardekooper & Papaloizou 2008,
mπ2Φ′2 d Σ formorediscussionofthisaspect).Forlowmassprotoplan-
Γc,m =−2dΩ/dmr dr ω , (4) ets, most of the corotation torque is produced in a region
„ «
4 S.-J. Paardekooper and J.C.B. Papaloizou
Here j =rvϕ is the specific angular momentum and jp is j
evaluated at theorbital radius of the protoplanet.
Wenowfollow Ward(1991)andassumethestreamline
Γ pattern is symmetric on the leading and trailing sides such
3.20 2
that a streamline entering on Γ2 can be identified with a
III corresponding streamline leaving on Γ1 at the same radial
location. But note that on account of the horseshoe turns,
I II these streamlines will enter R at different radii. Since, in
3.15 the barotropic case, potential vorticity or vortensity Σ/ω
φ is conserved along streamlines, and these streamlines are
disconnected,thiswill differon them.Accordingly wewrite
thetorqueas
IV
3.10
Σ
Γ Γc,hs=− ∆ ω ω(j−jp)(Ω−Ωp)rdr. (7)
1 ZΓ2 „ «
Here
3.05
0.96 0.98 1.00 1.02 1.04 ∆ Σ =2 Σ − Σp (8)
r „ω« „ω ωp«
Figure 1.Schematic overview of the horseshoeregion, withthe isthevortensitydifferenceonthecorrespondingstreamlines
planet denoted by the black circle. Solid grey lines indicate ap- which have been assumed to enter R at the same radial
proximate streamlines, withthe separatrices inblack. The sepa- distancefromtheplanetonoppositesides,ωp =Ωp/2isthe
ratices divide the region into four parts, labelled by Roman nu- vorticityat theplanetsorbital radiusand wehaveassumed
merals: I and II, the inner and the outer disc, respectively, and thevortensitytobeanevenfunctionofradialdistancefrom
III and IV the leadingand trailingpartof the horseshoe region,
theprotoplanet.
respectively.Twohorizontaldashedlinesindicatetheboundaries
The above integral is easily done if one adopts a first
ofthetorquecalculation,Γ1 andΓ2.Notethatinthisfigure,the
order Taylor expansion for the quantities in brackets and
locationofΓ1andΓ2ischosenarbitrarily.Inthemodel,theyare
chosen to be far enough from the planet so that the corotation integratesfromrp−xsrp torp+xsrp wherethedimension-
torqueisdeterminedwithin. less width ofthehorseshoe region is xs.Oneobtains (Ward
1991):
3 3
close to the planet with a length scale expected to be the Γc,hs= 4 2 −α x4sΣprp4Ω2p. (9)
larger of brp or H. It is important to note that it takes „ «
some time for the streamline structure described above to Note that as is apparent from the above discussion and
develop on this scale after insertion of a protoplanet into that given in the Appendix and also the results of nu-
a disc. As this structure represents a finite deviation from merical simulations to be presented later, the horse-
theinitialform,thistimedependsonthemassoftheproto- shoe drag, given by equation (9), occurs as a non
planet.Thecalculation ofthehorseshoe dragappliestothe linear effect that has no counterpart in linear theory
situation when this streamline structure has developed on (see also Paardekooper & Papaloizou 2008, in addition to
this scale. It is important to note that the time required is Paardekooper & Papaloizou 2009). We also note that nu-
shorterthanthatrequiredtodevelopitonascalecompara- mericalhydrodynamicalcalculationsnecessarilyhaveb>0,
ble to rp which is when issues of torque saturation need to and, two-dimensional simulations, generally adopt b ≈ h
beconsidered. to account for three-dimensional effects in an approximate
FollowingWard(1991)weconsiderthetorqueproduced way. We will see that this strongly affects the width of the
bymaterial on streamlines undergoing horseshoe turns. We horseshoe region xs. A detailed analysis of the horseshoe
consider a region R interior to the two separatrices sepa- region is presented in Paardekooper & Papaloizou (2009).
rating these from the first and fourth group of streamlines Here, we adopt a simple estimate for xs that has proved
that is bounded by two lines of constant ϕ, Γ1 and Γ2 on to be reasonable for smoothing lengths comparable to hp
thetrailingandleadingsidesoftheprotoplanetrespectively (Paardekooper & Papaloizou 2008):
(seeFig.1).Theseboundariesaresupposedtobesufficiently
2q
farfromtheprotoplanetthatthecorotationtorqueisdeter- xs= . (10)
3b
mined within. Assuming a steady state this torque may be r
obtainedbyconsideringtheconservationofangularmomen- Ingeneral, thisdependenceistobeexpectedwhenthesep-
tum within R written in the form aratrix streamline passesthroughthelocation oftheplanet
withbpossiblybeingreplacedbyhp forsmallsoftening(see
alsoMasset et al.2006).Wewillcomparethehorseshoedrag
Γc,hs = Σ ∂Φ∂ϕGp rdϕdr= resulting from xs to the linear corotation torque in Section
ZR „ « 5. Note that, since Γc,hs is proportional to x4s, the horse-
Γ2 shoedragwould thenbeproportionaltoq2/h2 in thesmall
− Σ(j−jp)(Ω−Ωp)rdr . (6) softeningcaseandq2/b2 inthelargesofteningpcase,just as
»Z –Γ1
Corotation torques and horseshoe drag 5
Γc,lin. This means that the dependence of the total torque 4
on q, b and hp would not be a good indication of linearity,
since both the linear and the non-linear corotation torque Re(W)
scale in the same way. The only ways to distinguish these
Im(W)
are through their magnitudes, and bytheir time evolution.
2
Thelinearcorotation torqueissetupinapproximately
anorbitaltimescale(seetheAppendix),similartothewave
torque. This means that, when thebackground state of the
disc is stationary, any evolution in the torque after a few
dynamicaltimescalesisduetonon-lineareffects.Thefinite W 0
width of the horseshoe region gives rise to a libration time
scale
8π
τlib= 3xsΩ−p1, (11)
-2
which is basically the time it takes for a fluid element at
orbital radius rp(1+xs) to complete two orbits in a frame
corotating with theplanet. Therefore, this is thetime scale
on which the corotation torque will saturate due to phase
mixing(seeWard2007),unlesssomeformofviscosityoper- -4
ates on smaller time scales (Masset 2002). Note that satu- 0.6 0.8 1.0 1.2 1.4 1.6 1.8
r
rationisanon-linearprocess(Ogilvie & Lubow2003),since
it hinges on the finitewidth of thecorotation region.
Figure 2. Real (solid) and imaginary (dotted) part of the en-
We will see that there are basically three time scales
thalpy perturbation W,for an isothermal,constant surfaceden-
in the problem, two of which are closely related to non- sitydiscwithhp=10−3/2 andasofteningparameter asusedin
linearity. First, there is the orbital time scale, on which all KP93(b=10−4)plotted asafunctionofradiusinunitsofrp.
linear torques are set up. This is the shortest time scale in
the problem. The longest time scale is τlib, on which satu-
rationoperates.Athirdtimescalegovernsthedevelopment
m-th Fourier coefficient. Here the velocity perturbation is
ofthehorseshoedrag,whichwewillseetakesafractionofa v′ = (v′ ,u′ ), W′ = c2Σ′ /Σ is the enthalpy perturba-
lainbdraltiieosnintimbeetw(seeeenatlhsoe tPimaaerdseckaoleopfoerrt&hePlainpeaaloriztoourqu20e0d8e)-, timon,σ¯ =mm(Ωm−Ωp)m,BissthemsecondOortconstant,andΦ′m
isthem-thFouriercomponentof Φp,beingareal quantity.
velopment and the saturation time scale.
The term proportional to β is not present in the equations
ofKP93,becausetheyconsideredastrictlybarotropicequa-
tionofstate.Eliminatingu ,weobtainasystemofordinary
m
4 LINEAR CALCULATIONS
differential equations (KP93):
Since we are interested in departures from linear theory, it dv′ 2mB v′
is necessary to first firmly establish what linear theory pre- m =− 1−α− m +
dr σ¯ r
dicts.The2DlinearcalculationsofKP93andTTW02usea „ «
m2 σ¯ im2Φ′
vanishinglysmallvalueforb,andtheseresultsaretherefore i − W′ + m (15)
notdirectlycomparabletoourhydrodynamicalsimulations. „σ¯r2 c2s« m σ¯r2
Our customised linear calculations, as briefly outlined be- dW′ σ¯2−κ2
m =−i v′ −
low, with b comparable to h, can be directly compared to dr σ¯ m
„ «
hydrodynamicalsimulations. 2mΩ(W′ +Φ′ ) dΦ′ βW′
m m − m − m, (16)
σ¯r dr r
4.1 Governing equations where κ2 =4BΩ is thesquare of theepicyclic frequency.
Linearisingequations(1)and(2),withfvisc=0,andassum- We have solved equations (15) and (16) using a sixth
ing a Fourier decomposition such that perturbation quanti- order Runge-Kutta method and outgoing wave boundary
tiesaretherealpartofasumofterms∝exp(imϕ−imΩpt), conditions(seeKP93).InFig. 2,weshowtheresultingWm′
m=0,1,2....,oneobtainsthefollowing systemofequations for m=10, for an isothermal, constant surface density disc
(Goldreich & Tremaine1979) for thecorresponding Fourier with hp = 10−3/2 and a very low value of the softening
coefficients, parameter b = 10−4. The same case was shown in KP93
(theirfigure 2), and the results agree verywell.
iσ¯vm′ −2Ωu′m+ dWdrm′ + ddΦr′m + βWrm′ =0, (12) The imaginary part of Wm′ is directly related to the
torquedensity:
imW′ imΦ′
iσ¯u′m+2Bvm′ + r m + r m =0, (13) dΓm =πmΦ′ ΣIm(Wm′ ), (17)
iσ¯Wm′ + dvm′ +(1−α)vm′ + imu′m =0, (14) dr m c2s
c2 dr r r
s and the total torque on the planet Γ can be found by in-
m
where primes denote perturbation quantities, and the sub- tegratingoverthewholedisc.Goldreich & Tremaine(1979)
script m, being the azimuthal mode number, indicates the showed thattheLindblad torqueiscarried away bydensity
6 S.-J. Paardekooper and J.C.B. Papaloizou
b (h) b (h)
0.01 p 0.10 1.00 0.01 p 0.10 1.00
0
1000
-1000 TTW02 (3D)
0
-2000
TTW02 (2D)
q2 q2
Γ/ Γ/
-3000
-1000
KP93
-4000 h const h const α=1/2 Γ
isothermal isothermal -2000 α=-1/2 ΓC
h const (KP93) Γ
L
TTW02 (3D)
-5000
0.0001 0.0010 0.0100 0.0001 0.0010 0.0100
b (r) b (r)
p p
Figure 3. Total torque on the planet, obtained from linear cal- Figure4.Lindblad(dash-dottedcurves)andcorotation(dashed
culations, as a function of the softening parameter b. The sur- curves) torques obtained from an isothermal (hp = 0.05) linear
face density is proportional to r−3/2 (α = 3/2), which means calculationasafunctionofsofteningparameterbfortwodifferent
that the corotation torque vanishes in the barotropic case. The density profiles. Black curves: α = 1/2, grey curves: α = −1/2.
disc is either fully isothermal with hp = 0.05 (dashed line) or Thesolidcurvesdenotethetotaltorques,andthehorizontaldot-
has a constant aspect ratio h = 0.05 (solid line). Calculations tedlinesthe3DresultsofTTW02.
similar to those of KP93 with constant h are indicated by the
dash-dottedline.Alsoshownaretheresultsoflinearcalculations
in 2D by KP93 (for constant h) and the isothermal results of mal (h constant; solid line) and locally isothermal without
TTW02,for2Daswellas3D.Both2Dstudiesusedaverysmall the term proportional to β in equation (12) (dash-dotted
(KP93)orvanishing(TTW02)softeninglength.Filledcirclesde- line). The latter case is similar to the one considered in
note the torques measured fromhydrodynamical simulations for KP93, and approaches the result of KP93 for small soft-
aq=1.26·10−5 planetinadiscwithconstanth,anddiamonds ening.Similarly,thestrictlyisothermal caseapproaches the
theresultsforisothermalsimulationswithhp=0.05. 2D result of TTW02 for small softening. For appropriate
values of b≈h, differences between the strictly and locally
isothermal models can be as large as 30%. In the remain-
waves, resulting in an angular momentum flux
der of this paper, we will restrict ourselves to the strictly
F =− mπrΣ Im(W′ ) d (Re(W′ )+Φ′ )− isothermal equation of state, since the horseshoe drag has
m σ¯2−κ2 m dr m m been analysed for barotropic fluids only (Ward 1991)2. For
»
dIm(W′ ) this case, a smoothing length of b = 0.025, corresponding
(Re(W′ )+Φ′ ) m . (18)
m m dr to 0.5 h, results in the total torque being equal to the 3D
–
calculations of TTW02.
Therefore, the Lindblad torque is given by ∆F =
m We now introduce a corotation torque in the calcula-
Fm(rout)−Fm(rin),whererin androut denotetheinnerand tionsbyconsideringdifferentvaluesofα.InFig.4,weshow
the outer radius of the disc, respectively. We can therefore
thelinear Lindblad and corotation torques for α=1/2 and
calculate thecorotation torque as (KP93):
α = −1/2. Note that in the latter case, the surface den-
Γc,m =Γm−∆Fm, (19) sitygradientispositive,whichmaybeunrealisticexceptfor
special locations in the disc (see for example Masset et al.
and the total torques can be found by summation over all 2006),butitservesasagood exampleofacasewithstrong
m. corotation torques. For α = 1/2, the corotation torque is
much smaller than the Lindblad torque, and is increased
by a factor of 3 going from b = h to b = 0, similar to the
4.2 Results
Lindbladtorque.Therefore,thecorotation torqueisasmall
Westartbyconsideringadiscwithconstantspecificvortic- fractionoftheLindbladtorqueforallvaluesofbconsidered
ity, which means that the corotation torque vanishes. This here.Thatisnotthecaseforα=−1/2,whereforsmallsoft-
allowsustolookinsomemoredetailattheLindbladtorque
alone. In Fig. 3, we show the total torque1 for three differ-
ent cases: strictly isothermal (dashed line), locally isother- 2 In Paardekooper &Papaloizou (2008), horseshoe drag was
studied in adiabatic flows, where entropy (but not vortensity)
isconserved.Foralocallyisothermalequationofstate,neitherof
1 In all figures, the torque is given in units of Σprp4Ω2p, and is these is conserved, with the consequence that it is unclear what
dividedbyq2 tomakeitindependent ofthemassoftheplanet. thehorseshoedragisinthiscase.
Corotation torques and horseshoe drag 7
b (h)
p
0.2 0.4 0.6 0.8 1.0
model
h=0.1
100 Γ
h=0.05 800 ΓC
h=0.025 HS
Γ 3D
C
10 600
q2
/C
Γ
q2
Γ/
400
1
200
0.1 1.0
b (r)
p
0.01 0.02 0.03 0.04 0.05
b (r)
Figure 5.Corotationtorque, obtained fromalinearcalculation p
with α = 1/2 and hp = 0.025 (solid line), hp = 0.05 (dashed Figure 6.Corotation torque, obtainedfromalinearcalculation
lpirneed),icatinodnhopfe=qu0a.1tio(nda2s4h,-dvoatlitdedfolrinhe)≪. Tbh≪e d1o.tted linegives the withα=1/2andhp=0.05(solidline),togetherwiththehorse-
shoeorbitdrag(dashedline)foroursimpleestimateofthewidth
of the horseshoe region. Also shown is the 3D result of TTW02
ening parameters the corotation torque is almost the same (dash-dotted line).
as the Lindblad torque. The sign of the total torque will
changeforα<−1/2,whichisexpectedfromthe2Dresults
of KP93 and TTW02. For 3D calculations, the situation is
fromlargem∼1/b,wereplacethesumbyanintegral,thus
different(TTW02),whichisillustratedbytheresultsinFig.
4forlarger softening.ForallvaluesofαconsideredinFigs.
3 and 4 we find that for a softening length of b= 0.5 h we ∞ 1 2 ∞ 1 1 2
canmatchour2Dcalculationswiththe3DworkofTTW02. mK0(mb)2 → b xK0(x)2dx= 2 b . (23)
m=1 „ « Z0 „ «
X
Thuswe obtain
4.3 The limit of large softening
∞ 3 4q2
Wenowrelateournumericalresultstothepredictionsmade Γc,m =Γc,lin= 2 −α 3b2Σprp4Ω2p.. (24)
usingthetorqueformula(4)givenbyGoldreich & Tremaine m=1 „ «
X
(1979). This is applicable to the situation where the soft-
We remark that Ward (1992) obtained a corresponding ex-
eningparameter is significantly larger thanthescale height
pression expression for a vertically averaged potential, that
andweuseittoevaluatethecorotationtorqueinthatlimit.
To do this we adopt the approximation for Φ′ introduced has the same scaling when b is replaced by hp in the above
m
equation.
by Goldreich & Tremaine (1980) for m>0 in theform
Acomparisonofthetorquesgivenbyequations(5)and
Φ′m = π1 Z02πΦpcos(mϕ)dϕ∼−π2GrMppK0(mζ), (20) (c2a4n)tifnodribca>te∼s 1th.4ahtpt.hBeuetffnecotteofthsaotftienniandgdsithioonuldtobreeqsuiginriinfig-
b/h ≫ 1 (24) also formally requires b ≪ 1. This means, as
where K0 is the standard Bessel function and
weshallseebelow,thatthelimitwhere(24)appliesrequires
rathersmall h<∼0.025.
ζ = (r−rp)2 +b2. (21) We now directly compare linear calculations with the
s rp2 prediction of equation (24). The results are illustrated in
Fig. 5, for three different values of hp. Note that in the
Thisisvalidintheimportantdomainofinterestwherem∼
derivation of equation (24) it was assumed that b ≪ 1, so
1/b,and|rp−r|<∼brp.Summingthetorquesgivenby(4) we expect the torque to be given by equation (24) for h≪
overm, using (20) we obtain
b≪1.FromFig.5,weseethatforhp =0.025wecannicely
∞ 3 8q2 ∞ reproduce equation (24) for 0.05 < b < 0.2. However, as
Γc,m = 2 −α 3 Σprp4Ω2p mK0(mb)2, (22) hp increases there are increasing deviations from equation
mX=1 „ « mX=1 (24). For hp = 0.05 and 0.1 < b < 0.2, the maximum amd
where Σp denotes the unperturbed surface density at the minimumdeviationsarebyafactorof2and1.5respectively,
planets location. Noting that b is small and that the domi- while for hp = 0.1 the deviation always exceeds a factor of
nant contribution to the sum on the right hand side comes two.
8 S.-J. Paardekooper and J.C.B. Papaloizou
5 A COMPARISON OF THE EXPECTED
HORSESHOE DRAG WITH THE LINEAR
α=5/2
COROTATION TORQUE α=3/2
Inthissectionwecomparecorotationtorquesobtainedfrom 0 α=1/2
α=-1/2
linear calculations to the expected horseshoe drag. We ob-
tainthelatterfromequation(9)whichrequiresanestimate
ofxs whichweobtainfromequation(10)whichsimulations
have indicated gives a good estimate for b ∼ hp. We shall
obtain estimates for thehorseshoe drag directly from simu- q2 -500
Γ/
lations and make furthercomparisons below.
In Fig. 6, we show the linear corotation torque (solid
line), together with the horseshoe drag, obtained as indi-
catedabove,(dashedline)foradiscwithα=1/2,fordiffer-
entvaluesofthesmoothingparameterb.Forreasonableval-
-1000
uesofb,thehorseshoedragisstrongerthanthelinearcoro-
tation torquemakingthetorqueon thedisc more negative.
Thusthehorseshoedragontheplanetispositiveandlarger
than the linear corotation torque. For smoothing parame-
0 5 10 15 20 25 30
ters 0.02<b<0.03, this non-linear torqueis a factor 2−3 t (orbits)
larger than the linear torque. For smaller values of b, equa-
tion(10)predictsavalueofxs thatistoolargeandaccord- Figure 7.Totaltorqueontheplanetforanisothermal,inviscid
inglyahorseshoedragthatistoolarge.Inrealityweexpect hp=0.05discwithb=0.03andq=1.26·10−5,fordifferentsur-
both the horseshoe width and therefore the horseshoe drag facedensity profiles.Thedotted lines indicatethe linear torque,
forincreasingαfromtoptobottom.
to reach limiting values for small b and these values should
be larger than those we obtain here using values of b for
which equation (10) applies. A more complete discussion of
0.0015rp in both directions. Testshaveshown thatthisres-
theseaspectsisgiveninPaardekooper & Papaloizou(2009).
olution is sufficient to capture the horseshoe dynamics for
Also shown in Fig. 6 is the corotation torque obtained for
xs>0.004.Weincludealldiscmaterialinthetorquecalcu-
the corresponding three-dimensional disc (TTW02) which
lation (not excludingany material close to theplanet), and
has zero softening. For b = 0.02, we can match our 2D lin-
theplanet is introduced with its full mass at t=0. Forlow
earcalculationtothe3DresultofTTW02.However,forthe
mass planets, thisdoes not affect theresults.
same softening thehorseshoe drag is threetimes as large.
First, in Section 6.1, we study the development of the
For a large softening with b=0.7hp =0.035, thehorseshoe horseshoedraganditsrelationshiptolineartheory.Then,in
dragisequaltothe3Dunsoftenedlinearcorotation torque.
Section 6.2, we study the long-term behaviour (saturation)
However, numerical results show that for this smoothing
of thecorotation torque.
length, equation (10) actually slightly underestimates the
truevalueof xs, with theconsequencethat the intersection
of the horseshoe drag with the 3D linear corotation torque 6.1 Development of the non-linear torque
occursforapproximatelyb=0.8hp.Lackingafull3Dmodel
Westartbyconsideringaplanetwitharelativelyhighmass,
of the horseshoe region, it is difficult to say what value of
b to choose so that 2D results match 3D results. However, q = 1.26·10−5, which corresponds to 4 M⊕ around a So-
lar mass star. Although this is the planet with the largest
numerical results suggest that b<0.02 may be appropriate
massweconsider,ithasbeenexpectedtobewellwithinthe
(Masset et al. 2006). We come back to this issue in Section
linear regime (Masset et al. 2006). In Fig. 3, we show that
7, but it is good to keep in mind that the effects discussed
for a disc with constant specific vorticity,we can reproduce
inthenextsectionsmaywellbestrongerin3Dcalculations.
theexpectedlineartorqueforallreasonablevaluesofb.For
verysmallvaluesofb,onemayexpectadeparturefromlin-
eartheory,sinceatthispoint,anenvelopemayformthatis
gravitationally bound tothe planet,which probably should
6 HYDRODYNAMICAL SIMULATIONS
be excluded from the torque calculation. We do not con-
In the previous section, we have established that horseshoe siderthisregimehere,sinceweexpecttoadoptavalueofb
dragispotentiallymuchstrongerthanthelinearcorotation comparable to h in order to account for vertical averaging.
torque. We now turn to numerical hydrodynamic simula- Theimportant point is that we can match ourlinear calcu-
tionstoshowthatindeedstrongdeviationsfrom linearthe- lationstothehydrodynamicalsimulationsintheabsenceof
oryareencounteredinpractice.WeusetheRODEOmethod corotation torques(linear as well as non-linear).
(Paardekooper & Mellema2006b)intwospatialdimensions, This is furtherillustrated in Fig. 7, where we show the
onaregulargridextendingfromr=0.5rp tor=1.8rp and timeevolutionofthetotaltorqueontheplanetfordifferent
which coversthewhole 2π in azimuth.Sincewe want tore- surface density profiles. We see that the case with α=3/2
solve the horseshoe region for even the smallest planets we nicelyfallsonthecorrespondinglinearresult.Notealsothat
consider,arelativelyhighresolutionisrequired:1024cellsin this torque is set-up in approximately one orbital period.
theradialand4096 cellsintheazimuthaldirection,making This isto beexpected for both theLindblad andthelinear
the resolution at the location of the planet approximately corotation torque(see theAppendix).
Corotation torques and horseshoe drag 9
11000000
400 Γ +Γ
L C,lin
Γ +Γ
Γlin+ΓHS 550000 ΓLL+ΓHC+SΓHS
200
ν=0
ν=10-4
0 ν=10-3 00
q2 qq22
Γ/ ΓΓ//
-200 Γlin-ΓC,lin+ΓHS --550000
-400
--11000000
Γ
lin
-600
--11550000
0 5 10 15 20 25 30 --11..00 --00..55 00..00 00..55 αα 11..00 11..55 22..00 22..55
t (orbits)
Figure 9. Total torque on a q = 1.26·10−5 planet embedded
Fdiisgcuwreith8.αTo=tal−t1o/r2qufeoronbt=he0p.l0a3neatnidnqan=iso1t.2h6er·m1a0l−h5pf=or0d.0if5- in an isothermal hp = 0.05 inviscid disc for b = 0.03 after 30
orbits.Thesolidlineindicates thelineartorque,andthedashed
ferent magnitudes of the viscosity. The horizontal lines indicate,
line shows the total linear torque plus the estimated horseshoe
frombottom totop,thetotal lineartorque,theLindbladtorque
drag. The dotted line indicates the linear Lindblad torque plus
withthehorseshoetorque(foundusingtheprocedureoutlinedin
the non-linear horseshoe drag. Symbols denote results obtained
Section 5)added, andthetotal lineartorque withthe horseshoe
from hydrodynamical simulations for different values of the sur-
torqueadded.
facedensitypowerlawindexα.
Differentsurfacedensityprofilesgiveremarkablydiffer-
entresults.Allcasesexceptα=3/2showadeparturefrom large α viscosity parameter of αvisc = 0.4. This is because
linear theory, the sign of which is dictated by the vorten- onlythencantheviscositydirectlyaffectthehorseshoeturn.
sity gradient. This indicates that the corotation torque is In Section 6.2, we will see that lower viscosities are in fact
enhancedwithrespecttoitslinearvalue.Thisenhancement sufficient to keep thetorque unsaturated.
takes approximately 20 orbits to develop after which the Returning to inviscid discs, we show in Fig. 9 the to-
torques attain steady values for the remainder of the sim- tal torque on the disc for different values of α, confirming
ulations. This time period, as we show below, can be un- that the horseshoe drag indeed replaces the linear corota-
derstood as beingafraction of thelibration timescale, and tion torque. Note that we expect a torque reversal around
is therefore related to a non-linear effect. Note that in all α=−1,whilethelinearcalculationspredictthiswouldonly
cases, linear theory is only valid at early times (less than happenaroundα=−3.Notealsothatsteepdensityprofiles
about twoorbits) duringwhichasexpectedthelinearcoro- resultinanaccelerationofinwardmigrationwithrespectto
tation torque is set up on a dynamical time scale. At later thelinear estimate.
times, it gets replaced bythe non-linear horseshoe drag. Linear theory predicts that the torque should scale as
We take a closer look at the α=−1/2 case illustrated q2.Sinceinoursimplemodel,thehorseshoewidthscalesas
inFig.8.Thehorizontaldottedlinesindicate,from bottom q1/2, the non-linear contribution to the torque (the horse-
totop,thelineartorque,thelinearLindbladtorquewiththe shoedrag)alsoscalesasq2.Therefore, itcouldbemisinter-
horseshoe drag added, and the total linear torque with the preted as a linear effect. However, the time scale on which
horseshoe drag added. From the solid line, we see that the the horseshoe drag term is set up depends on the mass of
total torque can be understood as the linear torque, with the planet. We will see below that is takes a fixed fraction
thelinearcorotationtorquereplacedbythehorseshoedrag. of a libration time scale. This means that although the full
Wehavefound thistobetruefor all valuesof αconsidered non-linear torque scales as q2, this will not be the case at
(seeFig. 9).It is also expected from thediscussion given in intermediate stages. This is illustrated in Fig. 10, where we
Section 3. show thetime evolution of the torquefor different mass ra-
Themagnitudeofthenon-lineartorquedependsonthe tios. For linear torques, all curves would fall on top of each
magnitude of the viscosity in the disc. Recall that when- other. This is indeed true at early times, when the horse-
ever we include viscosity in the model, we do not allow for shoe drag has not yet developed. When the horseshoe drag
large scale mass flow with respect to the planet (see also takes over, however, planets of different mass give different
Section 7). It is easy to understand this dependence, since results; lower-mass planets take longer to develop the non-
thehorseshoedragmodelhingesonvortensityconservation linear torque.
whilethefluidexecutesahorseshoeturn.Forstrongenough Ifwerescalethetimeaxistothelibrationtime(seeFig.
viscosities, we can recover the linear torque, but note that 11),thecurvesfallontopofeachotheragain.Thisisbecause
the large values ν >∼ 10−3 required correspond to a very thetimescaletosetupthehorseshoedragisafixedfraction
10 S.-J. Paardekooper and J.C.B. Papaloizou
0
0
q=1.26e-5 ν=0
q=9.45e-6 ν=10-7
-100 q=6.34e-6 ν=10-5
q=3.15e-6
-500
-200
Γ/q2 -300 Γ/q2
-1000
-400
-1500
-500
-600 0 100 200 300 400 500
0 2 4 6 8 10 t (orbits)
t (orbits)
Figure 12.Total torqueonaplanetinanisothermalhp=0.05
Figure 10.Totaltorqueonaplanetinanisothermalhp=0.05 disc for q = 1.26·10−5, b = 0.02 and α = −1/2, for different
disc for b = 0.03 and α = −1/2, for different values of q. The valuesoftheviscosityν.
dottedlineindicatesthelineartorque.
0
200
-100
0
-200
q2
Γ/
q2 -300 b=0.03
Γ/ -200 b=0.02
q=1.26e-5
-400 q=9.45e-6 ν=0
q=6.34e-6 ν=10-5
q=3.15e-6 -400
-500
0 5 10 15 20 25 30
t (orbits)
-600
0.00 0.05 0.10 0.15
t (τlib) Figure 13.Total torqueonaplanetinanisothermalhp=0.05
discforq=1.26·10−5andα=−1,fordifferentvaluesofsoftening
Figure11.Totaltorqueontheplanetembeddedinanisothermal parameterbandviscosityν.
hp=0.05discforb=0.03andα=−1/2, fordifferentvalues of
q.Thetimeisinunitsofτlib,whichscalesasq−1/2. executed a horseshoe turn at constant entropy (replacing
specific vorticity applicable to thiscase).
ofthelibration time(approximately 15%, accordingtoFig.
6.2 Long-term evolution
11). One might expect that, as the scale of the region con-
tributing the torque is of order H, the time required would Lindbladtorquesgiverisetowavesthatcarryawayangular
beontheorderofafactorhpsmallerthanthelibrationtime. momentum.Therefore, theplanetcancontinuetoexchange
Thustheresultspresentedhereareconsistentwiththistime angular momentum with the disc at Lindblad resonances.
beingabout 2hpτlib.Thistypeofphenomenonwasalso dis- However,thereisnowavetransportinthecorotationregion,
cussed in Paardekooper & Papaloizou (2008), where it was andthereforeintheabsenceofanyotherformoftransport,
shown that thedevelopmentof thenon linear torqueis due only a finite amount of angular momentum exchange with
to an high density ridge resulting from material that has the planet can occur before the structure of the coorbital