Particle Trapping and Streaming Instability in Vortices

12.7 MB·English

by  Natalie Raettig

DRAFTVERSIONJANUARY23,2015
PreprinttypesetusingLATEXstyleemulateapjv.5/2/11
PARTICLETRAPPINGANDSTREAMINGINSTABILITYINVORTICESINPROTOPLANETARYDISKS
NATALIERAETTIG1,HUBERTKLAHR1,ANDWLADIMIRLYRA2,3,4
DraftversionJanuary23,2015
ABSTRACT
We analyse the concentration of solid particles in vortices created and sustained by radial buoy-
ancyinprotoplanetarydisks,i.e. baroclinicvortexgrowth. Besidesthegasdragactingonparticles
wealsoallowforback-reactionfromdustontothegas. Thisbecomesimportantwhenthelocaldust-
to-gasratioapproachesunity.Inour2D,local,shearingsheetsimulationsweseehighconcentrations
5
ofgrainsinsidethevorticesforabroadrangeofStokesnumbers, St. Aninitialdust-to-gasratioof
1
1:100caneasilybereversedto100:1forSt=1. Theincreaseddust-to-gasratiotriggersthestreaming
0
instability,thuscounter-intuitivelylimitingthemaximalachievableoverdensities. Wefindthatpar-
2
ticletrappinginsidevorticesopensthepossibilityforgravity-assistedplanetesimalformationeven

n forsmallparticles(St=0.01)andlowinitialdust-to-gasratios(1:104).
a
J Subject headings: planets and satellites: formation, protoplanetary disks, accretion disks - hydrody-
namics-instabilities
2
2

