Three-dimensional Calculations of High and Low-mass Planets Embedded in Protoplanetary Discs

18 Pages·0.51 MB·English

by  M. R. Bate

Mon.Not.R.Astron.Soc.000,000–000 (0000) Printed2February2008 (MNLaTEXstylefilev2.2)
Three-dimensional calculations of high and low-mass
planets embedded in protoplanetary discs
M. R. Bate,1⋆ S. H. Lubow,2 G. I. Ogilvie,3 and K. A. Miller.4
1School of Physics, Universityof Exeter, StockerRoad, ExeterEX4 4QL
2Space Telescope Science Institute, 3700 San Martin Drive, Balitmore, MD 21218, USA
3
3Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA
0
4Astronomy Department, University of Maryland, College Park, MD 20742, U.S.A.
0
2

n
Accepted byMNRAS
a
J

9
ABSTRACT

1
v Weanalysethenon-linear,three-dimensionalresponseofagaseous,viscousproto-
4 planetarydisctothepresenceofaplanetofmassrangingfromoneEarthmass(1M⊕)
5 to one Jupiter mass (1 MJ) by using the ZEUS hydrodynamics code. We determine
1 the gas flow pattern, and the accretion and migration rates of the planet. The planet
1 is assumed to be in a fixed circular orbit about the central star. It is also assumed
0 to be able to accrete gas without expansion on the scale of its Roche radius. Only
3 planets with masses M ∼>0.1 M produce significant perturbations in the disc’s sur-
p J
0
face density. The flow within the Roche lobe of the planet is fully three-dimensional.
/
h Gas streams generally enter the Roche lobe close to the disc midplane, but produce
p muchweakershocksthanthestreamsintwo-dimensionalmodels.Thestreamssupply
- materialtoacircumplanetarydiscthatrotatesinthesamesenseastheplanet’sorbit.
o
Much of the mass supply to the circumplanetary disc comes from non-coplanar flow.
r
t The accretion rate peaks with a planet mass of approximately 0.1 M and is highly
s J
a efficient, occurring at the local viscous rate. The migration timescales for planets of
v: masslessthan0.1MJ,basedontorquesfromdiscmaterialoutsidetheplanets’Roche
lobes,areinexcellentagreementwiththelineartheoryofTypeI(non-gap)migration
i
X for three-dimensional discs. The transition from Type I to Type II (gap) migration
is smooth, with changes in migration times of about a factor of 2. Starting with a
r
a core whichcan undergo runawaygrowth,a planet cangainup to a few M with little
J
migration. Planets with final masses of order 10 M would undergo large migration,
J
which makes formation and survival difficult.
Key words: accretion, accretion discs – hydrodynamics – planets and satellites:
formation – planetary systems: formation – planetary systems: protoplanetary discs.
1 INTRODUCTION riedoutintwodimensionsinordertounderstandbetterthe
dynamics of a circular orbit planet embedded in a gaseous
Young planets interact with the surrounding discs from
disc (e.g. Bryden et al. 1999; Kley 1999; Lubow, Seibert &
which they form by accreting mass and exerting torques.
Artymowicz1999;Nelsonetal.2000).Thesestudiesconcen-
Inthecore-accretion modelforplanet formation, theinitial
trated on the interaction between a 1 MJ planet and a gas
growth proceeds as solids accumulate to form a planetary
disc that orbit a 1 M⊙ star. The tidal forces caused by the
core(Lissauer1995,Wuchterl,Guillot&Lissauer2000).The
planet create a gap in the disc. In spite of the presence of
coremaybecomeaterrestrialplanetormayfurtherdevelop
the gap, disc mass flow on to the planet continues through
tobecomeagasgiantthroughaccretionofgas.Thetorques
thegapwithhighefficiency.Nearlyalltheflowthroughthe
resultingfromplanet-discinteractionscauseaplanettomi-
gap is accreted by the planet at a rate comparable to the
grate typically inwards (Lin et al. 2000; Ward and Hahn
rate at which mass accretion would occur in thedisc in the
2000).
absence of a planet.
Numerical hydrodynamical calculations have been car-
The flow within the Roche lobe of the planet is highly
non-axisymmetric and involves shocks produced by collid-
⋆ E-mail:[email protected] ing gas streams (Lubow et al. 1999; D’Angelo, Henning &
2 M. R. Bate et al.
Kley 2002). Torques on the planet are exerted by circum- the usual α prescription. The origin of the coordinate sys-
stellar disc material that lies outside the gap. In addition, tem is taken to be the centre of the star. (We ignore the
torques are exerted locally by the material that flows close slight centre-of-mass shift caused by the planet.) The disc
to the planet. The net torque results in inward migration. self-gravity is ignored. The flow is modelled in the orbital
Gasaccretioncancontinuetoplanetmassesoforder10MJ, frame that rotates with the angular speed of the planet
at which pointtidalforces aresufficientlystrong toprevent Ωp = GM∗/rp3,whereM∗ isthemassofthestar,rp isthe
flow into thegap. orbitaplradiusoftheplanet,Gisthegravitational constant,
These earlier studies were limited by their neglect of andwehaveneglectedthemassoftheplanet.Inthisframe,
effects in the vertical direction (perpendicular to the orbit the flow achieves a near steady state. We adopt spherical
plane). The Roche lobe radius of a 1 MJ planet orbiting a coordinates(r,θ,φ)withassociatedflowvelocitiesinthero-
1 M⊙ star is comparable to the local disc thickness. The tating frame u = (ur,uθ,uφ). The equations of motion for
effectsin theverticaldirection are evenmoreimportant for thedisc are
lower-mass planets. In addition, sufficiently low-mass plan- ∂ρ
etsdonotopenagapinthedisc.Analyticmodelsofplanet +∇·(ρu)=0, (1)
∂t
migration inthenon-gap case, sometimes called TypeI mi-
gTrhaetiroens,ulwtseriemcpalirerdiedthaotutpliannetwtaorydmimigenrastioionnst(iWmeasrcdal1e9s9a7r)e. ∂∂str +∇·(sru) = ρrsin2θ rsuiφnθ +Ωp 2+ ρur2θ
(cid:16) (cid:17)
much shorter than the disc lifetimes. Planets would then
∂p ∂Φ
be accreted by the central star. Alternatively, the planets −∂r −ρ∂r +fr, (2)
could reside at the circumstellar disc inner edge, if there is
acentral hole in thedisc. This situation poses problems for ∂s u 2
planet survival and for gas giant planet formation. Recent ∂tθ +∇·(sθu) = ρr2sinθcosθ rsiφnθ +Ωp
(cid:16) (cid:17)
two andthreedimensional analytic calculations byTanaka,
∂p ∂Φ
Takeuchi, and Ward (2002) of the migration rates of low- −∂θ −ρ∂θ +rfθ, (3)
mass planets obtained smaller values (by about an order of
and
magnitude), making survivalmore plausible.
The two-dimensional numerical results also indicated ∂sφ +∇·(s u)=−∂p −ρ∂Φ +rsinθf , (4)
that accretion within the Roche lobe of a 1 MJ planet is ∂t φ ∂φ ∂φ φ
drivenbyshocks(Lubowet al. 1999; D’Angeloet al. 2002).
where ρ is the gas density, p is the gas pressure, sr = ρur
It is not clear whether the same flow pattern would persist
is the radial momentum per unit volume, s = ρru is the
θ θ
inthreedimensions.Recently,Kley,D’Angelo,andHenning
meridional momentum per unit volume, s = ρrsinθ(u +
φ φ
(2001) computed the three-dimensional flow for 0.5 and 1
Ω rsinθ)istheazimuthalangularmomentumperunitvol-
p
MJ planetsandfoundminordifferencesintheaccretionand ume,Φisthegravitational potentialduetothecentralstar
migrationratescomparedtothetwo-dimensionalcase.How- and the planet, and f = (fr,fθ,fφ) is the viscous force per
ever, theflow within theRoche lobe was unresolved.
unitvolumethatdescribestheeffectsofdiscturbulence.We
In this paper, we analyze the non-linear interactions usean unsoftened gravitational potential
between a planet and a gaseous disc by using global three-
dimensional numerical simulations that resolve the flow Φ(r)=−GM∗ − GMp , (5)
within the Roche lobe of a 1 MJ planet. We are interested r |r−rp|
in determining the flow patterns, the accretion rates, and where Mp is the planet mass. This is possible because the
migration rates for planets whose masses range between 1 location of the planet is such that all gravitational force
M⊕ (Earth mass) and 1 MJ (Jupiter mass). Each planet evaluations are made at a finite distance from the planet.
is assumed to be in a circular orbit about the central star. Wealso performed sometest calculations with gravity soft-
Themass andorbital radiusof each planet arefixedduring ened on the length scale of twice the grid resolution near
the simulation; consequently, the results do not include the the planet, but found no significant difference between the
effects on theflow of planetary migration. softened and unsoftened results.
The outline of the paper is as follows. Section 2 de- Equations (1), (2), (3), and (4) express conservation
scribesthecomputationalprocedure.Section3providesthe of mass, radial momentum, meridional momentum, and az-
results. In Section 4, we discuss the implications of the re- imuthal angular momentum, respectively. Equation (4) is
sultsforgiantplanetformation.Section5containsourcon- written in terms of the total azimuthal angular momentum
clusions. s ,ratherthanthatintherotatingframe.Thereasonisthat
φ
the s equation provides better numerical stability (Kley
φ
1998).
The equation of state istaken to belocally isothermal,
2 COMPUTATIONAL METHOD p∝ρT,withthetemperatureexpressedasaspecifiedfunc-
tionofradius,T(r).Thisequationofstateisappropriatefor
2.1 Basic equations
a gas that radiates internal energy gained by shocks with
We use a computational method that is similar to that high efficiency. The viscosity force f is assumed to be the
of Lubow et al. (1999), except that they performed two- standard Navier-Stokes force (see eq. [15.3] of Landau &
dimensional vertically averaged calculations whereas we Lifshitz1975; Klahr,Henning&Kley1999). Thecoefficient
solve theproblem in threedimensions. of shear viscoity µ represents the effects of disc turbulence,
Weassumeaviscousmodelforthediscturbulence,with whilethebulkviscositycoefficientζ issettozero.Thevalue
Three-dimensional calculations of planets in discs 3
for thekinematicturbulentviscosity ν =µ/ρ isassumed to betweentwoneighbouringzonesistoolargeoriftherandφ
beconstant inspaceandtime.Itcan beexpressed interms dimensionsofzonesaretoodifferent,numericalinstabilities
of the usual α prescription of Shakura & Sunyaev (1973). occur.Theinstabilitiesareworstintwo-dimensionalcalcula-
Namely,foradiscwithlocal isothermal soundspeedcs and tions,butalsoappearinthree-dimensionalcalculationsthat
verticalscaleheightH,dimensionlessparameterαisdefined havethesame gridding in r and φ.Wealso performed two-
through dimensionalcalculationswithincreasingresolution inr and
ν φ to test for convergence. Lowering the resolution in either
α(r)= . (6)
csH r or φ tends to increase the accretion rates. Convergence
testing of three-dimensional calculations was not practical
The above equations are non-dimensionalized so that
due to the factor of 16 increase in computational time that
the unit of time is the inverse of the planetary orbital fre-
would havebeen required to doubletheabove resolution in
quency Ωp, the unit of distance is the orbital radius of the all threedirections. Theaboveresolution was chosen as the
planet rp, and G=1. minimumrequiredbothtoavoidnumericalinstabilities and
togivequasi-steady-stateaccretionratesthathadconverged
towithinafewpercentinthetwo-dimensionalcalculations.
2.2 Numerical method
The Roche lobe of the planet had a radius
The equations are solved using a three-dimensional spheri-
cal coordinate version of the ZEUS-2D code (Stone & Nor- rR = Mp 1/3rp. (7)
man 1992). The code was written by K. A. Miller and (cid:16)3M∗(cid:17)
J. M. Stone. It was modified to include a standard three- We modelled planets from 1 Earth mass (3×10−6 M⊙) to
dimensional Navier-Stokes viscous force term. The code al- 1 Jupiter mass (1×10−3 M⊙). With the mass of the star
lows for variably-spaced gridding, which permits us to ob-
equalto1M⊙,theplanets’RocheradiirangedfromrR/rp =
tain higher resolution in the vicinity of the planet. The
0.010 to 0.069. Thus, the Roche radius of the Jupiter-mass
timestepssatisfy theusualCourant condition,for which we
planetwasresolvedby24zones(i.e.theRochelobecontains
have adopted Courant number 0.3. The code provides an about 3×104 zones), while that of the Earth-mass planet
artificial viscosity term, but since we introduce a Navier-
was marginally resolved by 3.5 zones (i.e. the Roche lobe
Stokesviscous force, we suppressthis artificial viscosity. Of
contains about 80 zones).
course,thereissomeintrinsicnumericalviscositybecauseof
The planet was assumed to beable to accrete material
thefinitegridding.Forauniformlyspacedmesh,thecodeis
withoutsubstantialexpansionofitsradiusonthescaleofits
formally second-order accurate in space and first-order ac-
Roche lobe. Some models of planet evolution suggest that
curate in time. For variably spaced meshes, it is formally
thisassumptionisvalidonlyforplanetsofmassgreaterthan
first-order accurate in space. However, a high level of accu-
about 50 M⊕ (Pollack et al 1996). However, these models
racyandresolutioncanbeattainedbylimitingthefractional
are subject to number of major simplifying assumptions,
change in mesh spacing between adjacent cells to be small,
as discussed in Wuchterl et al. (2000). We return to this
of order 1%. Thecode uses van Leerinterpolation.
point in Section 4. In any case, the accretion assumption
providesasimpleprescriptionforhandlingtheaccretionflow
that allows us to compute torques on the planet due to its
2.3 Numerical grid and initial conditions
interaction with thedisc.
The planet was fixed at location (r,θ,φ) = (1,π/2,π). To simulate the accretion on to the planet, we nearly
We modelled the disc in the region r ∈ [0.3,4.0], θ ∈ fully removed material in the four grid zones that sur-
[π/2−4H/r,π/2], and φ ∈ [0,2π]. We imposed reflective roundedthelocationoftheplanetineachtimestep.Aresid-
boundary conditions at the radial and θ grid boundaries ual density was retained in these zones to avoid numerical
and periodic boundary conditions at the azimuthal bound- divergences. The removed material was assumed to be ac-
ary. The radial boundaries were sufficiently far from the cretedontotheplanet,althoughthemassoftheplanetwas
planet that the reflected waves were not noticeable. The not increased. There was a pressure force directed toward
grid was uniform in θ, but non-uniform in r and φ. Both the planet at the edge of the evacuated region. This force
ther andφgridswereuniforminthevicinityoftheplanet, wassmallcomparedwithotherdynamicalforcesformassive
but the grid spacing increased logarithmically away from planets, but was noticable for the lowest mass planets that
thisuniform region (i.e. each zonewas a certain percentage we modelled (see Section 3.2).
larger or smaller than the zone preceding it). In addition, AtemperatureprofileT(r)∝r−1wasusedatalltimes.
the last two zones inside φ = 0 and φ = 2π were uni- Itwasnormalised suchthat H/rp=cs/(Ωprp)=0.05. This
form to allow periodic boundaries to be implemented eas- profilegivesaconstantvalueofH/r throughouttheunper-
ily. Each region of the radial grid had [0.3,0.8965] = 28, turbeddisc.Itisnumericallyconvenientwhenusingaspher-
[0.8965,1.1035] = 72, [1.1035,4.0] = 80 zones. The φ grid icalgrid,sincetheverticalresolutionofthediscisconstant.
had [0,0.19] = [2π − 0.19,2π] = 2, [0.19,π − 0.1035] = We set the kinematic disc viscosity to ν = 10−5 in our di-
[π+0.1035,2π−0.19]=106, [π−0.1035,π+0.1035] =72 mensionless units, which corresponds to α = 4×10−3 at
zones in each region. The θ grid modelled 4 scale heights r=rp =1 (equation 6).
above the disc midplane using 36 zones. Thus, the zones in The underlying initial disc density profile was
the vicinity of the planet had dimensions of 0.002875 in r chosen to be axisymmetric and follow ρ(r,θ,φ) ∝
andφ,andtwicethisvalueinθ.Thisgridwasdecidedupon r−3/2exp[−(θ−π/2)2r2/(2H2)]. Thus, the disc surface
after extensivetestingusingtwoand three-dimensionalcal- densityΣ(r)∝r−1/2.ForplanetswithmassesMp ≥0.1MJ,
culations. In particular, we found that if the change in size we imposed an initial gap near the planet. The initial gap
4 M. R. Bate et al.
Figure 2. The final azimuthally averaged disc surface density
for 1 (long-dashed), 0.3 (dot-dashed), 0.1 (dotted), 0.03 (short-
dashed), and 0.01 (thin solid) MJ planets. Only planets with
masses Mp > 0.1 MJ (Mp > 30 M⊕) produce significant per-
turbations.T∼hethicksolidlin∼egivestheresultfora1MJ planet
fromthetwo-dimensionalcalculationsofLubowetal.(1999).
Figure1.Theaccretionrateversustimeforeachofthe6planet
calculations.Theaccretionrateisaveragedoverevery1/20ofan 0.0075 M⊙).Thisgivesan unperturbeddiscsurfacedensity
orbit and plotted in units of the disc mass (7.5 MJ) per orbit. of 75 g cm−2 at theradius of the planet.
Thecalculationswererununtiltheaccretionratesreachedquasi- The calculations were run until the accretion rate on
steadyvalues. Thelinesgivetheaccretion ratesforplanets with to the planet reached a quasi-steady state (Figure 1). The
masses of 1 (long-dashed), 0.3 (dot-dashed), 0.1 (dotted), 0.03
number of orbits required depended on the mass of the
(short-dashed),0.01(solid),and0.003(dot-long-dashed)MJ.Low planet.Lower-massplanetsreachedasteadystateafteronly
massplanetsreachsteadyaccretionratesafteronlyafeworbits,
whilethehighestmassplanetsrequireapproximately100orbits. about5orbits,whilethe1MJplanetrequiredapproximately
100 orbits.Thecalculations neverquitereach atruesteady
statebecausetheprotoplanetary discisevolvingduetothe
Navier-Stokesviscosity.Thus,theaccretionratesofsomeof
sizeandstructurewereestimatedbyanapproximatetorque
the high-mass planets can be seen to decrease slightly with
balance condition between viscous and tidal torques near
time after thequasi-steady-state accretion is reached.
the planet. For the Jupiter-mass planet, the density at the
The calculations were performed on the United King-
middle of the gap was 1% of the unperturbed density. For
dom Astrophysical Fluids Facility (UKAFF), a 128-
the0.3MJ and0.1MJ planetstheinitialgapdensitieswere
processor SGI Origin 3800 computer, and on GRAND, a
4.2% and 35% of theunperturbeddensity,respectively.
24-processorSGIOrigin2000computer.Thecomputational
In order to confirm that the presence of an initial gap
timerequiredforeachcalculationwasrelativelyindependent
didnotaffectthefinalresults,weperformedtwocalculations
of the planet’s mass and required approximately 160 CPU
ofthe0.1MJplanet,onewithaninitialgapandonewithout.
hours per orbit on UKAFF. The total CPU time required
Theaccretionandmigrationratesofthecalculationwiththe
forallthecalculationswasapproximately50000CPUhours.
initial gap reached steady values after about 30 orbits. The
calculation without an initial gap initially had muchhigher
accretion and migration rates, but these converged to the 3.2 Density structure and gas flow
valuesgivenbytheothercalculationafterapproximately100
orbits.Thus,thepresenceofaninitialgapforthehigh-mass The interaction of the planet with the disc alters the den-
planets does not affect the final results, but it significantly sity and flow of the gas in the vicinity of the planet. This
reduced thecomputational time. interaction was analysed by means of two-dimensional sim-
ulations for high-mass planets by Lubow et al. (1999). Re-
cently,D’Angeloetal.(2002)studiedthegasflownearlow-
mass planets by using two-dimensional simulations. Here,
3 RESULTS westudytheprobleminthreedimensionsforbothhighand
low-mass planets.
3.1 Calculations
We performed calculations of 6 planets with masses of 1,
3.2.1 Global features
0.3, 0.1, 0.03, 0.01, and 0.003 Jupiter masses, MJ (i.e. 330,
100, 33, 10, 3.3, and 1 Earth masses, M⊕). The results are Figure2plotstheazimuthallyaveragedsurfacedensityasa
scaledsothattheplanetisatadistanceof5.2AUfroma1 functionofradiusfortheplanetswithmassesrangingfrom1
M⊙ star.Thediscmass,Md,between1.56and20.8AU(the to0.01MJ.Alsoplotted(withthethicksolidline)isthetwo-
boundariesofthegrid)istakentobe7.5Jupitermasses(i.e. dimensional result for a 1 MJ planet after 150 orbits from
Three-dimensional calculations of planets in discs 5
Figure 1 of Lubow et al. (1999). The surface density in the densityperturbationsexpectedfromlinearresonancetheory
vicinity a 1 MJ planet is well modelled by two-dimensional (see Figure 6 of Tanakaet al. 2002).
calculations (see also Kley et al. 2001). Only planets with Figure 6 provides the density of the gas in the vicinity
masses Mp ∼> 0.1 MJ produce significant perturbations in oftheplanet,butatheightH abovethemidplane.Figure7
the disc’s surface density; 0.03 MJ (10 M⊕) is insufficient. givesthedensityandvelocityvectorsinanr−zslicethrough
Thisresultisconsistentwiththeexpectationthatagapwill thediscatthelocationoftheplanet(i.e.φ=π).Theresults
starttoopenwhentheRocheradiusoftheplanetiscompa- areplottedbyusingaCartesiancoordinatesystem(x,y,z),
rable to the disc scaleheight H. The Roche radii of planets such that the the planet is located at (−1,0,0) and disc
with masses 0.3, 0.1, and 0.03 MJ are rR/rp =0.046,0.032, midplane is definedby z=0.
and0.022, respectively.ComparingthesetoH/r=0.05, we Noticethatwehavenotplottedstreamlinesinthesefig-
findthatrR >H/2isrequiredforasignificant surfaceden- ures.Itisdifficulttoplotasetofthree-dimensionalstream-
sity perturbation. Even for H = rR (the 0.3 MJ case), the lines because they cross each other in projection on to a
gap only dropsto 15% of theunperturbeddensity. graph. Streamlines can be plotted at the disc midplane
In Figure 3, we plot the global surface density of the wheretheverticalvelocityiszero(bysymmetry).Theveloc-
disc at the end of the simulations for each of the 6 plan- ityvectorsinFigure7cannotbeconnectedtoformstream-
ets. The (partial) clearing of gaps by the planets is clearly lines because information about the velocity in the depro-
visible in the global surface density plots for planets with jected direction (y) is required. Consequently, one cannot
masses Mp ≥0.1 MJ. Spiral density waves launched by the conclude that material rapidly drops to the midplane near
planet propagate inwards and outwards from the planet’s x = 1. Instead, material off the midplane flows past a low-
radiusintothedisc.FormassesMp ≥0.1MJ,thesedensity mass planet in they−direction.
waves are strong enough to appear even in the azimuthally Thesefiguresallowustodeterminethestructureofthe
averaged surface density profiles (Figure 2). The strength shocksoutsidetheRochelobeoftheplanet.Comparingthe
of the surface density perturbations decreases with lower locations of the shocks in similar panels between Figures 5
massplanets,buttheirformisessentiallyindependentofthe and 6, we find that the locations are similar for the low-
planet’smass.Forthelowest-massplanet,thespiraldensity mass planetsbut,for 1 and0.3 MJ, theshock frontson the
perturbations are too small to be visible over the underly- midplane lead the shock fronts off the midplane. Thus, a
ingdensitygradient inthedisc.Thesewavesareessentially section perpendicular to the shock front for the high-mass
two-dimensionalinnatureandhavebeenstudiedinthepast planetswouldhavea‘bowshock’shape,whereasforthelow-
usingtwoandthree-dimensionalnonlinearcalculations(e.g. massplanetstheshocksareessentiallyverticalplanes.These
Artymowicz 1992; Kley 1999; Bryden et al. 1999; Lubowet shock shapes can be seen clearly in Figure 7 with the tran-
al. 1999; Nelson et al. 2000; Kley et al. 2001) as well as sition from nearly vertical (0.1 MJ) to bowshock-shaped (1
analytically and through linearized numerical calculations MJ). The discontinuity in the velocity vectors across these
(Ogilvie & Lubow2002). shocks is also clear in Figure 7.
Returning to the density structure and streamlines at
themidplane(Figure5),weseethatgasatthesameradius
3.2.2 The flow in the vicinity of the planet
as the planet but outside the planet’s Roche lobe moves
In Figure 4, we plot the surface density in the vicinity of on horseshoe orbits, as was reported in earlier studies (e.g.,
the planet. Overlaid on the plots are the Roche lobes of Bryden et al. 1999; Kley 1999; Lubow et al. 1999). For the
theplanetsand theirinner(L1) and outer (L2) Lagrangian high-massplanets,thesehorseshoeorbitsoccupypartofthe
points (crosses). In Figure 5, we plot the density of the gap in the disc. The radial extent of the horseshoe orbits
disc(greyscale)andstreamlinesonthediscmidplaneinthe decreases as themass of theplanet is decreased.
vicinity of the planet. From these figures, we see that the Betweentheouterdiscandtheportionofthehorseshoe
waves are present in the circumstellar discs, as was found orbits ahead of the planet, and between the inner disc and
in the previous two-dimensional simulations. They are ini- theportion ofthehorseshoe orbitsbehindtheplanet,there
tiated near the planet and propagate as shocks. Even with are two streams of material that enter the Roche lobe of
the 1 M⊕ planet, the density jump across the shock is of the planet. This is the material (along with material from
order 10%, and themore massive planets havegreater den- above and below the disc midplane that is harder to visu-
sityjumps.Thediscontinuitiesinthestreamlinesacrossthe alise) that is accreted by the planet (see also Lubow et al.
shocks are obvious for planets with masses Mp ≥10 M⊕. 1999). For planets with masses Mp ≤ 0.1 MJ, we notice
The shape of the shocks is independent of the planet’s thatthebreadthofthesestreams generally increases as the
masssolongasadeepgapisunabletoform(i.e.forplanets planet’smassincreases,withthebreadthofthestreambeing
withmassesMp ∼<0.1MJ).Forthehigher-massplanets,the of order the Roche radius of the planet. In Section 3.3, we
shocksbecomemorecurved,theoutershockmovesforward, usethisobservationtodevelopamodelfortheaccretionrate
and the inner shock moves back. Unlike the 1 MJ case, the of the low-mass planets. For the 1 and 0.3 MJ planets, the
shocksgeneratedbythelow-massplanetsdonotextendra- streams are significantly narrower than the planet’s Roche
dially all the way to the planet’s Roche lobe. Instead, their lobe. In fact, for the 0.3 MJ planet, no streamlines on the
radialextentisdeterminedbythediscthickness.Theshocks midplane enter the planet’s Roche lobe. Since this planet
reach no closer to the planet than a distance of approxi- has the second highest accretion rate among the cases we
mately H (i.e. theybegin at approximately at r=rp±H). consider, material is still being accreted, but the flow must
This fact will be important when we come to consider the beintrinsicallythree-dimensionalwiththeaccretioncoming
torqueexerted on the planet by the disc. Although the dis- fromaboveandbelowthemidplane.The0.3MJ planetmay
turbancesarenon-linear,theirformissimilartothesurface be a special case, since the its Roche radius rR = 0.046 is
6 M. R. Bate et al.
Figure 3.Discsurfacedensities,Σ,for1,0.3,0.1,0.03,0.01,and0.003MJ planets (top-lefttobottom-right).
Three-dimensional calculations of planets in discs 7
Figure 4. Surface density, Σ, in the vicinity of the planet for 1, 0.3, 0.1, 0.03, 0.01, and 0.003 MJ planets (top-left to bottom-right).
AlsoplottedaretheRochelobes(whitecurves)andtheinner(L1)andouter(L2)Lagrangianpoints(crosses).
8 M. R. Bate et al.
Figure5.Discdensity,ρ(greyscale),andstreamlines(dashedandsolidlines)onthediscmidplanefor1,0.3,0.1,0.03,0.01,and0.003
MJ planets (top-left to bottom-right). The white solid streamlines are the critical streamlines that mark the boundaries between the
outer/innerdisc,theaccretionstreamsintotheRochelobe,andthehorseshoeorbits.Noticethatthetop2figuresgivethelogarithmof
thedensity.
Three-dimensional calculations of planets in discs 9
Figure 6. Same as Figure 5, but at distance H from the disc midplane. The white line indicates the size of the Roche lobe at the
midplane.
10 M. R. Bate et al.
Figure7.Discdensity,ρ,andvelocityvectorsinanr z slicethroughthediscatthelocationoftheplanet(i.e.φ=π)for1,0.3,0.1,
−
0.03,0.01,and0.003MJ planets(top-lefttobottom-right). Noticethatthetop2figuresarescaledbythelogarithmofthedensity.