Table Of ContentDraftversionMarch27,2017
PreprinttypesetusingLATEXstyleemulateapjv.12/16/11
SAVETHEPLANET,FEEDTHESTAR:
HOWSUPER-EARTHSSURVIVEMIGRATIONANDDRIVEDISKACCRETION
JeffreyFung1&EugeneChiang
DepartmentofAstronomy,UniversityofCaliforniaatBerkeley,CampbellHall,Berkeley,CA94720-3411
DraftversionMarch27,2017
ABSTRACT
Two longstanding problems in planet formation include (1) understanding how planets survive migration,
and (2) articulating the process by which protoplanetary disks disperse—and in particular how they accrete
7
1 ontotheircentralstars. Wecangoalongwaytowardsolvingbothproblemsifthediskgassurroundingplanets
0 has no intrinsic diffusivity (“viscosity”). In inviscid, laminar disks, a planet readily repels gas away from its
2 orbit. Onshorttimescales,zeroviscositygasaccumulatesinsideaplanet’sorbittoslowTypeImigrationby
ordersofmagnitude. Onlongertimescales,multiplesuper-Earths(distributedbetween,say,∼0.1–10AU)can
r
a torqueinviscidgasoutofinterplanetaryspace,eitherinwardtofeedtheirstars,oroutwardtobeblownaway
M inawind. Weexplorethispicturewith2DhydrodynamicssimulationsofEarthsandsuper-Earthsembedded
ininvisciddisks,confirmingtheirslow/stalledmigrationevenundergas-richconditions,andshowingthatdisk
4 transport rates range up to ∼10−7M yr−1 and scale as M˙ ∝ ΣM3/2, where Σ is the disk surface density and
(cid:12) p
2 M is the planet mass. Gas initially sandwiched between two planets is torqued past both into the inner and
p
outerdisks. Insum,sufficientlycompactsystemsofsuper-Earthscancleartheirnataldiskgas,inadispersal
] historythatmaybecomplicatedandnon-steady,butwhichconceivablyleadsoverMyrtimescalestolargegas
P
depletionssimilartothosecharacterizingtransitiondisks.
E
Subject headings: accretion, accretion disks — methods: numerical — planets and satellites: formation —
.
h protoplanetarydisks—planet-diskinteractions—circumstellarmatter—stars: variables:
p TTauri,HerbigAe/Be
-
o
tr 1. INTRODUCTION ahead of a migrating planet exerts a “feedback” torque that
s slows and can even stall migration (Hourigan & Ward 1984;
a Protoplanetarydiskshavetwojobs: makeplanetsandfeed
Ward & Hourigan 1989). Rafikov (2002b, see his equation
[ theirhoststars. Thefirsttaskisfrustratedbymigration: disk
torques force planetary orbits to decay (e.g., Kley & Nel- 53) calculated that for inviscid disks in which planet-driven
2 wavesdissipatebysteepeningintoshocks,thecriticalplanet
son2012),evacuatingtheveryregionswhereplanetsareob-
v massabovewhichTypeImigrationshutsoffis:
servedinabundance(∼0.1–10AUfromthestar;e.g.,Clanton
1
16 &DqurieGrseassiunadgim&2e0cC1hh4aa;nriBbsomunrkntoeeateurta2na0sl.p1o25r0)t.1aF5wu;laCfiylhltrihnisegtidtahinseskes’nsecaeontngdaullt.aa2rsk0m1ro5e--; Mcr (cid:39)4(cid:32)0h.p0/3r5p(cid:33)3(cid:32)MM∗(cid:33)Σp1r0p2/−M3 ∗153 M⊕, (1)
8 (cid:12)
mentum. Magnetictorquesarepromisingbutdependonseed
0 whererandharethediskradiusandscaleheight,Σisthedisk
fieldsofuncertainprovenance(e.g.,Bai2016).
.
1 Goodman & Rafikov (2001) proposed that the two prob- gassurfacedensity,M∗isthecentralstellarmass,andthesub-
0 lems are actually one: that planets themselves—if they can scriptpindicatesevaluationneartheplanet’sposition.Lietal.
7 survive migration—can provide an effective source of disk (2009)andYuetal.(2010)haveperformednumericalsimula-
1 tionssupportingtheanalyticcalculationsbyRafikov(2002b),
viscositybyexcitingdensitywavesthattransportangularmo-
: and confirming that super-Earths in low-viscosity disks mi-
v mentum outward. Sari & Goldreich (2004) emphasized that
gratemuchmoreslowly(anderratically)thanispredictedby
i suchplanetsmustbemassiveenoughtoopengapsandavoid
X TypeI.
Type I migration. They focused on giant Jupiter-mass plan-
Aplanetofmass M atr drivesadiskmasstransportrate
r ets, a demographic that is now understood to be rare (e.g., p p
a Cumming et al. 2008). In this paper, we turn our attention M˙ atdistancerof
tcoovseurpeedr-bEyaKrtehps:lebrotdoibeseroeflamtiavsesly1–co1m0mMo⊕ntphlaatchea(vee.gb.,eFernesdsiisn- M˙(r)=−2F0r∂ϕ(r) (2)
l ∂r
etal.2013).
where
disCkasnarseupseurffi-EcaiertnhtslyaivnoviidscTidy.peTIhemdigerpaetniodne?nceYoefs—thiefTthyepier F =Σ r2l Ω (cid:32)Mp(cid:33)2(cid:32)hp(cid:33)−3 (3)
Idriftrateondiskviscosityisperhapsunder-appreciated, as 0 p p p p M r
∗ p
itisnotexplicitlycalledoutinthetypicallyquotedTypeIfor-
measuresthetotalangularmomentumcarriedawaypertime
mula (see, e.g., Kley & Nelson 2012). Crucially, without an
by planet-driven waves (a.k.a. the total integrated one-sided
intrinsicdiskviscositytosmoothawaytheplanet’sperturba-
Lindblad torque), Ω is the orbital frequency, and l = Ωr2
tions to the disk’s surface density, a pile-up of disk material
is the specific angular momentum. The dimensionless func-
tion ϕ(r) describes how waves, as they travel away from the
email:jeff[email protected] planet,dampwithdistance,depositingtheirangularmomen-
1NASASaganFellow tumtodiskgasandtherebypropellingmaterialradially.From
2 Fung&Chiang
Rafikov(2002b,seehisequationsB1and33), 2. NUMERICALMETHOD
We use the graphics processing unit (GPU) accelerated
(cid:32) Mp (cid:33)−12 (cid:32)|rp−r|(cid:33)−54 hydrodynamics code PEnGUIn (Fung 2015) to perform 2D
ϕ∼ , (4)
M h simulations of disk-planet interactions. It is a Lagrangian-
thermal p
remapshock-capturingcodethatusesthepiecewiseparabolic
valid for M (cid:46) M ≡ (h /r )3M and ϕ (cid:46) 1 (i.e., dis- method (Colella & Woodward 1984) to solve the continuity
p thermal p p ∗
tancesfarenoughfromtheplanetthatthewavesaredissipat- andmomentumequations:
inginweakshocks),andwherewehaveignoredorder-unity DΣ
constantsandallradialvariationsinh,Σ,andgassoundspeed. =−Σ(∇·v) , (7)
Itfollowsthat Dt
Dv 1 1
M˙(r)∼sign(r−r )Σ r2Ω (cid:32)Mp(cid:33)32 (cid:32)hp(cid:33)−52 (cid:32)|rp−r|(cid:33)−94 . (5) Dt =−Σ∇p+ Σ∇·T−∇Φ, (8)
p p p p M∗ rp hp whereΣisthegassurfacedensity, vthevelocityfield, pthe
verticallyaveragedgaspressure,TtheNewtonianstressten-
Atr<rp,M˙ asgivenbyEquation(5)isnegative(massflows sor, and Φ the combined gravitational potential of the star
inward),andviceversa;aplanettendstorepelmaterialaway and the planet(s). We use a globally isothermal equation
fromitself. Toavoidminussigns,wewillignorethisformal of state: p = c2Σ with a spatially constant sound speed
signconventionsothatallourreportedvaluesfor M˙ willbe c = 0.035v s(cid:39) 1 km s−1 where v is the Keplerian
s K,1AU K,1AU
positive and understood to be inward unless otherwise indi-
velocityat1AUarounda1M star. (Thisc correspondstoa
(cid:12) s
cated.
disktemperatureof300Kassumingameanmolecularweight
Notehow M˙ ∝ M3/2 andnot M2. AlthoughthetotalLind-
p p of2.34.)
blad torque scales as M2 (Equation 3), that torque is dis- Inapolarcoordinatesystem(radiusr,azimuthφ)centered
p
tributedoveradistancethatincreaseswithdecreasing M (as onthestar,
p
Mp−2/5,ascanbeseenbysolvingfor|rp−r|intermsof Mp at GM (cid:88)Np
fixedϕinEquation4). Thusatfixeddistanceawayfromthe Φ=− ∗ + Φ (9)
planet,M˙ increaseswithM withapowerlessthan2. r p,i
Inserting M = 10M anpd other nominal parameters (for i=1
p ⊕ GM GM rcosφ(cid:48)
rM<˙(rrp))∼in1to0−(85)1y0iΣe−lp3drMsp2 (cid:32)2π1/yΩrp(cid:33)−1(cid:32)3M×p/1M0−∗5(cid:33)23 (cid:32)0h.p0/3r5p(cid:33)−14 Φp,i =−(cid:113)r2+rp2,i−2rrpp,,iicosφ(cid:48)i +rs2,i + p,irp2,i (i10)
(cid:12)
whereGisthegravitationalconstant,M =1M isthestellar
(cid:32)(rp−0.r5)/rp(cid:33)−49 M(cid:12)yr−1, (6) mofapslsa,nthetes,suMbpscirsipthteiplalabneelstemaacshs,pΦlapnetht,eN∗pplainsetht’es(cid:12)tgortaavlintautmiobnearl
potential, r the planet’s radial coordinate, r the smoothing
p s
comparable to accretion rates measured for classical T Tauri lengthoftheplanet’spotential,andφ(cid:48) =φ−φptheazimuthal
stars (e.g., Calvet et al. 2005; Hartmann et al. 2006; Sicilia- separation from the planet. The stress tensor T is propor-
Aguilaretal.2010;Inglebyetal.2013). NotehowweaklyM˙ tional to the kinematic viscosity ν. Most of our simulations
dependsonh /r ,underscoringhowM˙ doesnotscalesimply are of inviscid disks with ν = 0. For our viscous disk sim-
p p
as the total Lindblad torque in Equation 3 (which scales as ulations, we use ν = αcsh, where the Shakura-Sunyaev pa-
(hp/rp)−3), but depends also on the distance over which that rameter(cid:112)α = 0.001, h = cs/ΩK is the local scale height, and
torqueisexerted,aswehavedescribedabove. Ω = GM /r3 is the Keplerian orbital angular frequency.
K ∗
The above considerations indicate that with super-Earths Atr =1AU,h/r =0.035. Wesetr =0.5h,asisappropriate
s
we might have our cake (survive migration) and eat it, too for2Dsimulations(Mu¨lleretal.2012).
(drive disk accretion). Of course, a single super-Earth is in- A given planet feels the gravitational force from the star,
sufficientbecauseitsreachistooshort(M˙ dropsas|r−rp|−9/4). the disk, and other planets. The disk force on the planet is
Multiplesuper-Earthsareneededtoshuttletheaccretionflow calculatedbydirectsummationoverallmasselementsinthe
fromdistancesofafewAUdowntothestellarradius. Real- disk, withthe“background”axisymmetriccomponentofthe
itywillbenon-steadyandlikelymessy(see,e.g.,Figure4of disksurfacedensitysubtractedoff. Becausethediskdoesnot
Rafikov2002a), withmaterialbetweenadjacentplanetshav- feel its own gravity at all (i.e., we ignore disk self-gravity;
ingafatethatisnotobvious: doesthesandwichedgasdrain seeequation9), eliminatingthisaxisymmetriccomponentin
inward,ordoestheinnerplanetholdbackmaterialpushedin- thedisk-planetforcingimprovesconsistencybetweenthemo-
wardbytheouterplanet? Andtowhatextentdosuper-Earths tionsoftheplanetsandthedisk. Planetmigrationshouldbe
migratewiththeaccretionflowtheydrive? minimally affected by this procedure, since the background
Hereweexplorethesequestionsusingfullynon-linear,2D componentofΣexertsnotorque.Spuriousforcesarisingfrom
hydrodynamicalsimulationsofsuper-Earthsembeddedinin- withintheplanet’sHillspherearesometimesaconcernifthis
visciddisks. Wemeasurethemigrationhistoriesrp(t)andac- regionisunder-resolved. TheHillradius,rH =(Mp/3M∗)1/3,
cretionratesM˙ insimulationscontaining1or2super-Earths, rangesfrom0.3to0.6h ,similartothesmoothinglengthr .
p s
experimenting with varying the disk surface density and the Wehaveverifiedthatthetorquegeneratedwithinaradiusof
planet mass to test Equation 6. Our numerical methods are 0.5r fromtheplanetisnegligible, andsowedonotexcise
H
given in Section 2. Results are presented in Section 3 and the Hill sphere in force calculations. The planets’ motions
placedintobroadercontextinSection4. are integrated using a kick-drift-kick leapfrog scheme, with
SavethePlanet,FeedtheStar 3
1.2
thedriftstepoccurringsynchronouslywiththehydrodynam-
icsstep;i.e.,theplanets’positionsarelinearintimewithina
hydrodynamicsstep. 1.1
]
U viscous disk: fast migration
A
[ p 1 inviscid disks: slow/stalled migration
r
n
o 0.9
2.1. Initialandboundaryconditions,andgridparameters siti
o
p
Table1liststheparametersusedbyour6models. Thedisk et 0.8 type I migration
isinitializedwithapower-lawsurfacedensity: an single planet; viscous
pl single planet; inviscid
Σ=Σ0(cid:18) r (cid:19)−32 . (11) 0.7 ionunteerr ppllaanneett;; iinnvviisscciidd
AU
0.6
Weconsiderbothgas-richdiskshavingΣ0 =8.5×103gcm−2 0 1000 2t0im00e t [ yea3r0s0 ]0 4000 5000
resembling the minimum-mass extrasolar nebula (Chiang &
Laughlin 2013), and gas-poor disks havinga surface density Fig.1.—TypeImigrationcanbedefeatedforsuper-Earthsininvisciddisks.
Shownherearetheorbitalradiivs.timeofplanetsinthreeofoursimulations:
1000× lower. The initial velocity field is axisymmetric and
Keplerian,withcorrectionsfromgaspressure: asingle10M⊕planetinagas-rich(Σ0=8.5×103gcm−2)viscous(α=10−3)
disk(blue,model#1);asingleplanetofthesamemassinagas-richinviscid
(cid:114) (α = 0)disk(black, model#2); andtwosuchplanetsinagas-richinvis-
1 dp ciddisk(red+magenta,model#3). Thebluedashedlineisthetheoretically
Ω= Ω2 + . (12) expectedtrajectoryfromTypeImigration(theintegralofEquation13with
K rΣ dr C=2).
One planet, whose mass is increased gradually over the first boundaries. Ourresultscanonlybetrustedtotheextentthat
10yrofthesimulationtothefullvalueof M (either1,3,or theseboundarypile-upsdonotinterferewithplanetmigration
p
10M ),isplacedinitiallyatr=r (either1or0.75AU)and and disk accretion. We therefore limit ourselves to studying
⊕ p,1
φ = π. Intwo-planetmodels, weplaceasecondplanetof only the first few thousand years of planet-disk interactions,
p,1
equalmasstothefirstatr = r (either1.2or1.05AU)and beforeboundaryeffectsbecometoolarge.
p,2
φ =πinitially.
p,2
Our simulation grid spans the full 2π in azimuth, and ex- 3. RESULTS
tends from an outer radius of 1.8 AU to an inner radius r We assess to what extent planets migrate in inviscid disks
in
thatequalseither0.4or0.3AUdependingonwhetherr =1 (§3.1), and study how planet-driven accretion rates evolve
p,1
AUor0.75AU(seeTable1). Griddimensionsare800(r)× with time and depend on disk and planet masses (§3.2). For
3200(φ)whenr =0.4AU,and960×3200whenr =0.3 planet migration, models #1–3 demonstrate differences be-
in in
AU.Cellsarespacedlogarithmicallyinradiusanduniformly tween viscous and inviscid disks, and between single-planet
in azimuth. Our choices yield a resolution of ∼18 cells per andtwo-planetsystems. Forplanet-drivenaccretion,wevary
scale height h in both directions at r = 1 AU (similar to the diskandplanetmassesinmodels#3–6totestEquation6.
resolution of Li et al. 2009). Simulations at twice our stan-
dard resolution did not produce significant changes in either 3.1. Planetmigration
planetmigrationordiskaccretionrateforthefirst100yr. We Figure1plotstheorbitalevolutionofplanets,eachofmass
alsotestedourinvisciddiskmodelwithoutaplanet,andfound 10M , in the gas-rich disk models (#1–3). Overplotted for
thatthenumericalnoisein|M˙|wasabout3ordersofmagni- comp⊕arison is the trajectory expected from integrating the
tude below planet-driven disk accretion rates, corresponding TypeImigrationrate,
toanumericalviscosityofα<10−5.
stuRdayd.iAalftbeoruenxdpaerryimceonntidnigtiownisthraeqfueiwrewsapyescitaolmceaaresurinedthisiks r˙p,TypeI =−2CrpΩpΣMprp2(cid:32)MMp(cid:33)2(cid:32)hrp(cid:33)−2 , (13)
p ∗ p
accretionrates,wefoundthatthemoststablemethodwasto
trackthetotaldiskmasswithinacylinderofradius0.6AU— using the unperturbed surface density law in Equation 11 to
adistanceintermediatebetweentheinnermostplanetandthe evaluate Σ (using the actual surface density in the viscous
p
innerdiskboundary—whilepreventingmassfromleavingthe disk simulation #1 would give practically identical results,
grid. We adopt “zero flux” boundary conditions where mass since gaps do not form in that model). Three-dimensional
and momentum fluxes across the inner and outer disk edges simulations suggest that C ∼ 2–3 for our given disk profile
are always zero. In PEnGUIn, this is achieved by solving a (D’Angelo&Lubow2010;Fungetal.2017);thebluedashed
special Riemann problem at the boundaries, one where no curveinFigure1usesC = 2. Figure2displaysazimuthally
wavetravelstowardthesimulationdomain,andwherethera- averagedsurfacedensityprofilesatvariousepochsingas-rich,
dial velocity outside the domain is always zero. This imple- invisciddiskmodels#2and#3,withplanetlocationsmarked.
mentationconservedthetotalmasswithinthesimulationdo- In agreement with the simulations of Li et al. (2009) and
maintonumericalaccuracy. Theaccretionrate M˙ atr = 0.6 Yuetal.(2010),theplanetmigrationrateinourviscousdisk
AU is calculated by following over time the disk mass en- (model #1) is similar to the Type I rate (punctuated by what
closed, M (t). Because the function M (t) fluctuates appeartobeepisodesofevenfasterTypeIIImigration;Mas-
0.6AU 0.6AU
strongly, we fit independent lines to segments of data each set&Papaloizou2003;Peplinski2008),andmuchslowerin
lasting20yr,takingM˙ fromthebest-fittingslopes. invisciddisks(models#2and#3). Theinitiallyrapidmigra-
Asaplanetrepelsmaterialawayfromitsorbit,ourbound- tion seen in the inviscid simulations at t (cid:46) 200 yr is a tran-
ary conditions result in gas piling up at the inner and outer sient that decays after disk surface densities adjust to plan-
4 Fung&Chiang
TABLE1
ModelParameters
Model# Mp(M⊕) rp,1(AU) rp,2(AU) Σ0(gcm−2) α rin(AU) tend(years)
1 10 1 – 8.5×103 10−3 0.4 700
2 10 1 – 8.5×103 0 0.4 5000
3 10 1 1.2 8.5×103 0 0.4 5000
4 10 0.75 1.05 8.5 0 0.3 2000
5 3 0.75 1.05 8.5 0 0.3 2000
6 1 0.75 1.05 8.5 0 0.3 2000
Note. —tendistheendtimeofasimulation,inunitswheretheKeplerianorbitalperiodat1
AUis1year. Also,rp,1andrp,2aremerelytheinitialplanetlocationsatt=0;theplanetsare
completelyfreetomigrateinthesimulations.
Fig.2.—Left:Surfacedensityprofiles,azimuthallyaveragedandnormalizedagainsttheinitialpower-lawprofile,measuredat1000-yrintervalsforoursingle-
planet,inviscid,gas-richdiskmodel(#2).Right:Analogousprofilesforourtwo-planet,inviscid,gas-richdiskmodel(#3).Circles(triangles)markthelocations
oftheinner(outer)planet.
etary Lindblad torques, i.e., after the surface density pile-up
ahead of the planet attains a fractional amplitude on the or-
der of unity. After this initial adjustment period, migration
slowsandevenstallsattimes,withradialpositionschanging
by∼10–30%,orless,overkyrtimescales.
In the single-planet, gas-rich simulation (black curve in
Figure 1, model #2), the planet journeys slowly inward for
the first ∼2000 yr and practically stops in the mean from
t (cid:39) 2000–3500 yr, as disk gas that the planet has pushed in-
ward to r (cid:39) 0.7–0.8 AU piles up (blue and green curves in
theleftpanelofFigure2)andstymiesfurthermigration. At
t (cid:39) 3500 yr, the planet experiences a sudden drop in orbital
radius; we traced this drop to a close encounter between the
planet and a vortex formed at its outer gap edge at r (cid:39) 1.1
AU (see Figure 3). Some time after the encounter, the vor-
texgraduallydisperses,completelydecayingawaybytheend
of our simulation at 5000 yr (Figure 4, left panel). Similar
planet-vortex interactions were found in simulations by Lin
& Papaloizou (2010) and Yu et al. (2010). Thereafter, the
planet’smigrationreturnstoitsnear-zeromeanpace.
Inourtwo-planet,gas-richsimulation(#3),theouterplanet
stalls for the first ∼3000 yr (magenta curve in Figure 1), ap-
parently trapped at a local surface density maximum created
Fig.3.—Surfacedensitysnapshotofoursingle-planet,gas-richsimulation
(model#2)at3500years.Greendashedcirclesareat0.5,1,and1.5AU.The bytheinnerplanetwhoseforcingdominates:seethemagenta
bluedotindicatestheplanet’sposition. Thissnapshotistakennearthetime curve in the right panel of Figure 2, and note how similar it
oftheplanet-vortexcloseencounter. Thevortex,ofapproximatemass100 is to the corresponding magenta curve in the left panel for
M⊕,isseenataround1.1AU.Aftertheencounter,theplanetisperturbed thesingle-planetcase. Gaspushedinwardbytheouterplanet
radiallyinward(Figure1).Thevortexultimatelydecaysaway(Figure4).
strengthensthetorqueontheinnerplanetandforcesthelatter
to migrate inward by ∼20% over the same time period (red
curveinFigure1);contrastthisbehaviorwiththestallingob-
SavethePlanet,FeedtheStar 5
Fig.4.—Surfacedensitysnapshotsofourinviscidgas-richsimulations,foroursingle-planetmodel(#2,left)andtwo-planetmodel(#3,right)cases,takenat
theirendtimeof5000years.Greendashedcirclesareat0.5,1,and1.5AU.Thebluedotsindicatetheplanets’positions.Inthesingle-planetsnapshot,notethat
thevortexseeninFigure3hascompletelydecayedaway.Nocomparablystrongvortexwasfoundinthetwo-planetsimulationatanytime,andthedensityfield
appearsmoreaxisymmetricthaninthesingle-planetcase.
servedinthesingle-planetcase(blackcurve). Eventually, at Thecomparisonisbestmadeatearlytimesofthesimulation,
t ∼ 4000 yr, the outer planet disperses the surface density t(cid:46)1000yr,beforeradialsurfacedensityprofilesbecometoo
maximum in its vicinity, and proceeds to migrate slowly in- distorted. Forthesingle-planet,gas-richmodel#2,Equation
ward, slowing down near t ∼ 5000 yr as it runs into sand- 6yieldsM˙ ∼3×10−8M yr−1.Thispredictioniswithinafac-
(cid:12)
wichedgas. Meanwhile,theinnerplanetultimatelycomesto tor of ∼2 of the simulated result at early times (left panel of
a near halt in much the same way that it does in the single- Figure5). Asforthecorrespondingtwo-planetmodel#3,we
planet simulation, having run into material that has piled up expectfrom(6)thatM˙ shouldbeonlyfractionallylargerthan
justinteriortoitsorbit. TherightpanelofFigure4showsthe forthesingle-planetcase;thesecondplanetisfartherremoved
finalsurfacedensitydistribution.Weemphasizethatthispile- fromwherewemeasureM˙ (r=0.6AU),andsomakesonlya
up is physical as it is located at r (cid:39) 0.6 AU, away from the ∼50%contributiontotheaccretionflowthereascomparedto
innergridboundaryofthesimulationatr =0.4AU.Thelat- theinnerplanet. ThisisapproximatelyconsistentwithFigure
in
terlocationhasitsownseparatepile-up,whichdoesnotgrow 5.
to significance over the limited duration of our simulations. Forourgas-poormodel#4,Equation6predicts M˙ ∼ 1.5×
Thesamestatementappliestotheoutergridboundary. 10−10M yr−1,againwithinafactorof2ofthesimulatedre-
(cid:12)
Wenotethatinnoneofthetwo-planetsimulationsdidwe sult (right panel of Figure 5). Scaling the planet mass M
p
observe the formation of a vortex like the one seen in our downbyafactorof10frommodels#4to#6should,accord-
single-planet simulation. This difference might be physical, ing to Equation 6, reduce M˙ by a factor of 103/2 (cid:39) 30. By
anddeservesattentioninfuturestudiesofplanet-vortexinter- comparison,Figure5showsafactorof∼20decreasebetween
actions. thesetwomodels;weconsiderthisacceptableagreementwith
Allotherfactorsbeingequal,lowerdiskmassesshouldlead theanalyticexpectation.Insummary,oursimulationssupport
to even slower planetary migration rates. This is confirmed the variousfunctional dependencies predictedby Equation 6
in our gas-poor simulations (models #4–6) which exhibit no towithinafactorof2.
measurablechangeinplanetmeanradialpositions. Thusthe Atlatertimes,t(cid:38)1000yr,weobservetimevariabilityinM˙
gas-poor simulations can be used to diagnose disk accretion causedbythedeepeningofplanetarygaps,andbyplanetmi-
rateswithoutthecomplicatingeffectsofplanetarymigration, gration. Thesevariationsarelimitedtofactorsofafew. The
aswediscussinthenextsubsection. simulations easiest to interpret are models #4–6 (right panel
ofFigure5)whichhavetoolittlediskgastodriveplanetmi-
3.2. Diskaccretion gration. Theinitialgradualdeclinein M˙ seeninmodel#4is
causedbythedeepeningofgapsopenedbyits10-M planets;
Figure5showsdiskaccretionratesasfunctionsoftimefor ⊕
over the course of 2000 yr, the gas density in the immediate
all our inviscid models. As described in §2.1, the accretion
vicinityoftheplanetsdecreasesbyafactorof∼5fortheinner
rate M˙ ismeasuredbytrackingthebuildupofdiskmassin-
planetandbyafactorof∼3fortheouterone. Models#5and
sider =0.6AU,alocationinteriortotheplanetsatalltimes.
#6exhibitsteadieraccretionrates,astheirplanetshavelower
Later,in§3.2.1,wetrackthemovementofmassinitiallybe- masseswhicharelesseffectiveatopeninggaps.
tweentwoplanets.
More complicated behavior is seen in the gas-rich simula-
Webeginbycheckingwhetheroursimulationsarecompat-
tionswhereplanetsmigratemoreappreciably.Comparisonof
iblewiththeanalyticexpectationfor M˙ givenbyEquation6.
6 Fung&Chiang
Fig.5.—Diskaccretionratesmeasuredatr =0.6AU,alocationinteriortoallplanetsatalltimes. Eachdatapointrepresentsanaverageover20years(see
§2.1).Left:Gas-richα=0simulations(model#2inblack,#3inred).Right:Gas-poorα=0models(#4–6inblack,red,andblue,respectively).Att(cid:46)1000
yr,accretionratesagreewithpredictionsfromEquation6towithinafactorof2.VariationsinM˙ atlatertimesreflectplanetmigration(comparewithFigure1)
anddeepeningofgaps.
Fig.6.—Trackingthematerialinitiallysandwichedbytwoplanets. Fromlefttoright,weshowresultsfrommodels#4,5,and6. Blackcirclesdenotethe
planets’locationswhichdonotchangeinthesegas-poorruns.Inallcases,thegasescapesovertimetoeithersideoftheplanetarypair.
Figures 1 and 5 reveals that increases in M˙ can be traced to containing the highest mass planets, some gas does concen-
planets moving inward, either gradually, as in the first 2000 tratealongthemidlinebetweentheplanetsina“shepherded”
yr of models #2 and #3, or suddenly, as in the planet-vortex ring,resultinginlessmassleakingoutofthesandwichedre-
encounteratt (cid:39) 3300yrinmodel#2. Decreasesin M˙ corre- gion. Looking at models #4–6 in Figure 6, we see no clear
spondtoplanetsopeninggapsuponmovingtonewlocations. trendbetweentherateatwhichsandwichedgasescapesand
planet mass. There seems to be a complicated confluence
3.2.1. Thefateofgasinitiallysandwichedbetweenplanets of effects in the sandwiched region. Lindblad torques act to
shepherd some gas while also opening gaps that reduce the
Disk accretion driven by planets would be impractical if
local gas density; and co-orbital torques allow gas to escape
materialresidingbetweenplanetswereunabletoescape. We
viahorseshoeorbits,whoselibrationtimesscaleonlyweakly
trackthissandwichedgasinsimulations#4–6,eachcontain-
ingapairofplanetswhichmigratenegligibly.Weassigneach with planet mass (tlib ∝ Mp−1/2; Paardekooper & Papaloizou
gas parcel a “passive scalar” η that equals 1 for gas initially 2009). A detailed analysis is deferred to another paper; for
locatedbetweenr =0.8and1AU(betweenthetwoplanets), now, we conclude that, at least for comparable mass planets
andis0everywhereelse. Gaselementscarryηasaconserva- withorbitalspacingsliketheonewehaveassumed,theinner
tivequantity. planetpresentsaporousbarriertomaterialpushedinwardby
Figure6followstheη-taggedgasbyplotting the outer planet. Apparently gas that is pushed by the outer
planet toward the inner planet can be shuttled past the latter
∂Mη =(cid:90) 2πηΣrdφ, (14) onhorseshoeorbits(andviceversa).
∂r
0
vs. r at various times. We find that the sandwiched gas is
torquedbothinwardandoutward, escapinginroughlyequal 4. SUMMARYANDDISCUSSION
amounts to the inside of the inner planet and to the outside Usinghydrodynamicalsimulations, wehavedemonstrated
oftheouterplanet. Theopposingtorquesfromthetwoplan- that super-Earths in inviscid disks can simultaneously avoid
ets do not in general cancel. The opposing torques from the typeImigration(Figure1)andpromotediskaccretion(Fig-
two planets do not in general cancel, although in model #4 ure 5) by driving density waves. Disk accretion rates mea-
SavethePlanet,FeedtheStar 7
suredfromoursimulationsverifyanalyticpredictions(Equa- hoststar. Oritmaybetorquedbyplanetssofaroutwardthat
tion 6) to within a factor of 2 . We observed gap opening itescapesfromthesystemaltogetherinaphotoionizedwind
and planet migration in inviscid disks to be modest and to (e.g.,Alexanderetal.2014).
introduce order-unity effects on the disk accretion rate. We Howmassiveagasdiskcanasetofsuper-Earthsdrain? A
alsofoundinourtwo-planetsimulationsthatmaterialinitially first consideration is that the disk surface density can not be
sandwiched between two planets leaks past both into the in- so large that the embedded planet mass M < M ∝ Σ5/13,
p cr p
nermostandoutermostdisks(Figure6). lest the planet migrate away. Based on the typical parame-
Our models omit a number of effects. Many of these are terslistedinEquation1,adiskcontainingΣ r2 ∼ 10−3M ∼
p p ∗
not overly concerning. Although our simulations are 2D, 300M ofgascouldbeevacuatedby∼6super-Earthsweigh-
no substantive difference between 2D and 3D treatments of ing a ⊕total of ∼30M . In such an initially gas-rich environ-
⊕
planet-diskinteractionsinviscousdiskshasbeenreportedvis-
ment, we anticipate the planets would migrate to and fro by
a`-visgapopening(Fung&Chiang2016)orplanetarytorques
a few tens of percent in orbital distance (Figure 1). As the
(Fungetal.2017). Ourneglectofdiskself-gravityshouldbe
disk drains in the long term, whatever slow and erratic mi-
an excellent approximation, as the Toomre Q-values of our
gration the planets undergo diminishes. Most of the angular
disks greatly exceed unity. Our planets are not allowed to
momentumofinterplanetarygaswouldbetransportedtothe
accrete gas, but super-Earths/sub-Neptunes are inferred ob-
outermostdisk,exteriortoalltheplanets,eitherbyLindblad
servationallytohaveonlymodestamountsofgas—lessthan
torquesorbydirectadvection.
10% by mass—acquired gradually over the entire disk life-
The large cavities of transitional disks may have been ex-
time(Lee&Chiang2016).
cavated over time by families of super-Earths. Observation-
More interesting frontiers to pursue include incorporat-
ally,gasdensitiesinsidecavitiescanbesuppressedrelativeto
ing disk thermodynamics, as radiative cooling and differen-
theirvaluesoutsidebytwotofourordersofmagnitude(e.g.,
tial heating across gap walls are thought to materially affect
Carmona et al. 2017; Dong et al. 2017). To reproduce these
planet-disk interactions (e.g., Kley & Nelson 2012; Tsang
strong depletions, appeal is commonly made to giant plan-
et al. 2014). Of course, extending the durations of the sim-
etsinviscousα-disksthatcanopendeepgaps(e.g.,Dong&
ulations,andincludingmoreplanetswithdifferentorbitalar-
Dawson2016;Dong&Fung2017). Butanalternativeinter-
chitectures, would also be welcome, for greater realism and
pretationisthatthecavitieshavebeenerodedgraduallyover
to enable more direct connections with observations. Closer
timebymuchsmallermassplanetsininvisciddisks.Sofaras
studyofgapdepthsiswarranted;weobservedsurfacedensity
wehavemeasuredinoursimulations,suchplanetsopengaps
contrastsonlyontheorderofunity(Figure2),incleardevi-
having only order-unity surface density contrasts. Neverthe-
ationfromscalingrelationsderivedfromviscousdisks(e.g., less,givensufficienttime,theycandraininterplanetarygasby
Fung et al. 2014), and surprising insofar as less viscous gas ordersofmagnitude.Becauseplanet-drivenaccretionratesM˙
should be less effective at diffusively back-filling gaps. The scale linearly with gas surface densities Σ, we have Σ˙ ∝ −Σ
shallowness of the gaps is due partly to the planets migrat-
whichimpliesexponentialdecayofthegascontent.Ifittakes
ingandre-startingthegap-openingprocessateachnewradial
t ∼ 105 yr to reduce a total disk mass of 10−3M by a
location (Malik et al. 2015). But how much of it is due to e−fold (cid:12)
factor of e (this assumes a contemporaneous mass transport
thelimiteddurationofoursimulations(≤ 5000yr),ortohy-
rate of 10−8M yr−1),2 then it takes 7t ∼ 7×105 yr to
drodynamicalinstabilitiesliketheRayleighinstability(Fung (cid:12) e−fold
reduceitbyafactorof1000. Anotherfeatureofthispicture
&Chiang2016),remainstobeworkedout. Finally,survival
isthatifsuper-Earthsdraindiskmassfasterthanittakesfor
against planetary migration is not guaranteed: planets with
theirnascentatmospherestocoolandacquiremoremass,they
mass M > M (Equation 1) stall but less massive planets
p cr maybeabletoforestallrunawayaccretion(Leeetal.2014).
donot. Thequestioniswhetherrockyplanetscancoagulate
Depleting the local disk density by a factor of (cid:38) 100 (rela-
fastenoughtocrosstheM thresholdbeforetheysuccumbto
cr tivetotheminimum-massextrasolarnebula)overtimescales
migration.
of∼1Myrsufficestokeepthegas-to-solidsmassfractionof
Our results support the proposal by Goodman & Rafikov
super-Earths(cid:46)10%,inaccordwithobservations(Lee&Chi-
(2001) that predominantly rocky planets—super-Earths and
ang2016).
Earths—can solve, or at least help to solve, the problem of
Onepotentialproblemwiththisscenarioisthatitpredicts
howprotoplanetarygasdisksultimatelydisperse.Giventhata
massaccretionratestolowerinproportiontodiskgasdensi-
singleplanetcanpushgasoveralengthscaleofapproximately
ties. Althoughsometransitionaldisksdohavelowaccretion
half its orbital radius, shuttling gas from 5 AU down to 0.1
rates (e.g., Dong et al. 2017), others do not, with a few hav-
AUwouldrequireabout6super-Earthsdistributedinroughly
ing M˙ as high as 10−7M yr−1 (e.g., Rosenfeld et al. 2014;
equal logarithmic intervals across this distance. Such planet (cid:12)
Carmona et al. 2017; Wang & Goodman 2017). Even so,
multiplicities are reasonable, given the profusion of super-
Earths/sub-Neptunes discovered by Kepler (e.g., Pu & Wu system-to-system variations in orbital architectures, particu-
larlyinthemassesoftheorbitingcompanions,couldhelpto
2015). Tobesure,adiskaccretionflowdrivenbyaplanetary
resolve this problem. We have focused here on super-Earths
systemwillbeunsteady, changingnotonlyonsecular, Myr-
becausetheyarecommonplace,butinprinciplegasgiantsand
longtimescales,butalsoonmuchshorterones,withmassal-
perhaps brown dwarfs or even low-mass stars may serve in
ternatelyaccumulatinganddispersingininterplanetaryspace,
theirstead.3
aswehaveseeninoursimulations(Figure2). Nomatterhow
complicated the history, however, all gas must ultimately be
torquedoutofsufficientlycompactplanetarysystems. Itmay 2Thistransportcanbeinwardoroutward—itshouldnotmatteraslongas
theregionoccupiedbytheplanetsismonotonicallydrainedofgasovertime.
betorquedbyplanetssofarinwardthatturbulencedrivenby
3 The case of the transitional disk HD 142527 is especially intriguing:
trhegeiomnasgwnehtioc-hroataretiosunffialciinesnttalbyiltihtye,rmacatlilvyatieodniiznedth(ee.ign.n,eDrmesocsht i2t0s0h6o)s,tanstdaritsacccarveitteyscaotntaairnasteao0f.2MM˙(cid:12)∼co1m0p−a7nMio(cid:12)nyhri−g1hly(GinarccliinaeLdotopetzheedtiaslk.
&Turner2015),takesoverthejobofdiskaccretionontothe (Casassusetal.2015;Lacouretal.2016).
8 Fung&Chiang
Thatdiskgasisintrinsicallyinviscid(laminar)issuggested recentdevelopments,inadditiontoourpresentwork,suggest
on other, independent grounds. Large-scale asymmetries in thatthereasonthecommunityhasnotdiscoveredarobustex-
transitional disks (e.g., Casassus et al. 2013; van der Marel planationforanon-zeroαinprotoplanetarydisksisthatnone
et al. 2013; Pinilla et al. 2015) have been interpreted as vor- exists: that in fact such disks are for the most part inviscid,
tices(e.g.,Zhu&Baruteau2016;Baruteau&Zhu2016),but and accrete primarily by the action of gravitational torques,
these vortices are spawned only in low-viscosity disks. Zhu exertedeitherbydiskgasitselfatearlytimes(e.g., Gammie
& Baruteau (2016) found that the Shakura-Sunyaev viscos- 2001;Cossinsetal.2009),4orbyplanetsatlatetimes.
ity parameter α needed to be 10−4 or lower before vortices
could grow from sharp density gradients. In other news, at-
We thank Pawel Artymowicz, Ruobing Dong, Eve Lee,
tempts to detect turbulence in the outer portions of disks us-
Hui Li, Zhi-Yun Li, Fre´de´ric Masset, Norm Murray, Ruth
ing molecular line observations have so far come up empty-
Murray-Clay,Sijme-JanPaardekooper,RomanRafikov,Yan-
handed (Flaherty et al. 2015; Flaherty et al. 2017, in prepa-
qin Wu, Zhaohuan Zhu, and an anonymous referee for en-
ration). AndRafikov(2017),inasystematicanalysisofdisk
couraging discussions and helpful feedback. This work was
accretion rates and masses, suggests that accretion may not
proceed viscously (i.e., diffusively), but may be enabled in- performedundercontractwiththeJetPropulsionLaboratory
stead by spiral density waves and/or disk winds. All these (JPL) funded by NASA through the Sagan Fellowship Pro-
gramexecutedbytheNASAExoplanetScienceInstitute. EC
isgratefulforfinancialsupportfromNASAandNSF.
REFERENCES
Alexander,R.,Pascucci,I.,Andrews,S.,Armitage,P.,&Cieza,L.2014, Hartmann,L.,D’Alessio,P.,Calvet,N.,&Muzerolle,J.2006,ApJ,648,484
ProtostarsandPlanetsVI,475 Hourigan,K.,&Ward,W.R.1984,Icar,60,29
Bai,X.-N.2016,ApJ,821,80 Ingleby,L.,Calvet,N.,Herczeg,G.,etal.2013,ApJ,767,112
Baruteau,C.,&Zhu,Z.2016,MNRAS,458,3927 Kley,W.,&Nelson,R.P.2012,ARA&A,50,211
Burke,C.J.,Christiansen,J.L.,Mullally,F.,etal.2015,ApJ,809,8 Lacour,S.,Biller,B.,Cheetham,A.,etal.2016,A&A,590,A90
Calvet,N.,D’Alessio,P.,Watson,D.M.,etal.2005,ApJL,630,L185 Lee,E.J.,&Chiang,E.2016,ApJ,817,90
Carmona,A.,Thi,W.F.,Kamp,I.,etal.2017,A&A,598,A118 Lee,E.J.,Chiang,E.,&Ormel,C.W.2014,ApJ,797,95
Casassus,S.,vanderPlas,G.,M,S.P.,etal.2013,Nature,493,191 Li,H.,Lubow,S.H.,Li,S.,&Lin,D.N.C.2009,ApJL,690,L52
Casassus,S.,Marino,S.,Pe´rez,S.,etal.2015,ApJ,811,92 Lin,M.-K.,&Papaloizou,J.C.B.2010,MNRAS,405,1473
Chiang,E.,&Laughlin,G.2013,MNRAS,431,3444 Malik,M.,Meru,F.,Mayer,L.,&Meyer,M.2015,ApJ,802,56
Christiansen,J.L.,Clarke,B.D.,Burke,C.J.,etal.2015,ApJ,810,95 Masset,F.S.,&Papaloizou,J.C.B.2003,ApJ,588,494
Clanton,C.,&Gaudi,B.S.2014,ApJ,791,91 Mu¨ller,T.W.A.,Kley,W.,&Meru,F.2012,A&A,541,A123
Colella,P.,&Woodward,P.R.1984,JournalofComputationalPhysics,54, Paardekooper,S.-J.,&Papaloizou,J.C.B.2009,MNRAS,394,2297
174 Peplinski,A.2008,PhDthesis,UniversityofStockholm
Cossins,P.,Lodato,G.,&Clarke,C.J.2009,MNRAS,393,1157 Pinilla,P.,vanderMarel,N.,Pe´rez,L.M.,etal.2015,A&A,584,A16
Cumming,A.,Butler,R.P.,Marcy,G.W.,etal.2008,PASP,120,531 Pu,B.,&Wu,Y.2015,ApJ,807,44
D’Angelo,G.,&Lubow,S.H.2010,ApJ,724,730 Rafikov,R.R.2002a,ApJ,569,997
Desch,S.J.,&Turner,N.J.2015,ApJ,811,156 —.2002b,ApJ,572,566
Dong,R.,&Dawson,R.2016,ApJ,825,77 —.2017,ApJ,837,163
Dong,R.,&Fung,J.2017,ApJ,835,146 Rosenfeld,K.A.,Chiang,E.,&Andrews,S.M.2014,ApJ,782,62
Dong,R.,vanderMarel,N.,Hashimoto,J.,etal.2017,ApJ,836,201 Sari,R.,&Goldreich,P.2004,ApJL,606,L77
Dressing,C.D.,&Charbonneau,D.2015,ApJ,807,45 Sicilia-Aguilar,A.,Henning,T.,&Hartmann,L.W.2010,ApJ,710,597
Flaherty,K.M.,Hughes,A.M.,Rosenfeld,K.A.,etal.2015,ApJ,813,99 Tsang,D.,Turner,N.J.,&Cumming,A.2014,ApJ,782,113
Fressin,F.,Torres,G.,Charbonneau,D.,etal.2013,ApJ,766,81 vanderMarel,N.,vanDishoeck,E.F.,Bruderer,S.,etal.2013,Science,
Fung,J.2015,PhDthesis,UniversityofToronto,Canada 340,1199
Fung,J.,&Chiang,E.2016,ApJ,832,105 Wang,L.,&Goodman,J.J.2017,ApJ,835,59
Fung,J.,Masset,F.,Lega,E.,&Velasco,D.2017,AJ,153,124 Ward,W.R.,&Hourigan,K.1989,ApJ,347,490
Fung,J.,Shi,J.-M.,&Chiang,E.2014,ApJ,782,88 Yu,C.,Li,H.,Li,S.,Lubow,S.H.,&Lin,D.N.C.2010,ApJ,712,198
Gammie,C.F.2001,ApJ,553,174 Zhu,Z.,&Baruteau,C.2016,MNRAS,458,3918
GarciaLopez,R.,Natta,A.,Testi,L.,&Habart,E.2006,A&A,459,837
Goodman,J.,&Rafikov,R.R.2001,ApJ,552,793
4Inself-gravitatingdisks,characterizingtransportintermsofanon-zero
forceandnotnaturallycapturedwithinalocaltheoryliketheonedefiningα.
αiscommonlydone,butmostlyforconvenience. Gravityisalong-range