1. INTRODUCTION most independently on distance from the star a max-
]
P Aquantitativeprescriptionofplanetesimalformation imum Stokes Number of about 0.01−0.1 which corre-
E isoneofthekeyissuesofplanetformationtheory.Mod- sponds to 1−10 cm at 1AU, 0.3−3 cm at 10 AU and
0.3−3 mm at 100 AU (see figs. 5 & 6 in Birnstiel et al.
. elsofsimplecollisionalstickingarecontroversiallydis-
h 2012).
cussed, both conceptually and also in the explanation
p To form large planetesimals these difficulties need to
ofcurrentpropertiesofasteroidandKuiperbeltobejcts
- becircumvented.Oneproposedmethodisgravitational
o (Weidenschilling2011;Morbidellietal.2009;Nesvorný
instability (Safronov 1972; Goldreich & Ward 1973; Jo-
r etal.2011;Bottkeetal.2005).Whethersmalldustgrains
t hansenetal.2006a,2007): Whenasufficientamountof
s stick to one another, bounce or fragment depends on
particles is close enough together, their mutual attrac-
a theirsizeandtheirrelativevelocities. Ingeneral,theory
[ predictsthatcollisionalvelocitiesriseasparticlesgrow, tion can trigger gravitational collapse, rapidly forming
which holds for particles with Stokes number smaller large planetesimals that then sweep up small particles
1 than1.TheStokesnumberisdefinedasSt≡Ωτ,where (Lambrechts & Johansen 2012). Different methods to
v f capturedusthavebeenstudied,suchaszonalflowsby
4 Ω is the Keplerian frequency, and τf is the stropping Johansenetal.(2011)andDittrichetal.(2013),pressure
6 time of a particle; both drift velocity and turbulence-
rings around stars (Whipple 1972; Klahr & Lin 2001,
3 induced relative velocities have a maximum for St=1,
2005), convection cells (Klahr & Henning 1997) or vor-
5 i.e., when the stopping time is of the order of an or-
ticeseithernumerically(Tangaetal.1996;Johansenetal.
0 bitalperiod(seee.g.reviewbyWeidenschilling&Cuzzi
2004; Lyra et al. 2008a, 2009a,b; Meheut et al. 2012a,b)
. 1993). At 5AU, St=1 corresponds to meter-size objects.
1 oranalytically(Barge&Sommeria1995;Chavanis2000;
0 Fragmentation occurs at velocities of only a few ms−1, Klahr&Bodenheimer2006;Chang&Oishi2010;Lyra&
5 which limits particle sizes to mm-cm (Wurm & Blum Lin2013).
1 2000; Brauer et al. 2008; Birnstiel et al. 2010). Güttler
From these studies it has become clear that dust can
: et al. (2010) and Zsom et al. (2010) introduced another
v easily concentrate in anti-cyclonic vortices. However,
boundary, the so called bouncing barrier where parti-
i feedbackofdustonthevorticalflow,whichcanpoten-
X cles hit one another and bounce without mass transfer.
tiallyfrustrateparticleconcentration,hasnotbeenstud-
At even smaller size scales, Okuzumi et al. (2011a,b)
r iedindetail. Thisisthetopicofthepresentstudy.
a foundthechargebarrier,wheresmallparticlesarepre-
Besides Lyra et al. (2008a, 2009a,b); Meheut et al.
vented from approaching one another due to the elec-
(2012b); Ataiee et al. (2013); Zhu & Stone (2014), and
tric charges built up through collisions with free elec-
Zhuetal.(2014),whoanalyzedparticletrappinginvor-
trons. Birnstiel et al. (2012) determined the sizes that
ticesexcitedbytheRossbywaveinstability,otherworks
particlescanobtainlocallyastheirgrowthislimitedby
did not choose a particular vortex formation mecha-
radial drift and collisional destructions. They find al-
nism. In our study, vortices are naturally produced by
the radial stratification of disks (Klahr & Bodenheimer
1Max-Planck-InstitutfürAstronomie,Königstuhl17,69117,Hei- 2003; Petersen et al. 2007a,b; Lesur & Papaloizou 2010;
delberg,Germany;[email protected],[email protected]
Lyra & Klahr 2011; Raettig et al. 2013). The growth of
2Jet Propulsion Laboratory, California Institute of Technol-
the vortices occurs proportionally to the radial buoy-
ogy, 4800 Oak Grove Drive, Pasadena, CA 91109, USA;
[email protected],[email protected] ancy frequency (aka Brunt-Väisälä frequency) squared
3DivisionofGeological&PlanetaryScience,CaliforniaInstitute N2 (Raettigetal.2013)andthelocalthermalrelaxation
of Technology, 1200 E California Blvd MC 150-12, Pasadena CA, timescaleasarguedbyLesur&Papaloizou(2010),nu-
91125USA
4SaganFellow mericallyconfirmedbyRaettigetal.(2013). Thebuoy-
2 Raettigetal.
ancyfrequency,whichisafunctionoftherelativescale bysolvingtheequationsofmotionsforLagrangianpar-
height H/r, logarithmic radial entropy and pressure ticles(seeSection3).
gradientsβ ,β andtheadiabaticindexγ: In general the gas and dust feel the same external
K P
forcesexceptthepressureforce f =−ρ−1∇p, where ρ
N2=−1 (cid:18)H(cid:19)2β β . (1) is the gas density and p the prespsure. gThis term onlyg
K P
γ r affects the gas. For instance the global pressure gradi-
entinthegasleadstoasub-Keplerianorbitalgasveloc-
A linear theory to explain this behaviour has been put ityu.Thecorrespondingbuoyancyforceforaparticleis
forward by Klahr & Hubbard (2014) who identify the f = −ρ−1∇p, where ρ is the material density of the
instabilityasaradial"convectiveoverstability"inaccre- p,s s s
solidmaterial,canbeneglectedbecauseρ (cid:29)ρ . Since
tiondisks. Lyra(2014)performed3Dsimulationsofthis s g
the particlesdo not feelthis global pressure, they need
"convective overstability" and showed that in the non-
toorbitatKeplerianvelocitiesinordertobeincentrifu-
linear phase indeed vortices were emerging from the
gal balance with stellar gravity. The resulting veloc-
flow.
ity difference between gas and large particles acts as a
In this paper we set out to assess how efficiently
headwind on the particles which decreases their angu-
dust can be concentrated in vortices enforced by real-
lar momentum leading to radial drift inward. Smaller
istic values for the radial stratification in temperature
particlesgetdraggedalongbythegasandthereforfeela
and density, also with a plausible value for the ther-
netaccelerationtowardsthestaralsomakingthemdrift
mal relaxation time. We do this as a first step via two-
inward.
dimensional, local simulations. Ultimately, only three-
Thetimeonwhichthedustparticlesadjusttothegas
dimensional stratified runs will be able to include all
velocity is the friction time. For subsonic Epstein drag
relevantphysicsfromvortexstability(Barranco&Mar-
cus2005;Lesur&Papaloizou2009;Lyra&Klahr2011) this is τs = ρss(ρgcs)−1 (Weidenschilling 1977) which
to sedimentation of dust. Nevertheless, here we find depends on particle size s and local sound speed cs.
even in 2D that particle concentration in the vortices TheEpsteinregimeconsidersparticlessmallerthanthe
issufficientlyintensethatback-reactionoftheparticles mean free path of gas. Larger particles have to be
ontothegasmotionhastobeconsidered. Wefindthat treatedintheStokesregime. Butaslongasrelativeve-
thestreaminginstability(Youdin&Goodman2005;Jo- locities between dust and gas are small, it is possible
hansenetal.2006a)istriggered,counter-intuitivelylim- to describe the coupling force to be linear with respect
itingthemaximumdust-to-gasratio, andseverelyper- totherelativevelocity f =mδu independentonthede-
turbing the gas flow inside the vortex. In some cases τs
tailedcalculationof τ. Especiallyitisnotnecessaryto
s
we see vortices getting disrupted and later reforming,
determine the shape and density of the particles or the
startinganewcycleofparticleconcentration.
gasdensityandsoundspeeduntiloneasksforthepar-
In simulations similar to ours but three-dimensional
ticlesizecorrespondingtoagivenfrictiontime.
and more numerically expensive, where particles were
This particle gas coupling is usually expressed in
trappedinzonalflowsofmagneto-hydrodynamicalori- terms of the dimensionless Stokes number St ≡ Ωτ.
s
ginandincludingself-gravity,theoverdensitiesthatwe
Particles of different size, shape and density, but ulti-
foundinthepresentpaperdidalreadyleadtogravita-
mately with the same Stwill behave the same as far as
tional collapse (Johansen et al. 2006b, 2007) and to the
aerodynamicsareconcerned.
formationofplanetesimals. Ontheotherhand2Dsim-
The friction time and thus the Stokes number is a
ulationsincludingparticlefeedbackareideallysuitedto
function of the local gas density. We nevertheless as-
studyawideparameterrangetolearnwhetherstream-
sume the Stokes number to be constant, which is justi-
inginstabilitycanbetriggeredandhowitaffectsvortex
fiedasthegasdensityfluctuatesonlybyafewpercent
stability.
in the sub-sonic turbulence we consider in our simula-
Thepaperisstructuredasfollows.Wefirstreviewthe
tions.
underlyingphysicsofdust-motioninagasdisk,specif-
Since our simulations are only two-dimensional we
icallyinavortexincludingdust-gasinteractions. Then
can not consider vertical settling of particles. The fact
wedescribethenumericalsetup.Theresultsofoursim-
thatthereal3Dgasdiskisverticallystratifiedandthus
ulationsaregiveninSection4. Hereweespeciallylook
particlesofagivensizehaveaverticallyvaryingfriction
at the reached dust-to-gas ratios compared to the aver-
time,hastobediscussedinmoredetail. Onecaninter-
ageinitialdust-to-gasratio. InSection5wediscussthe
pretthisapproximationasfocussingonanarrowverti-
prospects to form planetesimals via gravitational frag-
calrangearoundthemidplane,inwhichneitherthegas
mentation of the dust enhancements in vortices and in
nor the particle stratification are too strong. Then one
Section 6 we present our measurements of collisional
can assume that also the friction time and the dust to
speeds among dust grains. Finally, we summarize and
gasratiohasnostrongverticalvariations. Thisassump-
concludeinSection7.
tionisalittlelessstringentthansayingthatallparticles
areactuallyinthemidplane.
2. PHYSICALBACKGROUND Thegeneralbehaviorofparticledriftinanon-laminar
Theevolutionofthegascomponentofthediskinthe yet geostrophic flow5 can be understood easily from a
PencilCodeisgivenbytheNavierStokesEquationcon- a simplified treatment of the Lagrangian equations of
taining stellar gravity as well as the virtual forces of
the rotating and shearing box as explicit terms (Lyra & 5 In a geostrophic flow the pressure gradients are in equilibrium
Klahr2011). Dustgrainsontheotherhandareevolved withthecentrifugalandCoriolisacceleration.
ParticleTrappingandStreamingInstabilityinVortices 3
motionforgasanddust:
TABLE1
SIMULATIONSETUPANDεmax
run St Σd,0/Σg,0 feedback εmax FractionofMd FractionofMd
1 ρ withε≥0.01(%) withε≥1(%)
∂ u=F− ∇p− d (u−v) (2) NF1 0.01 1:100 no 6.26 86.64 2.88
t
ρ ρ τ
g g s NF2 0.05 1:100 no 636.00 92.31 75.63
1 NF3 1 1:100 no 843.70 79.71 1.26
∂tv=F− (v−u), (3) NF4 20 1:100 no 46.81 99.90 99.80
τ
s F1 0.01 1:100 yes 1.07 70.40 0
F2 0.05 1:100 yes 3.86 95.24 2.13
F3 1 1:100 yes 77.33 98.73 83.74
F4 20 1:100 yes 0.73 79.82 8.09
wherevistheparticlevelocityand F collectstheterms DG1 0.05 1:1000 yes 1.15 86.74 0
that are the same for both dust and gas, e.g., gravita- DG2 0.05 1:10000 yes 0.70 78.30 0
tional force from the central star. The second term of DG3 1 1:1000 yes 11.53 98.46 56.37
theequationsdescribethepressureforce,andthethird DG4 1 1:10000 yes 4.17 97.30 76.15
termisthefrictionforcebetweengasanddustparticles.
WithsomemathitcanbeshownthatbysubtractingEq.
(2) and (3) and assuming |∂ u−∂ v| (cid:28) |u−v|/τ par-
t t s TABLE2
ticles will move towards regions with higher pressure COLLISIONALVELOCITIES
v=u+τ ∇pevenifthepressuremaximumistinyand
s
thTehperovfiolertriecleastiwveelycsomnsoiodtehr(saereeea.gn.tiK-clyachlro&nicLivno2r0ti0c1e)s. St ε0 α vccosll √urαmcss √vrαmcss √vcαoclls ∆√vαcocsll vc√olαl,cOsC
and have a slightly lower epicyclic frequency than the (/10−2) (/10−3) (/10−2) (/10−2)
Keplerianorbitalfrequency. Sincethereisageostrophic 1 10−2 1.03 2.77 2.84 1.37 2.73 2.56 1
balance between the Coriolis and the pressure force a 1 10−3 1.56 3.17 4.36 3.61 2.54 3.34 1
highpressureregioninsidethevortexiscreated. Parti- 1 10−4 1.54 1.25 3.88 2.27 1.01 4.01 1
0.05 10−2 1.56 2.86 3.27 2.41 2.29 0.59 0.31
cle accumulation inside a vortex then basically works
0.05 10−3 1.57 2.71 4.29 3.80 2.16 1.28 0.31
the same way as the radial drift in an accretion disk
0.05 10−4 1.54 2.08 3.88 3.04 1.68 0.09 0.31
works. The particles do not feel the pressure gradient
insidethevortexandthereforetheirepicyclicfrequency
equals the Keplerian frequency. Yet, the rotation fre-
quency of a pressure supported vortex is smaller than dust-to-gasratioε=Σ /Σ
d g
the epicyclic frequency. Thus the headwind from the
gas causes the particles to lose eccentricity and forces DΣg +(u·∇)Σ =−Σ ∇·u+ f (Σ ), (4)
the particles to spiral towards the center of the vortex. Dt g g D g
Foranindepthanalysisofvortexcapturingincludinga Du 1
comparison of all forces acting on the gas and dust in- +(u·∇)u=− ∇p−2Ω (zˆ ×u)
Dt Σ 0
side the vortex see Adams & Watkins (1995); Barge & g
(cid:18) (cid:19)
Sommeria (1995); Tanga et al. (1996); Chavanis (2000) 3 βp 1 1
andJohansenetal.(2004). +2Ω0uxyˆ + R0 Σ − Σ xˆ + fν(u,Σ)
0 g 0
If the dust density becomes comparable to the gas
Σ
density the drag forces from particles onto the gas can − d (u−v), (5)
nolongerbeneglectedaslongastheStissmallerthan Σ τ
g s
≈ 10. This back-reaction can alter the motion of the Ds 1 (cid:26) (T−T )
gas and also leads to even higher dust concentrations +(u·∇)s= ∇·(K∇T)−Σ c 0
through the streaming instability (Youdin & Goodman Dt ΣgT g v τcool
2005). Thelasttermofequation(2)representstheback- (cid:27)
βp u
reactionfromdustgrainsontothegas. + 0 x + f (s). (6)
r (γ−1) K
0
Here, Σ is the gas density, u is the deviation of the
g
gas velocity from the Keplerian value, s the entropy, T
thetemperature,c thespecificheatatconstantvolume,
v
τ isthethermalrelaxationtimescale,andKtheheat
cool
3. NUMERICALSETUP conductivity. Thesymbol
We perform two dimensional shearing sheet simula-
D ∂ ∂
tions with the PENCIL CODE where we consider the ≡ +u(0) (7)
vertically integrated densities Σ instead of the three- Dt ∂t y ∂y
dimensionaldensitiesρ.TheEuler-equationsforthegas
are solved on a Cartesian grid, identical to the setup represents the Keplerian derivative where u(y0) =
of (Raettig et al. 2013), but now augmented by a term −3/2Ω x and Ω is the Keplerian frequency at the or-
0 0
fortheparticlefeedbackonthegasε(u−v)/τ withthe bit of the shearing box R . The radial deviation from
s 0
4 Raettigetal.
R is given by x and the linearized azimuthal direc- positionofthedustparticle.Theindex∗isomittedfrom
0
tion is now measured in terms of y. Boundary condi- nowon.
tions are periodic in y and shear-periodic in x. Fur- For our 2D simulations we choose η = −0.01 which
ther terms in the equations are diffusion terms to en- corresponds to a pressure gradient of β=2 and a disk
sure numerical stability of the finite difference Pencil aspect ratio of h = 0.1. The terminal velocities of the
Code f (Σ ),f (u,Σ),f (s) and radial stratification of particles for a given Stokes number in a laminar disk
D g ν K
the disk in terms of pressure and entropy for radially arenowintermsofη(Weidenschilling1977)
constant density: β = −dlogp = −dlogs. The stratifica- 2ηu
dlogR dlogR v = K , (12)
tiontermβoccursinthelinearizedcomponentofradial x St+St−1
buoyancy (Eq. (5)) and in the term for radial transport ηu
of entropy (Eq. (6)) as derived in Lyra & Klahr (2011). vy= St2K. (13)
Thesetermsincombinationwiththermalrelaxationare
thedriverofvorticesasfoundinPetersenetal.(2007a). Thesearethebestassumptionsonecandofortheinitial
As described in their paper gas that is moving radially velocitiesofparticlesinoursimulations.
outwardisgenerallywarmerandoflowerdensitythan Ourphysicaldomainspans4diskscale-heights,±2H,
gas moving inward. The density mismatch across the around the co-rotational radius in radial-direction (x-
vortexleadstoamismatchinbuoyancyandthusavor- axis)and16H inazimuthal-direction(y-axis). Thegrid
texwillfeelatorqueacceleratingit.Themaximalampli- itself consists of 2882 grid cells6. We choose β =2 for
fication occurs when thermal relaxation is on the same theentropyandpressuregradientswhichisarelatively
orderastheinternalrotationperiodofthevortex. stronggradientcomparedtogradientsweexpectinpro-
Thedustgrainsaremodeledwithaparticleapproach. toplanetary accretion disks are between β = 0.5 and
For each individual particle we solve the equation of β = 1. However, Raettig et al. (2013) found that the
motion including gas drag. We do not allow for self- general behavior of vortices is the same for weak and
gravitysofar. Inprinciple, foranindividualparticleit strongentropygradients. Thedevelopmentofthevor-
doesnotmatterifthereareotherparticlesinthesimula- ticesismerelyfasterwithstrongergradientsandwecan
tionornot.However,inthesimulationswhereweallow scanalargeparameterspacewithareasonableamount
forback-reactionoftheparticlesonthegas, thereisan ofcomputingtime. Forafirstestimationthestronggra-
indirect influence of one particle on the other particles dientissufficient,butinourfuture3Dstudiesincluding
duetothealteredgasvelocity. self-gravity we will use β=1 because then we aim for
Thegasdiskisco-rotatingwithKeplerianvelocityat quantitativeresultsonplanetesimalformationratesand
the co-rotational radius r0. But for simulations involv- massspectrum.
ingbothdustandgasweneedtoincludeavelocityoff- We first evolve the disk for 200 local orbits without
setduetothepressuregradient. Thepressuregradient particles. This way we make sure that a fully grown,
isbalancedbytheCoriolisforcewhichisthelinearized long lived vortex has developed before we put in par-
expression for the radial centrifugal acceleration in the ticles. Physically,thiscorrespondstovariousscenarios:
localco-rotatingsystem A)avortexdevelopingintheouterpartsofthediskand
thenmigratinginwardsthroughregionswithparticles,
1 ∂p H∂lnp
2Ωu = =c Ω . (8) B)thegrowthofparticlesintoasizeregimewherethey
y Σg ∂r s r ∂lnr get trapped by the nearest vortex or C) a radial flux of
particlesclosetobouncingand/ordriftbarrierthaten-
The deviation by which the gas velocity is lower than
tersthevicinityofthevortex.
theKeplerianvelocityu =Ωris(Nakagawaetal.1986)
K Aftertheinitial200orbits,wedistribute400000parti-
(cid:18)H(cid:19)2∂lnp clesrandomlyinthedisk. Thiscorrespondsto4-5par-
η=0.5 , (9) ticles per grid cell initially. That way we minimize nu-
r ∂lnr mericaleffectsthatcanariseiftherearenotenoughpar-
ticlesinthecomputationaldomain,e.g. theeffectasin-
whichleadstoasub-Kepleriangasvelocityandadevi-
gle particle has can be extremely overestimated if not
ationfromtheKeplerianvelocityu of
K enoughparticlesareconsidered. Wheneverwereferto
(cid:18)H(cid:19)2∂lnp timesinthispaper,wemeantimeelapsedsinceparticles
u =0.5 Ωr=ηu . (10) wereputintothesimulation. Atthetimewhenthedust
y r ∂lnr K
to gas ratio approaches unity and streaming instability
is triggered, 500 or more particles are in one grid cell.
Thisvelocitydeviationηisaddedtothesimulationsar-
FollowingJohansen&Youdin(2007)thesenumbersare
tificially.
safely beyond the minimum particle number per grid
For each particle i the equation of motion is solved
cellneededtoachievenumericalconvergence.
individuallyvia
Starting our simulations with particles included
would havelittle changedour resultsas theinitial and
dv∗(i) =2Ωv∗(iy)xˆ− 1Ωv(∗ix)yˆ − 1 (cid:16)v(∗i)−u(cid:16)x(i)(cid:17)(cid:17), (11) athveerdaugsetdtoushtatvoegiamspraactitoosnartheemguacshdtyonoalmowicsinriogrhdtefrrofomr
dt 2 τ
s
wherev∗(i)=v−ηuKisth(cid:16)epar(cid:17)ticlevelocitycorrectedby com6 Tphuetantiuomnablecrluosfte2r88wghriicdhcneellesdcsoammesulftriopmletohfe1a2rtcohipteacrtaullreeliozfeothuer
thevelocityoffsetand u x(i) itthegasvelocityatthe PENCILCODEefficiently.
ParticleTrappingandStreamingInstabilityinVortices 5
0.2
0.1
cs 0.0
> / y
uy −0.1
<
−0.2
−0.3
8

