The interaction of planets with a disc with MHD turbulence III: Flow morphology and conditions for gap formation in local and global simulations

1.1 MB·English

by  John C. B. Papaloizou

Mon.Not.R.Astron.Soc.000,000–000 (0000) Printed2February2008 (MNLATEXstylefilev1.4)
The interaction of planets with a disc with MHD
turbulence III: Flow morphology and conditions for gap
formation in local and global simulations
John C.B. Papaloizou, Richard P. Nelson & Mark D. Snellgrove
Astronomy Unit, Queen Mary, Universityof London, Mile End Rd, London E1 4NS
4
Received/Accepted
0
0
2
ABSTRACT
n
We present the results of both global cylindrical disc simulations and local shearing
a
box simulations of protoplanets interacting with a disc undergoing MHD turbulence
J
with zero net flux magnetic fields. We investigate the nature of the disc response and
0
conditions for gapformation. This issue is animportant one for determining the type
2
andnature ofthe migrationofthe protoplanet,withthe presenceof adeep gapbeing

believed to enable slower migration.
2
v For both types of simulation we find a common pattern of behaviour for which
1 the main parameter determining the nature of the response is MpR3/(M∗H3), with
5 Mp,M∗,R, and H being the protoplanet mass, the central mass, the orbital radius
3 and the disc semi–thickness respectively. We find that as MpR3/(M∗H3) is increased
8 to ∼ 0.1 the presence of the protoplanet is first indicated by the appearance of the
0
well known trailing wake which, although it may appear to be erratic on account of
03 the turbulence, appears to be well defined. Once MpR3/(M∗H3) exceeds a number
around unity a gap starts to develop inside which the magnetic energy density tends
/
h to be concentrated in the high density wakes. This condition for gap formation can
p beunderstoodfromsimpledimensionalconsiderationsoftheconditionsfornonlinear-
- ity and the balance of angular momentum transport due to Maxwell and Reynolds’
o
stresses with that due to tidal torques applied to the parameters of our simulations.
r
t An important result is that the basic flow morphology in the vicinity of the pro-
s
a toplanet is very similar in both the local and global simulations. This indicates that,
: regardless of potentially unwanted effects arising from the periodic boundary con-
v
ditions, local shearing box simulations, which are computationally less demanding,
i
X capturemuchofthephysicsofdisc–planetinteractions.Thusthey mayprovideause-
r fultoolforstudyingthelocalinteractionbetweenformingprotoplanetsandturbulent,
a protostellar discs.
Key words: accretion,accretiondisks—MHD,instabilities,turbulence–planetary
systems: formation, protoplanetary discs
1 INTRODUCTION planets are likely to form at significantly larger radii than
those observed, requiring a migration mechanism to bring
The recent and ongoing discovery of extrasolar giant plan- themclosertothecentralstar.Disc–protoplanetinteraction
ets has stimulated much current investigation of theories providesa natural mechanism.
of planet formation (e.g. Mayor & Queloz 1995; Marcy, Inthestandardpictureaprotoplanetexertstorqueson
Cochran, & Mayor 1999; Vogt et al. 2002). At the present a protostellar disc through the excitation of spiral density
time giant planets are believed to originate in protostellar waves(e.g.Goldreich &Tremaine1979). Thesewavescarry
discs.Possibleformation scenarios areeitherthroughdirect anassociated angularmomentumfluxwhichisdepositedin
gravitationalfragmentationofayoungprotostellardisc(e.g. thediscwherethewavesaredamped.Thisprocessresultsin
Boss2001)orthroughtheaccumulationofasolidcorewhich a negative torque acting on the protoplanet from the outer
thenundergoesmassbuildupthroughgasaccretiononcethe discandapositivetorqueactingonitfromthediscinterior
core mass reaches 15 Earth masses (e.g. Bodenheimer & toitsorbit,withtheouterdisctorquenormallydominating.
≃
Pollack 1986; Pollack et al. 1996). In either scenario giant Forlowmassplanetsthatinducelinearperturbationsinthe
(cid:13)c 0000RAS
2 J.C.B.Papaloizou, R.P. Nelson & M.D. Snellgrove
disc, the resulting migration is usually referred to as type I estimated but most other features of theglobal simulations
migration. are reproduced.
Previousworkontheinteractionbetweenaprotoplanet In this paper we focus on the morphology of the disc
andalaminardiscwithNavierStokesviscosity(Papaloizou response as a function of the protoplanet mass leaving the
& Lin 1984; Lin & Papaloizou 1993) indicates that a suf- analysis of the disc–protoplanet torques resulting from the
ficiently massive protoplanet can open up an annular gap interaction to a companion paper (Nelson & Papaloizou
in the disc centred on its orbital radius. For typical pro- 2003b – hereafter paper IV).
tostellar disc models the protoplanet needs to be approxi- The plan of thispaper is as follows:
mately a Jovian mass for gap formation to occur. Recent In 2wedescribeboththeinitialglobaldiscmodelsandthe
§
simulations (Bryden et al. 1999; Kley 1999; Lubow, Seib- shearing box models used for thesimulations.
ert, & Artymowicz 1999; D’Angelo, Henning, & Kley 2002) In 3 we discuss the nature of the waves that can propa-
§
examined the formation of gaps by giant protoplanets, and gate and how their excitation leads to the existence of the
also estimated the maximal gas accretion rate onto them. prominentwakethatisfrequentlyseeninthesimulationsof
The orbital evolution of a Jovian mass protoplanet embed- protoplanetdiscinteractionevenwhenthediscisturbulent.
ded in a standard laminar viscous protostellar disc model In 4 we describe thenumerical procedureused.
§
was studied by Nelson et al. (2000). They found that that In 5 we present the results of the local and global simula-
§
gap formation and accretion of the inner disc by the cen- tions. Finally, we summarize our results in 6 .
§
tral mass led to the formation of a low density inner cavity
in which the planet orbits. Interaction with the outer disc
resulted in inward migration on a time scale of a few 105
×
yr. Here the viscous evolution of the disc is responsible for 2 SET UP OF INITIAL MODELS
drivingtheinwardmigrationofthegiantplanet.Thismode
Inthispaperwehaveperformedbothglobalandlocalsimu-
of migration is normally referred to as type II migration.
lationsofprotoplanetsinteractingwithmodeldiscsinwhich
Recent work examining gap formation, mass accretion, and
MHD turbulence driven by the MRI operates. We consider
migrationtorquesinthreedimensionshasbeenpresentedby
simulationsforwhichthemagneticfieldhaszeronetfluxand
D’Angelo, Kley, & Henning (2003) and Bate et al. (2003).
thereforeaninternallygenerateddynamoismaintained.We
The disc models in these studies all adopted an anomalous
determinethediscresponseasafunctionoftheprotoplanet
discviscositymodelledthroughtheNavier–Stokesequations
massandasfarasispossiblestudyitscharacteristicsforthe
without consideration of its origin.
differentkindsof model focusing on theconditions that de-
Themost likely origin oftheviscosity isthrough MHD
terminewhethertheresponseisweak(linear)orstrong(non
turbulence resulting from the magnetorotational instability
linear).Wenowdescribethemodelsetupfortheglobaland
(MRI) (Balbus & Hawley 1991) and it has recently become
local simulations.
possible through improvements in computational resources
to simulate discs in which this underlying mechanism re-
sponsible forangular momentumtransport isexplicitly cal-
2.1 Global Simulations
culated.Thisisnecessarybecausetheturbulentfluctuations
do not necessarily result in transport phenomena that can The governing equations for MHD written in a frame ro-
bemodelled with theNavier Stokesequation. tating with arbitrary uniform angular velocity Ωpkˆ with kˆ
To this end Papaloizou & Nelson (2003) and Nelson & beingtheunit vectoralong therotation axisassumed tobe
Papaloizou (2003a) (hereafter papers I and II) developed in the vertical (z) direction are the continuityequation:
models of turbulentprotostellar accretion discs and consid-
∂ρ
ered the interaction with a giant protoplanet of 5 Jupiter + ρv=0, (1)
∂t ∇·
masses. The large mass was chosen to increase the scale of
the interaction so reducing the computational resources re- theequation of motion
quired. This protoplanet was massive enough to maintain
∂v p ( B) B
a deep gap separating the inner and outer disc and exert +v v+2Ωpkˆ v= ∇ Φ+ ∇× × ,(2)
∂t ·∇ × − ρ −∇ 4πρ
torques characteristic of type II migration. A study of gap
formationinturbulentdiscshasalsobeenpresentedbyWin- and theinduction equation
ters, Balbus, & Hawley (2003a)
∂B
= (v B). (3)
It is the purpose of this paper to extend our previous ∂t ∇× ×
work to a wider range of protoplanetary masses and disc
where v,P,ρ,B and Φ denote the fluid velocity, pressure,
parameters. In particular we study the nature of the disc
density and magnetic field respectively. The potential Φ =
protoplanet interaction as the protoplanet mass increases
and Tinovedsotigthatisewtheeuctoilnisdeitbioonthfogrlogbaapl afonrdmlaotciaoln.simulations. Φcernottr+ifuΦgGalcponotteanintsiacloΦnrtoritb=ut−io(n1s/d2u)Ωe2pt|okˆg×rarv|i2t.y,ΦG andthe
The gravitational contribution is taken to be due to a
The global simulations are more realistic but too computa-
centralmassM∗ andplanetwithmassMp.Thusincylindri-
tionally demanding for significant extended parameter sur-
calcoordinates(r,φ,z)withtheplanetlocatedat(rp,φp,0)
veys.Ontheotherhandlocalshearingboxsimulationsbeing
and thestar located at (r∗,φ∗,0) ΦG =Φc+Φp, where
less computationally demanding enable a wider parameter
sduergvreeyeoaftshyimghmeertrryestohluetyiodno.nOoftcaolluorwsembiegcraautisoenorfattehsethoigbhe Φc =−|rG−Mr∗∗| and Φp =− r2+rp2−2rGrpMcops(φ−φp)+b2.(4)
p
(cid:13)c 0000RAS,MNRAS000,000–000
Protoplanets and MHD turbulence 3
Here,asinpapersIandII,forreasonsofcomputationalex- G5 it was placed at (rp r∗ ,φp) = (2.5,π/4). All simu-
| − |
pediency, we have neglected the dependence of the gravita- lations were performed in a rotating frame whose angular
tionalpotentialsonzalongwiththeverticalstratificationof velocity equalled that of the circular Keplerian planetary
thedisc.Thustheglobalsimulationsareofcylindricaldiscs orbit.The unitof time adopted when discussing theresults
(e.g.Armitage1998,2001;Hawley2000,2001;Steinacker& is Ω−1 evaluated at r = 1, the inner boundary of the com-
Papaloizou 2002). To model the effects of the reduction of putational domain. We note that the orbital period at this
the planet potential with vertical height, we have incorpo- radius is P(r = 1) = 2π, the orbital period of a planet at
rated a softening length b in its potential. rp = 2.5 is P(r = 2.5) 24, and the orbital period of a
≃
Weuse a locally isothermal equation of state in theform planet at rp =3 is P(r=3) 32.
≃
In order to explore the disc–planet interaction over a
P =ρ c(r)2, (5)
· wider range of parameter space and with higher resolution
withcdenotingthesoundspeedwhichisspecifiedasafixed inthevicinityoftheplanet,wehaveadoptedalocalshearing
functionofthecylindricalradiusr.Theabovegivesthebasic box approximation. This we describe below.
equations used for theglobal simulations.
2.2 Local Shearing Box Simulations
2.1.1 Initial and Boundary Conditions
In the shearing box limit (Goldreich & Lynden-Bell 1965)
The global simulations presented here all use the same un- we consider a small Cartesian box centred on the point at
derlying turbulent disc model. This model has a constant which the Keplerian angular velocity is Ωp. Wedefine local
aspect ratio H/r c(r)/(rΩ) = 0.07, where Ω is the disc Cartesian coordinates (x,y,z) with associated unit vectors
≡
angularvelocityasmeasuredintheinertialframe.Theinner (i,j,k).Thedirectioniistakentobeoutwardalongtheline
radial boundaryof thecomputational domain is at Rin =1 joiningtheorigin tothecentralobjectwhilekpointsinthe
andtheouterboundaryisatRout=8.Theboundarycondi- vertical direction. The direction j is along the unperturbed
tions employed are very similar to those described in paper rotational shear flow.
I. Regions of the disc in the close vicinity of the inner and The equation of motion may then bewritten as
outerboundariesweregivennon–Keplerianangularvelocity
∂v p ( B) B
pvarolufielsesofththaet daerenssittyabilneotrodetrhetoMmRaIin,taanindrwadhiiachlhhyadvreosltaartgiec ∂t+v·∇v+2Ωpkˆ×v−3Ω2pxi=−∇ρ −∇Φp+ ∇×4πρ× .(6)
equilibrium. These regions act as buffer zones that prevent Theinductionandcontinuityequationsretainthesame
the penetration of magnetic field to the radial boundaries, formasequations(3)and(1)butareofcourseexpressedin
thusmaintainingzero-netfluxinthecomputationaldomain. component form in the local box centred Cartesian coordi-
The inner buffer zone runs from r = 1.14 to r = Rin, and nates.Werecallthattheterm xinequation(6)isderived
theouter bufferzone runsfrom r=7.0 to r=Rout. fromafirstorderTaylorexpan∝sionaboutx=0ofthecom-
Thediscwasinitiatedwithazeronetfluxtoroidalmag- bination of thegravitational acceleration duetothecentral
netic field in a finite annulus in the disc between the radii massandthecentrifugalacceleration (Goldreich&Lynden-
2 r 6.Theinitialfieldvariedsinusoidallyoverfourcom- Bell 1965). However, as for the global disc simulations, we
≤ ≤
plete cycles within this radial domain, and the initial ratio haveneglectedthez dependenceofthegravitationalpoten-
ofthevolumeintegratedmagnetic energytothevolumein- tialandthereforeverticalstratification.Thusinthisrespect
tegrated pressurewas 1/β =0.032. The initial disc model thesimulations areof unstratifiedboxesofthetypeconsid-
h i
had azimuthal domain running between φ = 0 to π/4, and ered by Hawley, Gammie & Balbus (1995).
vertical domain from zmin = 0.21 to zmax = 0.21. Peri- Here,whenintroduced,theplanetislocatedinthecen-
−
odicboundaryconditionswereemployedintheverticaland treof thebox, thus
azimuthal directions.
TheunperturbeddiscmodelwasevolveduntiltheMRI Φp = GMp . (7)
produced a turbulent disc model with a statistically steady − x2+y2+b2
turbulence.This model had a volumeaveraged valueof the p
Inthesteadystateboxwithnoplanetormagneticfield,
ratio of magnetic energy to mean pressure of 0.011 and
≃ theequilibrium velocity is dueto Keplerian shear and thus
avolumeaveraged valueof theShakura–Sunyaevstress pa-
rameter α 7 10−3. In order to generate the full 2π disc v=(0, 3Ωpx/2,0). (8)
≃ × −
modelsusedtostudytheinteractionbetweenturbulentdiscs
We adopt an isothermal equation of state P = ρc2,
and protoplanets, eight identical copies of this relaxed disc
with c being thesound speed which is taken to beconstant
sector werejoined together. Themodel described in table2
throughout thebox.
with azimuthal domain π/2 (run G5) contained only two
copies of the π/4 sector, and periodic boundary conditions
were employed in theazimuthal domain.
2.2.1 Initial and Boundary Conditions
Gravitational softening of the protoplanet was em-
ployed in all simulations. In runs G1–G4 described in ta- Theshearingboxispresumedtorepresentalocalpatchofa
ble 2, the softening parameter b=0.33H. For run G5, this differentially rotating disc. Thus the appropriate boundary
was reduced to b = 0.1H. For runs G1, G2, and G3, the conditions on the bounding faces y = constant = Y and
±
planet was placed at a radius such that rp r∗ = 3 and z =constant = Z derive from the requirement of period-
| − | ±
at azimuthal location φp = π in the disc. For runs G4 it icity in the local Cartesian coordinate directions normal to
was placed such that (rp r∗ ,φp) = (2.5,π) and for run theboundaries.Ontheboundaryfacesx=constant= X,
| − | ±
(cid:13)c 0000RAS,MNRAS000,000–000
4 J.C.B.Papaloizou, R.P. Nelson & M.D. Snellgrove
the boundary requirement is for periodicity in local shear- nitude of the disc density may be arbitrarily scaled. The
ing coordinates. Thus for any state variable F(x,y,z), the models simulated here all start with a uniform density ρ
0
condition is that F(X,y,z) = F( X,y 3ΩpXt,z). This which, together with H, may be used to specify the mass
− −
meansthatinformationononeradialboundaryfaceiscom- scaling. Then we define a dimensionless density and mag-
municated to the other boundary face at a location in the neticfieldthroughρ′=ρ0ρandB′ =B/(ΩpH√4πρ0).The
azimuthalcoordinateyshiftedbythedistancethefaceshave equationofmotion maybewrittenin theabovedimension-
sheared apartsincethestart ofthesimulation,t,(seeHaw- less variables in the form
ley,Gammie & Balbus 1995).
paraFmoertearllwmasodtealksenadaosptbin=g 0a.3sHhe.aTrinhge pbootxe,ntthiaelswoaftsenfliantg- ∂∂vt′′ +v′·∇v′+2kˆ×v′−3x′i=
tened(madetoattainaconstantvalueinacontinuousman-
ner)atlargedistancesfromtheprotoplanetinordertosat- (ρ′c′2) ′ ( B′) B′
∇ Φ + ∇× × , (9)
isfy the periodic boundary conditions. In simulations Ba0– − ρ′ −∇ p ρ′
Ba4(describedintable1)thiswasdoneoutsidethecylinder
where
ofradius3H centredontheplanet.InsimulationsBb1–Bb4
′
this was done only close to theboundary. Φ′ = Mp , (10)
ThesimulationsBa0–Ba4occupiedatotalwidth8H in p − x′2+y′2+b′2
xand4πH iny.Theseboxeswerefoundtobelargeenough p ′
with dimensionless softening parameter b = b/H and of
toaccomodatethescaleofthediscresponsetoforcingbythe
coursethespatialderivativesarewithrespecttothedimen-
protoplanet and enable gap formation. But notethat inter-
sionless coordinates. Similarly the dimensionless induction
ferencefrom protoplanetsinneighbouringboxesimplied by
equation becomes
theperiodic boundaryconditions is not eliminated entirely.
In the small protoplanet mass regime waves propagate be- ∂B′ ′ ′
= (v B). (11)
tweenboxesandinthenonlinearregimetorquesexertedby ∂t′ ∇× ×
protoplanets in neighbouring boxes may influence the gap
and thedimensionless continuity equation is
formation process. However, comparison between local and
global runs indicates that such effects are not major and ∂ρ′ ′ ′
+ (ρv)=0. (12)
do not prevent the box runs from capturing the essential ∂t′ ∇ ·
physics. We remark that in dimensionless coordinates the
To examine the effect of varying box size we ran simu- boundaries of the box are x′ = X/H, y′ = Y/H,
lationsBb1–Bb4whichwereofthesameresolution buthad z′ = Z/H. Also if the magnetic fi±eld has zero n±et flux
twice the extentin y. and is±dynamo generated and thus a spontaneous prod-
Allofthesemodelswereinitiatedbyimposingavertical uct of the simulation, given that the only parameter oc-
field,withzeronetfluxthroughthebox,varyingsinusoidally curring in equations (9 - 12 ) is the dimensionless mass
in x with a wavelength of H. The amplitude was such that Mp′ = MpR3/(M∗H3), we should be able to consider the
theinitial valueof theratio of thetotal magneticenergy to dependenceofthetimeaveragedoutcomeofasimulationto
volume integrated pressure was 0.0025. The initial velocity beonly on X/H,Y/H,Z/H, and Mp′ =MpR3/(M∗H3).
in the x direction at each grid point was chosen to be the For a planet embedded in a box with natural verti-
productofarandomnumberbetween 1.0and1.0and0.1c. cal scale H with sufficiently large X/H to incorporate the
−
Theinitialvelocityinthey–directionisgivenbyequation8, tidalinteraction without interferencefrom other boxescen-
andthez componentofvelocitywasinitializedtozero.The tred on different radii, the only parameters determining
unitoftimeforallsimulationswastakentobeΩ−p1.Inthese long term averaged properties of a simulation should be
unitsthe orbital period at the centreof the box is 2π. Y/H and M′. Thus for such a quantity Q, we should have
p
Q=f(Y/H,MpR3/(M∗H3)).
Note that the above discussion implies that for a box
2.3 Dimensionless Variables and Scaling with
of fixed dimension, the only distinguishing parameter is
Planet Mass
Mp′ = MpR3/(M∗H3). This parameter may also be inter-
It is helpful to write the governing equations for the shear- preted as the cube of the ratio of the Hill radius to disc
ing box in dimensionless form. In this way the dependence scale height and the condition that M′ > 1 leads to the
p
of theoutcomes of the simulations on planet mass and disc thermal condition of Lin & Papaloizou (19∼93) that the Hill
parameters is made clearer (see also Korycansky & Pa- radiusshouldexceedthediscscaleheightforgapformation
paloizou 1996). To do this we adopt dimensionless coordi- to occur. From Korycansky & Papaloizou (1996) this con-
nates (x′,y′,z′)=(x/H,y/H,z/H), where the length scale ditionisalsorequiredinorderthattheperturbationdueto
H = c/Ωp would correspond to the disc scale height were theprotoplanet be non linear.
it vertically stratified and the dimensionless time t′ = Ωpt. Theabovediscussionindicatesthatwerewetovarythe
We also note that the Keplerian relation Ω2p = GM∗/R3, boxsizeintheazimuthal(y)direction,theconditionforgap
′
may be used to relate the rotation rate of the centre of the formation should be of the form Mp > fg(Y/H). As gap
box to a putative central mass M∗ and orbital radius R. formation by the same planet in a longer box is expected
We now introduce the dimensionless velocity, sound speed to be more difficult, we expect fg to be a non decreasing
andplanetmassthroughv′=v/(ΩpH),c′ =c/(ΩpH),and function of its argument. In order to correspond to a full
Mp′ =MpR3/(M∗H3)respectively.Aswedonotincludethe circle of radius R,we requirethat Y =πR.
self-gravity of the disc in calculating its response, the mag- Itisperhapssuggestivetopointoutthatiffg(πR/H)=
(cid:13)c 0000RAS,MNRAS000,000–000
Protoplanets and MHD turbulence 5
Model z domain xdomain y domain t1 t2 (MpR3/(M∗H3) nz nx ny
Ba0 (−H/2,H/2) (−4H,4H) (−2πH,2πH) 0.0 637 0.0 35 261 200
Ba1 (−H/2,H/2) (−4H,4H) (−2πH,2πH) 354 621 0.1 35 261 200
Ba2 (−H/2,H/2) (−4H,4H) (−2πH,2πH) 354 644 0.3 35 261 200
Ba3 (−H/2,H/2) (−4H,4H) (−2πH,2πH) 354 619 1.0 35 261 200
Ba4 (−H/2,H/2) (−4H,4H) (−2πH,2πH) 354 650 2.0 35 261 200
Bb1 (−H/2,H/2) (−4H,4H) (−4πH,4πH) 216 392 0.2 35 261 396
Bb2 (−H/2,H/2) (−4H,4H) (−4πH,4πH) 216 431 0.5 35 261 396
Bb3 (−H/2,H/2) (−4H,4H) (−4πH,4πH) 216 411 1.0 35 261 396
Bb4 (−H/2,H/2) (−4H,4H) (−4πH,4πH) 216 485 2.0 35 261 396
Table 1.Parametersoftheshearingboxsimulations:Thefirstcolumngivesthesimulationlabel,thesecond,thirdandfourthgivethe
extent ofthe coordinate domainsconsidered. Thefirstninesimulationswereperformedinshearingboxes withthe x,y, andz domains
referringtotheCartesiandomains.Thelasttwosimulationswereglobalsothatinthesecasesthex,y,andzdomainsrefertother,φ,and
zdomainsforthecylindricalcoordinatesystemused.Thefifthandsixthcolumngivethestartandendtimesifthesimulationmeasured
indimensionlessunits.Simulationswithprotoplanetswerestartedbyinsertingtheprotoplanetintotheturbulentmodelgeneratedfrom
thesimulationBa0attimeindicated. Theseventh columngivesthevalueof(GMp/(Ω2H3).Theeighth, ninthandtenth columns give
the number of computational grid points in the x,y, and z coordinates respectively. These numbers include any ghost zones used to
handleboundaryconditions.
40α(R/H),theconditionforgapformation indicatedabove ally becoming a spiral wave propagating away from it (e.g.
becomes Mp/M∗ >40αH2/(R2) which is of the same form Kley1999;Nelsonetal2000;Lubow,Seibert&Artymowicz
as the viscous criterion of Lin & Papaloizou (1993) which 1999; Ogilvie & Lubow 2002). Such structures, though dis-
wasobtainedbyrequiringangularmomentumtransportdue torted, are also seen in our simulations of discs with MHD
to excited waves to exceed that due to viscosity. However, turbulence in which there are strong and highly time vari-
it is important to note that if we introduce such an α pa- able fluctuations, even for weakly perturbing protoplanets.
rameter here, it cannot correspond to an imposed viscosity In view of their importance and persistence we here give a
asinLin&Papaloizou(1993).Insteaditmustberelatedto briefdiscussion oftheirinterpretationasbeingdefinedbya
the turbulent transport arising from MHD turbulence. We characteristicrayemanatingfromtheprotoplanetthatsep-
take it to correspond to a time average of the volume av- aratestheflowshearingpastitintopartscausallyconnected
eraged stress BxBy/(4π) expressed inunitsofthevolume andnotcausally connectedbywaves.Thewakesthusnatu-
h i
averagedpressure P ,thuscorrespondingtotheShakura& rallycorrespondtothelocationastrongwavesentoutfrom
h i
Sunyaev(1973) parameterization.Thisprocedureisreason- theprotoplanet into theflow shearing past it.
ably well defined in the parts of the disc not strongly per-
turbedbyaprotoplanet(eg.PaperII).Forsimulationswith
zero net magnetic flux α 5 10−3 but it takes on larger 3.1 The Characteristic Rays
∼ ×
values when a net poloidal or toroidal flux is imposed (eg.
The set of basic equations defined by (1), (2) and (3) are
Hawley, Gammie & Balbus 1996; Steinacker & Papaloizou
well known to form a hyperbolic set which can be written
2002).
Adopting α = 5 10−3, suggests fg(Y/H) in the standard form (Courant & Hilbert 1962)
× ∼
0.07(Y/H). We recall that should 0.07(Y/H) be less than
unity,therequirementofnonlinearperturbationbythepro- A(k)∂Ui +B(k) =0 (k=1,2,....,n). (14)
toplanet in order for gap formation to occur is Mp/M∗ > Xi,j (cid:18) ij ∂xj(cid:19)
Ct(H/R)3, whereCt is a constant of order unity(Korycan- Here the set of equations (1), (2) and (3) are arranged into
sky&Papaloizou 1996). Boththeseconditionscanbecom-
a set of n simultaneous first order PDEs for the n compo-
bined to give therequirement that
nents Ui of U which are in turn composed of ρ together
MpR3/(M∗H3)>max(Ct,0.07(Y/H)) (13) withthecomponentsofBandv.Thexj (j =1−4)arethe
three spatial coordinates and the time t respectively. The
In this way we may see how both the thermal and viscous
characteristic rays, applicable to both the global and local
conditionsofLin&Papaloizou(1993)mayberelatedtothe
formulation may be found by first obtaining the local dis-
simulations.
persion relation associated with (14). To do this one sets
Ui = Ui0exp(ikjxj iωt) under the assumption that the
−
magnitudes of the components of the wave vector k(l =
l
3 WAKES AND CHARACTERISTICS 1,2,3) ≡ k = (kx,ky,kz) and the wave frequency |ω| are
verylarge.Thedispersion relationisthengivenbythevan-
A well known phenomenon found in simulations of proto- ishing of a determinant of coefficients |A(ijk)kj−A(i4k)ω|=0.
planets interacting with laminar discs is a prominent wake As is very well known (eg. Campbell 1997) this dispersion
that is observed to trail away from the protoplanet eventu- relationleadstothreedistinctwavemodes.Thefirstarethe
(cid:13)c 0000RAS,MNRAS000,000–000
6 J.C.B.Papaloizou, R.P. Nelson & M.D. Snellgrove
Model φdomain H/r Mp/M∗ (MpR3/(M∗H3) nr nφ nz
G1 2π 0.07 1×10−5 0.03 450 1092 40
G2 2π 0.07 3×10−5 0.09 450 1092 40
G3 2π 0.07 1×10−4 0.30 450 1092 40
G4 2π 0.07 3×10−3 8.75 450 1092 40
G5 π/2 0.07 3×10−3 8.75 450 276 40
Table2.Parametersoftheglobalsimulations:Thefirstcolumngivesthesimulationlabel,thesecondgivestheextentoftheazimuthal
domain,thethirdgivestheH/rvalueofthedisc,andthefourthgivestheprotoplanet-starmasratio.ThefifthgivestheratioofMp/M∗
to(H/r)3.Thesixth,seventh,andeighthcolumnsdescribethenumberofgridcellsusedineachcoordinate direction.
Alfv´enwavespropagatinginoppositedirectionswhichhave theratioofthecomponentsandnotthemagnitudematters.
ω = k v k vA, with vA = vAB/B, where vA is the The natural choice is kz = 0, and ky/kx such that ω =
Alfv´e−n sp·eed±be·ing such that v2 = B|2/|(4πρ). The other 0. This is because we expect a zero frequency wave to be
A | |
modes are fast and slow magnetosonic waves which satisfy excited by the planet in a frame in which that is at rest.
However,notethatin aturbulenttimedependentsituation
(ω+kv)2= k2(c2+vA2) k4(c4+vA4)−2k2c2(k·vA)2).(15) ω would not be constant.
· 2 ±p 2 The simplest situation is one when the flow is time in-
Herek= k andwenotethatthefrequency,ωoccursadded dependent and the velocity is given by (8). In that case it
| |
to the quantity k v in these relations. This is due to the follows from equations (16) that both ω = 0 and ky are
·
Dopplershiftoccurringonaccountofthemotionofthefluid. constant.
For each of these waves we may specify ω = ω(k,r,t). It further follows from equations (16) that the charac-
Then the characteristic rays follow from integrating the teristic rays then satisfy
equations of geometric optics (e.g. Whitham 1999)
dy
dr ∂ω dk ∂ω = v2/c2 1. (18)
= , and = . (16) dx − y −
dt ∂k dt −∂r p
To be specific, we here concentrate on the fast mag- Usingequation(8)for vy andintegratingweobtainfor
netosonic mode which in our case leads to the most rapid thecharacteristic ray for x>x
0
communication. Further,as for the most part the magnetic
energyison theorderofonepercentofthethermalenergy, x 1
y/x = x2/x2 1+ ln x2/x2 1+x/x ,(19)
we shall neglect the field so that the mode becomes an or- 0 −2x0p 0− 2 (cid:16)p 0− 0(cid:17)
dinary sound wave such that
with x = 2H/3. To complete the picture we remark that
ω= k v+ck. (17) 0
− · thecharacteristicforx<0,canbefoundfromthecondition
We consider the general situation to be one for which y( x)= y(x).
− −
the flow streams past the protoplanet from y > 0 to y < 0 We comment that the computational domains in all of
for x > 0 and with the flow moving in the opposite sense our simulations are extensive enough to contain the char-
for x < 0. Considering without loss of generality the case acteristic associated with the waves emitted by an exciting
y > 0, when the flow streams past supersonically, if we ne- object, be it turbulent fluctuations or a protoplanet (i.e. x
glect dependence on z, it is divided into two parts by the exceeds2/3H bya significant margin).
characteristic originating close to theplanet from the point
with y = 0 and the smallest possible positive value of x at
whichtheflowspeedisequaltothewavespeed.Iftheplanet
had no gravity,the fluid that is upstream of this character-
istic would be unaware of thepresence of theplanet.
Thefluidthatfirstknowsabouttheplanet,namelythat
4 NUMERICAL PROCEDURE
onthecharacteristic isanaturalchannelforwavestoprop-
agatethatareexcitedbytheplanetwhichiswhythisregion The numerical scheme that we employ is based on a spa-
issoprominentinthesimulations.Notethattheconceptof tially second–order accurate method that computesthe ad-
characteristicraysisrobustinthatitcanapplytoageneral vection usingthemonotonic transport algorithm (VanLeer
time dependent fluid with turbulence. Also with z depen- 1977). The MHD section of the code uses the method of
dence restored we would expect the fluid to be divided by characteristicsconstrained transport(MOCCT) asoutlined
a surface made up from many rays rather than a single ray in Hawley & Stone (1995) and implemented in the ZEUS
whichisofthesameformindependentofz.Thiswouldlead code. The code has been developed from a version of NIR-
to some blurring of the wake. VANA originally written by U. Ziegler (see SP, Ziegler &
Acharacteristicrayisdefinedbyitspointoforigin(here Ru¨diger (2000), and references therein)
specifiedabove)andthevalueofk.However,notethatonly
(cid:13)c 0000RAS,MNRAS000,000–000
Protoplanets and MHD turbulence 7
Figure 1.Theratioof thetotal magnetic energyto thevolume Figure 3. Density contour plot in a (x,y) plane for simulation
integrated pressure, h1/βi, is plotted as a function of time for Ba0whichhasnoembeddedprotoplanet.Theextensionsare8H
modelsBa0,Ba1,Ba2,Ba3,andBa4.Theplotsbegintoseparate in the x direction and 4πH in the y direction. The plot applies
after the protoplanet is inserted at time 354.. At time 500 the neartheendofthesimulation.Therangeofvariationintheden-
curves correspond to simulations Ba4, Ba1, Ba2,Ba3 and Ba0 sityisaboutafactorofthree.Thisistypicalofwhatisobtained
with decreasing values of h1/βi respectively. Although there are inany (x,y) planeonce astatisticallysteady state turbulence is
large fluctuations, the maximum variation at any time is factor setup.Notethetrailingelongated structures.
of 3.This issignificantlysmallerthan thedensity depressionin
thedeepest gap(seefigure10).
4.1 Vertically and Horizontally Averaged Stresses
and Angular Momentum Transport
Inordertodescribeaveragepropertiesoftheturbulentmod-
els,following papersIandIIweusequantitiesthatarever-
tically and azimuthally averaged over the (y,z) domain for
the shearing boxes or the (φ,z) domain for the global sim-
ulations (e.g. Hawley 2000). Thus we may consider y φ
≡
and x r in these cases. Sometimes an additional time av-
≡
erage may be adopted. The vertical and azimuthal average
of some quantityQ is then defined through
ρQdzdy
Q= . (20)
R ρdzdy
R
The disc surface density is given by
1
Σ= ρdzdy. (21)
2Y Z
The vertically and azimuthally averaged Maxwell and
Reynoldsstresses, are given by:
Figure2.Thevolumeaveragedstressparameterhαiisplottedas
afunctionoftimeformodel Ba0whichhas noprotoplanet. The TM =2YΣ BxBy (22)
noisy full curve contains contributions from both the magnetic (cid:18) 4πρ (cid:19)
andReynolds’stresseswhilethedottedcurvegivesthemagnetic
and
contribution. The dashed curve gives the running time average
of the total stress. The large initial values are due to the initial TRe=2YΣδvxδvy. (23)
strong instability which produces a channel phase. However, af-
ter about four orbits at the centre of the box, the running time Here the velocity fluctuations δvx and δvy are defined
averageiswithinafactoroftwoofthefinalvalue. through,
δvx =vx vx, (24)
−
δvy =vy vy. (25)
−
ThehorizontallyandverticallyaveragedShakura&Sunyaev
(1973) α stress parameter appropriate to the total stress is
given by
(cid:13)c 0000RAS,MNRAS000,000–000
8 J.C.B.Papaloizou, R.P. Nelson & M.D. Snellgrove
Figure 5. The volume averaged total stress parameter hαi is
plottedasafunctionoftimeformodelsBa1,Ba2,Ba3,andBa4.
Theplotsbeginfromthemomenttheprotoplanetisinserted.The
noisyuppercurvescontaincontributionsfromboththemagnetic
andReynolds’stresses.Theycanbeidentifiedatearlytimesasfor
Ba4(uppermost fullcurve), Ba3(uppermost dotted curve), Ba2
(next highest full curve), Ba1 (next highest dotted curve). The
lowerandmuchlessnoisycurvesgivethemagneticcontribution.
These aresimilaruntil they separate at about time 500. At this
time they may be identified with increasing values of the stress
parameterasassociatedwithBa1,Ba3,Ba2andBa4respectively.
Thevaluesofthetotalhαishowasharpincreasewhentheproto-
planet is inserted in the cases with higher mass protoplanets on
accountoftheexcitationofdensitywaveswhichprovideangular
momentumtransport.
5 NUMERICAL RESULTS
We now present numerical results for the evolution of the
Figure 4. Contours of magnetic energy density corresponding simulations.
to the plot in figure 3. The upper plot corresponds to the (x,y)
plane illustrated in figure 3 while the lower plot is taken in a
characteristic(x,z)plane.TheextensionsareHinthezdirection 5.1 Shearing Box Simulations
and8H inthexdirection.Duetoregionsofverysmallmagnetic
All of the shearing boxes are derived from the simulation
field, the range of variation of magnetic energy density is much
larger than that of the density being as much as ∼ 107. These Ba0whichwascarriedoutwithoutaprotoplanet.Itsettled
plots are typical of what is found for the turbulence once it has into a quasi steady state with conserved zero net flux and
reached a statistically steady state. Note that trailingelongated dynamo maintained turbulence as has been been described
features are also apparent in the distribution of the magnetic byothers(e.g.Hawley,Gammie&Balbus1996).Thismodel
energydensity. was run for 100 orbits at the box centre. The ratio of the
total magnetic energy to the volume integrated pressure,
B2/(8π)dV/ PdV = 1/β , is plotted as a function of
h i
Rtime in figure 1R.
After running for 354 time units, this model was used
toprovideinitialconditionsformodelswithprotoplanetsof
α= TRe−TM , (26) varying masses. Simulations Ba1, Ba2, Ba3, and Ba4 were
2YΣ(P/ρ) begun in this way with MpR3/(M∗H3) = 0.1,0.3,1.0 and
2.0respectively.Wealsoplot B2/(8π)dV/ PdV = 1/β
The angular momentum flow across a line of constant x is as a function of time in figureR1 for modelsRBa1, Ba2,hBa3i,
given by andBa4.Theplotsbegintoseparatefromeachothershortly
after after the simulations with protoplanets have begun.
=ΣWαP/ρ. (27) Note that although 1/β is in general confined between
F h i
0.008 and 0.02 with or without a protoplanet, there are
Here W = R2 for the shearing box simulations or W = r2 strong fluctuations on a range of timescales down to the
for the global simulations. orbitaltimescale.Theexistenceofsuchfluctuationsandthe
(cid:13)c 0000RAS,MNRAS000,000–000
Protoplanets and MHD turbulence 9
Figure 6. Density contour plots in a typical (x,y) plane for Figure 7. As in figure 6 but in this case magnetic energy den-
simulations Ba1, top left, Ba2 top right, Ba3 bottom left, and sitycontoursareplotted.AsMpR3/(M∗H3)increasesandagap
Ba4 bottom right taken near the end of the simulations. These forms,themagneticenergygenerallydecreasesinthegapregion.
had MpR3/(M∗H3) = 0.1,0.3,1.0 and 2.0 respectively. As However,someremainsconcentratedinregionsnearthehighden-
MpR3/(M∗H3) increases, the wake becomes more prominent, sitywakes.Thequalitativestructureseeninthesesimulationsis
materialisaccretedbytheprotoplanetandispushedtowardthe maintainedonceaquasisteadystatehasbeenattained.
radialboundariesasagapisformed.Thelocationofthecharac-
teristicraygivenbyequation(19)isalsoplotted.Thequalitative
structure seen in these simulations is maintained once a quasi initial strong instability which produced a channel phase.
steadystatehasbeenattained. However, after about four orbits at the centre of the box,
therunningtimeaverageiswithinafactoroftwoofitsfinal
valueof0.008.Thisvalueisinagreementwiththatgivenin
difficulties they cause for defining mean state variables for Winters, Balbus & Hawley (2003b). But note that therun-
the turbulent flow has been emphasised recently by several ning time average still shows variations at the 10 percent
authors (e.g. papers I and II; Winters, Balbus & Hawley level after 100 orbits.
2003b; Armitage 2002). The pattern of turbulence in simulation Ba0 is estab-
Thevolumeaveragedstressparameter α isplottedas lished after several orbits. A typical density contourplot in
h i
a function of time for model Ba0 in figure 2. As for 1/β , a(x,y)planeispresentedinfigure3neartheendofthesim-
h i
erratic fluctuations are seen. The sum of the magnetic and ulation. The range of variation in the density in this plane
Reynolds’ stresses nearly always exceeds the magnetic con- is about a factor of three. The situation illustrated here is
tribution for which the short term fluctuations are signifi- typicalofwhat isobtainedinany(x,y)planeonceastatis-
cantly less severe. In these respects the behaviour of local tically steady state turbulence is set up. Elongated trailing
andglobalsimulationsaresimilar(Steinacker&Papaloizou structuresare noticeable. Weplot contours of themagnetic
2002,papersIandII,andsection5.2).Therunningtimeav- energy density in figure 4 which correspond to the plot in
erageof α isalsoplottedandthisshowsmuchsmootherbe- figure 3. The upper plot corresponds to the (x,y) plane il-
h i
haviour. Somewhat large initial values occurbecause of the lustrated in figure 3 while the lower plot is taken in a char-
(cid:13)c 0000RAS,MNRAS000,000–000
10 J.C.B.Papaloizou, R.P. Nelson & M.D. Snellgrove
Figure8.Plotscorrespondingtothoseinfigure6butformodels Figure9.Plotscorrespondingtothoseinfigure7butformodels
Bb1,Bb2,Bb3,Bb4.Asimilarsequence isindicated. Bb1,Bb2,Bb3,Bb4.Asimilarsequence isindicated.
acteristic(x,z)plane.Duetotheexistenceofregionswhere withH/R=0.1presentedinpaperIIandintheglobalsim-
the magnetic field is near zero, the range of the magnetic ulations described in section 5.2. After a sharp initial rise
energydensityismuchlargerthanthatofthemassdensity, thevolumeaveragedReynolds’stressheretendstodecrease
being as much as 107. The plots are typical of what is as a gap is opened and material is moved away from the
∼
foundoncetheturbulencehasreachedaquasisteadystate. protoplanet.Themuchlessnoisyvolumeaveragedmagnetic
stresses are also plotted in figure 5. These behave similarly
When a protoplanet is placed in the centre of the box foreachofthesimulationsBa1,Ba2,Ba3andBa4butthey
waves are excited that cause outward angular momentum then separate at about time 500. Although simulation Ba4
transport. These waves also have an associated Reynolds’ which has the highest protoplanet mass also tends to have
stress (e.g. Papaloizou & Lin 1984) and may affect the un- the highest volume averaged magnetic stress, there is no
derlying turbulence. other such correlation between magnetic stress and proto-
Thevolumeaveragedtotalstressparameter α isplot- planet mass and so thismay not be significant.
h i
tedasafunctionoftimeformodelsBa1,Ba2,Ba3,andBa4 Although, as we have indicated above, strong fluctu-
in figure 5. The plots commence from the moment when ations are indeed the rule, there are definite behavioural
the protoplanet was introduced. The introduction causes trendsastheprotoplanetmassincreases.Toillustratethese,
the excitation of a wave and an increase in the Reynolds’ densitycontourplotsinatypical(x,y)planeforsimulations
stress which in turn increases the total stress. This is par- Ba1, Ba2 Ba3 and Ba4 near the end of the simulations are
ticularly noticeable in simulations Ba3 and Ba4 which have showninfigure6.Itcanbeseen thatasMpR3/(M∗H3)in-
MpR3/(M∗H3)=1.0and2.0respectivelyandsohavelarge creases, the wake becomes more prominent. It can be seen
enoughprotoplanetstoformagap.Thisbehaviourwasalso that this wake approximately coincides with that given by
found in thesimulation of a 5 Jupitermass planet in a disc equation (19) derived from the discussion about character-
(cid:13)c 0000RAS,MNRAS000,000–000