FIG.1.—Maximumdust-to-gasratio,εmax,forsimulationswithout 6
particlefeedback. Thedifferentlinesrepresentthedifferentparticle >
sizes: dotted(black)line: St=20, dash-dotted(green)line: St=1, Σ0
solid(blue)line: St=0.05,anddashed(red)line: St=0.01. Almost <
allparticlesofintermediatesize(St=1andSt=0.05)accumulatein > / y 4
thevortex.Stronglycoupledparticles(St=0.01)accumulateonlypar- Σp
tially,becausetheyalsocoupletothegasoutsideofthevortex.Large <
particles(St=20)hardlycoupletothegasatallandthereforearenot 2
affectedbythevorticalmotion.
0
the beginning. Retrospectively this assumption was
−0.2 −0.1 0.0 0.1 0.2
later justified in simulations (see Section 4.2) in which
x
the vortex gets disrupted via the streaming instability
and then forms again, this time in the presence of dust FIG.2.—Azimuthallyaveragedgasvelocityandparticledensityfor
grains. This second growth phase does not differ from St=20particles.Theparticlesaccumulateinthezonalflowstructure.
theinitialgrowthphasewithoutdustgrainspresent.
Each of the 400000 particles represents one super-
feedbackisincluded.
particle,acollectionofseveralparticles,allofthesame
Figure 1 shows the maximum dust-to-gas ratio for
size, and a given mass according to the initial dust-to-
these simulations. We see from the almost constant
gasratio.
dash-dottedgreenlinethatalloftheavailablesolidma-
Generally we have an initial dust-to-gas ratio ε of
0 terial in the St=1 particle simulations accumulates in
1:100whichmeansthatthediskconsistsoutof1%solid
thevortex. Thiscorrespondstoaparticleconcentration
materialand99%gas. Notethatintwodimensionsthe for St=1 particles of c≈Σ /Σ ≈12000 and ε≈320.
PENCIL CODEassumessurfacedensitiesinsteadofvol- St=0.05eventuallyalsoaldlaccdu,0mulatewithinthevor-
ume densities. Therefore the dust-to-gas ratios we talk
about in this paper refer to ε = Σ /Σ rather than to texandthenarekepttrappedthere. Themaximumpar-
d g ticleconcentrationisreachedsoonerforSt=1particles
ε = ρd/ρg. To simulate different particle sizes we use than for St=0.05 particles, which can be attributed to
different Stokes numbers of St=0.01, St=0.05, St=1 theirhigherradialdriftvelocity. EvenforSt=0.01par-
and St=20. By only defining the Stokes Number we ticlesthedust-to-gasratioincreasestoε≈1.
are able to cover particles both in the Epstein as well We see two different behaviors. For St = 0.05 and
asintheStokesregimeasfrictionforcesarelinearwith St=1,basicallyallparticlesarecollectedwithinthevor-
respect to velocity. At 5 AU this corresponds to par- tex,althoughSt=0.05take10timeslonger,yetinboth
ticles between 3mm and 6m for Σ = 300gcm−2 and casesconcentrationsofabout4ordersofmagnitudeare
g
ρ =2gcm−3 . The parameters for all our simulations reached. ForSt=0.01concentrationinsidethevortexis
s
canbeseeninTable1. Aseparatesimulationiscarried onlyby2ordersofmagnitudeandSt=20particlescon-
outforeachSt. centrate by one order of magnitude, but interestingly
Particlesgrowthduetocollisionsisnottreatedinour outside the vortex in correlation with the zonal flow
simulation. Physically particles would have a growth related to the radial position of the vortex (see Fig. 2).
timeofseveralthousandyears,thereforeweneglectthis Thesezonalflowsshowupintheazimuthallyaveraged
effectforthesakeofkeepingthesimulationsimpleand rotation profile, both in magneto hydrodynamical sim-
easytobeevaluated. ulations(Lyraetal.2008b;Johansenetal.2009;Dittrich
etal.2013)andalsoinourhydrodynamicalsimulations
4. RESULTS ofvortices.Thezonalflowsaredeviationsfromthesub-
Keplerian mean rotation profile confined by long-lived
4.1. NoParticleFeedback
radialpressurevariations(geostrophicbalance).
Wefirstconsidersimulationswherewedonotinclude The different behavior of particle accumulation can
back-reactionfromtheparticles,effectivelysettingρ = be explained through the gas coupling. Particles with
d
0 in Equation (5). That means that gas drags particles St=0.01 are strongly coupled to the gas. As an effect
along,butevenhighdustdensitiesdonotaffectthegas it takes much longer to capture the particles into the
velocity in any way. The purpose of these simulations vortexandasasecondeffectthedeviationbetweenthe
is to serve as comparison to the later models in which actual vortex and a analytical solution like for instance
6 Raettigetal.
FIG.3.—Eachofthefourpanelsshowvorticityωz (inunitsofthelocalKeplerianfrequencyΩ0),anddust-to-gasratioε. Theleft-hand-side
panelsdonotincludeparticlefeedback,whereastheonesintheright-hand-sidedo.TheparticlesizescorrepondtoStokesnumbersSt=20(top)
andSt=1(bottom).
ParticleTrappingandStreamingInstabilityinVortices 7
FIG.4.—SameasFig.3,butforparticlesofSt=0.05(top)andSt=0.01(bottom).
8 Raettigetal.
InthelastSectionweshowedthattherearesignificant
particleoverdensitieswithinthevortex.Thismeansthat
we cannot neglect the momentum transfer from dust
onto gas. Because of this we now include dust back-
reactiononthegas,usingthefullEquation(5).
For the St=20 particles there is no change with re-
spect to the simulations without particle feedback be-
cause overdensities stay below ε=1 (compare the left
andrightsideofthetoppanelinFig.3). Fortheselarge
particles the streaming instability can not be triggered
at the present dust-to-gas ratio. All other tested parti-
clesstillaccumulateinthevortex,yetdifferently:
Particles with St = 1 concentrate more locally (in a
smaller area) than smaller particles (bottom panel of
Fig.3andFig.4).St=1particlesaregettingconcentrated
on time scales which are not longer than the dynami-
cal time scale of the vortex (Barge & Sommeria 1995),
thus they can easily follow the changes of the vortex
in shape, strength and location. Smaller particles have
much longer time scales to spiral into the center of the
vortex-longerthanthedynamicaltimescaleofthevor-
tex. Astheattractorinsidethevortexischangingitslo-
cation on dynamical time scales (change of amplitude,
shapeandlocationofthevortex), theparticleshaveno
chancetoevercatchuporfollowit.
An effect that occurs when ε approaches unity is the
streaminginstability(Youdin&Goodman2005;Youdin
& Johansen 2007). Once particles concentrate, their lo-
callyincreaseddust-to-gasratioleadstoaslowerradial
drift(Nakagawaetal.1986).Thisproducesafurtheren-
hancement of solids since faster particles from slightly
larger radii bump into the accumulation, like a traffic
jam, which eventually results in streaming dust struc-
tures. Thisisalsothecaseinoursimulations,although
withourresolutionitispossibletoresolveonlyapartof
theunstablewavelengths. Thereforewearenotableto
FIG.5.—Gasvorticity(leftpanel)anddust-to-gasratio(rightpanel) studythelinearandnonlinearbehaviorofthestreaming
forasimulationwithSt=0.01particlesandback-reactionontothe instability in all detail and model the correct growth-
gas. Theellipticalvorticalgasflowisdistinguishableinthevorticity
rates of the instability. Yet, for the wavelength we re-
plot. Therearestrongaccumulationsofparticleswithinthevortex.
Althoughmanyparticlesspreadoutovertheentirevortex,mostpar- solve we find the right threshold dust-to-gas ratio to
ticlesconcentrateinthecenterofthevortex. Thepositivevorticity trigger the instability. See also the appendix for addi-
valuesinthevortex(lightareasintheleftplot)showtheeffectthat tionalresolutiontestsonthestreaminginstability.
particleshaveonthegas.Wherethegasencountersobstacles,namely
Figure 5 shows the vertical gas vorticity ω in units
highparticleconcentrations,steepvorticitygradientsdevelop. z
ofthelocalKeplerianfrequencyΩ anddust-to-gasra-
0
tio ε=Σ /Σ for a simulation with St=0.01 particles
g d
a Kida vortex (see discussion in Lyra & Lin 2013) pre- after 300 local orbits (run F1). It is clear that the parti-
vents the particles from accumulating too densely. An clesaccumulateinsideofthevortexandfollowthevor-
additionaleffecttopreventtoostrongconcentrationsin tical motion. Where the concentration is highest they
3D vortices might be turbulence, triggered for instance createstrongmaximainthegasvorticity. Thisisanef-
by elliptical instability (Lyra & Lin 2013) and internal fectofthecyclonicrotationoftheparticleclumpswhich
circulation (Meheut et al. 2012a). The St=20 particles was already reported in earlier work (Lambrechts &
are the other extreme: they are hardly coupled to the Johansen 2012) and is a consequence of the conserva-
gasatall,andthereforeareonlyweaklyaffectedbythe tionofangularmomentumundercontractionwhichcan
vorticalmotionofthegas. Inastationaryanalyticvor- turnanticyclonicmotioninarotatingsystemtocyclonic
tex these particles would get trapped if one waits long motion.
enough. Yet, our numerical vortex is dynamically ac- Ingeneralduetotheback-reactionofthedustonthe
tiveinchangingitsstrengthandshapeandadditionally gas and the resulting streaming instability the initially
drifting in the radial direction. This vortex dynamics elliptical gas streamlines are bent into more complex
gives the large particles the chance to escape from the motionsthanincaseswithoutback-reaction.
vortex again, but still stopped from radial drift in the In case of the St = 1 particles where the local parti-
zonalflow. cle concentration, and therefore the back-reaction, are
strongest, the vortex structure is disrupted. Because
4.2. IncludingParticleFeedback ofthis, theparticletrappingmechanismslosestrength:
ParticleTrappingandStreamingInstabilityinVortices 9
FIG.6.—Vorticityωz inunitsofthelocalKeplerianfrequencyΩ0(1stand3rdcolumn)anddust-to-gasratio(2ndand4thcolumn)forSt=1
particlesfordifferentpointsintime.Atthefirstsnapshot,thevortexisstillclearlydistinguishable.Inthesecondsnapshot,ithasbeendisrupted
stronglybytheparticleaccumulation.Asthisparticleaccumulationspreadsout,thevortexcanslowlyregainitsshape(3rdsnapshot)andforma
large,yetstillperturbed,vortexagain(4thsnapshot).
10 Raettigetal.
FIG.7.—Maximumdust-to-gasratio,εmax,forsimulationswithpar-
ticlefeedback.Thedifferentlinesrepresentthedifferentparticlesizes:
dotted (black) line: St=20, dash-dotted (green) line: St=1, solid
(blue)line: St=0.05, anddashed(red)line: St=0.01. Particlesof
St=1reachthehighestdustenhancementsandconcentrateinavery
localareainthevortex, whilesmallerparticlesspreadoutoverthe
entirevortex.Thereforetheoverdensitiesreachedarelower.
100
Md 10−1
of gFasIGr.a1ti0o.s—foPraSrtti=cle1c(otonpce)natnradtiSotn=(s0o.l0id5l(ibnoetst)oamn)dpmaratxicilmesu.mThdeuscto-ltoor-
n 10−2 representthedifferentinitialdust-to-gasratios:1:100(black),1:1000
o
acti St=20 (turerded),ainndth1e:v1o0r0ti0c0es(gforerean)l.owMoinreitiianldεi,vbideucaaulsseuptheeri-rpabratcickl-erseaacrteiocnapis-
fr 10−3 St=1 lessefficient. Sinceeachofthesesuper-particlesislessmassivethan
St=0.05 withhigher ε, theoveralldust-to-gasratioforlowinitial ε islower
10−4 St=0.01 thanthatoflargerinitialε.
10−4 10−3 10−2 10−1 100 101 102
Σ/Σ
d g Fig.6oftheelectronicversionofthispaper).
The destruction of the vortex occurs on very short
FIG.8.—Fractionoftheentiredustmassinaspecificdust-to-gasra- time scales set by the streaming instability, whose
tioorhigher.ε≥1arereachedforallparticlesizes.ForSt=1particles
growth time is a function of the dust to gas ratio and
morethan80%ofallparticleshaveconcentratedinareaswithlargerε
than1. Stokesnumberofthedustthatcanbeasfastasoneor-
bitalperiod(Johansen&Youdin2007). Thenewvortex
amplification on the other hand takes hundreds of or-
1.0 bitslikeintheinitialgrowthasitdependsontheradial
stratificationofthediskandthethermalrelaxationtime
0.8 ofgas.
Md This vortex disruption and reforming process dis-
of 0.6 cussed here is not to be confused with the vortex in-
n stability discussed in Chang & Oishi (2010). They an-
o
cti 0.4 feedback alyze the stability of a three dimensional vortex based
a
r no feedback on the density contrast between the interior of a vor-
f 0.2
texandtheambientmediumanddeterminethatifthis
0.0 contrastishigherthanafew10%thenthevortexwillbe-
10−4 10−3 10−2 10−1 100 101 102 comeunstable. Althoughweseesuchhighdustdensity
Σ/Σ concentrationsinsidethevortextheyareverylocalized
d g andthereforeourvorticesremainstable.Forthemtobe-
come entirely unstable according to Chang and Oishi’s
FIG.9.— Fractionoftheentiredustmassinaspecificdust-to-gas
ratioorhigher. ε≥1arereachedforSt=0.05particleswith(dashed analysisthedensitycontrastneedstobeuniformlyhigh,
blackline)andwithoutfeedback(solidblueline). spread out over the entire vortex, which is never the
caseinoursimulations.
pressuregradientsacrossthevortexbecomeshallower, Figure7showsthemaximumdust-to-gasratioofour
Coriolisforcesinthevortexgrowweaker,andtheparti- different simulations with different St. We clearly see
clesbegintoescapethevortex. Becausethelocalparti- that St = 1 particles (dash-dotted green line) have the
cleconcentrationdecreases,thelargevorticitygradients highestconcentration.Astheparticlesizedecreases,the
flattenoutagain.Thiseventuallyleadstoanewamplifi- concentrationalsodecreases. Itisimportanttonotethat
cationofthevortexduetothebackgroundstratification St = 20 particles do not accumulate inside the vortex.
andtheprocessrepeatsitself(seeseriesofsnapshotsin Thehighestdust-togasratioisreachedforSt=1parti-