Gravoturbulent Planetesimal Formation: The Positive Effect of long-lived Zonal Flows

3.8 MB·English

by  Karsten Dittrich

Draftversion January10,2013
PreprinttypesetusingLATEXstyleemulateapjv.5/2/11
GRAVOTURBULENT PLANETESIMAL FORMATION: THE POSITIVE EFFECT OF LONG-LIVED ZONAL
FLOWS
K. Dittrich1, H. Klahr1, and A. Johansen2
1Max-Planck-Institutfu¨rAstronomie,K¨onigstuhl17,D-69117Heidelberg,Germany and
2LundObservatory,DepartmentofAstronomyandTheoreticalPhysics,Box43,SE-22100Lund,Sweden
Draft version January 10, 2013
ABSTRACT
3 Recent numerical simulations have shown long-lived axisymmetric sub- and super-Keplerian flows
1 in protoplanetary disks. These zonal flows are found in local as well as global simulations of disks
0
unstabletothemagnetorotationalinstability. Thispapercoversourstudyofthestrengthandlifetime
2
of zonal flows and the resulting long-lived gas over- and underdensities as functions of the azimuthal

n and radial size of the local shearing box. We further investigate dust particle concentrations without
a feedback on the gas and without self-gravity. Strength and lifetime of zonal flows increases with
J the radial extent of the simulation box, but decreases with the azimuthal box size. Our simulations

9 support earlier results that zonal flows have a natural radial length scale of 5 to 7 gas pressure scale
heights. This is the first study that combines three-dimensional MHD simulations of zonal flows and

] dust particlesfeeling the gas pressure. The pressurebumps trapparticleswith St=1 veryefficiently.
P
WeshowthatSt=0.1particles(ofsomecentimetersinsizeifat5AUinanMMSN)reachahundred-
E
fold higher density than initially. This opens the path for particles of St = 0.1 and dust-to-gas ratio
. of 0.01 or for particles of St 0.5 and dust-to-gas ratio 10−4 to still reach densities that potentially
h
≥
p trigger the streaming instability and thus gravoturbulentformation of planetesimals.
- Subject headings: magnetohydrodynamics(MHD) - planets and satellites: formation- protoplanetary
o
disks
r
t
s
a 1. INTRODUCTION sub-Keplerian gas (Weidenschilling 1977a). Thus, dust
[ has to grow very quickly from some centimeters to sev-
Planets form as a side product in star formation. The
eral kilometers in size in order to avoid drifting into the
3 general understanding on how planets in our solar sys-
inner region of the protoplanetary disk.
v tem form was detailed in Safronov (1969). Low-mass
Turbulence in protoplanetary disks around young
5 stars form out of molecular clouds which consist of 99%
9 hydrogen and helium (further referred to as gas), and stars provides promising mechanisms for rapid planetes-
0 1% dust and ices (Lodders 2003), i.e., everything that imal formation (Johansen et al. 2007, 2011). Shearing
2 has a higher complexity than hydrogen molecules or he- box simulations (Brandenburg et al. 1995) are a pow-
1. lium. Those molecular clouds have cores between less erful tool for analyzing the magnetorotational instabil-
ity (MRI; Balbus & Hawley 1991, 1998) as a source of
1 than0.1M⊙andmore than10M⊙(Krumholz et al.2012)
turbulence. These simulations consider a local, coro-
2 that are gravitationally unstable. Most parts of the
tating box, representing a small part of a Keplerian
1 mass will collapse into a newborn star. The remaining
disk. Johansen et al. (2009a) reported long-lived ax-
: 1% of the total mass will form an accretion disk with
v ∼ isymmetric sub- and super-Keplerian flows, zonal flows,
pressure supported, sub-Keplerian gas (Weidenschilling
i in shearing box simulations of turbulence caused by
X 1977a;Cassen & Moosman1981)aroundtheyoungstar.
the MRI. These zonal flows have been seen in several
Those disks have lifetimes in the order of a few million
r other local (Fromang & Stone 2009; Stone & Gardiner
a years (Haisch et al. 2001; Fedele et al. 2010). Dust par-
2010; Simon et al. 2012) and global (Lyra et al. 2008;
ticles grow due to coagulation (Weidenschilling 1997).
Dzyurkevich et al. 2010; Uribe et al. 2011; Flock et al.
However, coagulation models show that there are sev-
2011, 2012) simulations using a wide variety of codes.
eral barriers to overcome to grow dust large enough to
Zonal flows are a product of large-scale variations
become gravitationally bound in kilometer-sized plan-
in the magnetic field that transport momentum differ-
etesimals, such as the bouncing barrier (Zsom et al.
entially, creating regions of slightly faster and slightly
2010; Windmark et al. 2012a,b), the fragmentation bar-
slower rotating gas. Large-scale pressure bumps are ex-
rier(e.g.,Beitz et al.2011;Birnstiel et al. 2012, andref-
cited through geostrophic balance. This creates long-
erencestherein),andthekilometer-sizebarrier(Ida et al.
lived over-densities that potentially trap dust particles.
2008; Cuzzi et al. 2008). Dust growth mechanisms are
Amorethoroughdescriptionofzonalflowsandtheircre-
summarized in Dominik et al. (2007) and the review of
ationputforwardinJohansen et al.(2009a)foundzonal
Blum & Wurm (2008) gives an overview on the men-
flows always populating the largest radial mode avail-
tioned barriers.
able in the local box approximation. Their largest box
This paper addresses the fragmentation barrier or
was simulating 10.56 pressure scale heights (H). More
meter-size barrier. Pebbles of several decimeters in size
recentlySimon et al.(2012)foundamorecomplexstruc-
will drift very fast inward due to the headwind from the
ture in their largest simulation with L = 16H. They
x
furtherstudiedtheautocorrelationfunction(Guan et al.
[email protected]
2 Dittrich, Klahr, & Johansen
2009) of the magnetic field and the gas density. Both velocity and by shear flow, respectively. The right-hand
have a two-component structure. The first is tilted with sidecontainstheCoriolisforce,theverticalcomponentof
respect to the azimuthal axis and highly localized. The the stellar gravity,the Lorentz force, the pressure gradi-
second component is seen at the largest scales and can ent,andtheviscosityterm. Here,u(0) = (3/2)Ωxisthe
y
be associated with the (predominantly toroidal) back- Keplerian orbital velocity. The magnetic−field B as well
ground magnetic field. Simon et al. (2012) measure the as the current density J are calculated from the vector
radial length scale of the zonal flows to converge at 6H. potentialAusingB =∇ AandJ =µ−1∇ (∇ A),
In this paper we consider even larger physical extents × 0 × ×
respectively. Here, µ is the vacuum permeability. The
for zonal flow structures. This gives us the opportunity 0
viscosity term f is explained in Section 2.2.1.
to measure physical properties such as size and lifetime ν
We evolve the magnetic potential with the uncurled
independent of the simulated domain. Further, we in-
induction equation
vestigate properties of the zonal flows in radially and
azimuthallystretchedboxes. We alterthe radialandaz- ∂A ∂A 3
+u(0) =u B+ ΩA xˆ+f (A) . (2)
imuthal domain up to 20 gas pressure scale heights. ∂t y ∂y × 2 y η
∼
Additionally, we study the behavior of dust in zonal
The terms on the right-handside express the electromo-
flows. Whipple (1972) was the first to suggest that ax-
tiveforce,thestretching(creationofazimuthalmagnetic
isymmetric pressure bumps can trap gas. Pinilla et al.
field from radial field) by Keplerian shear and the resis-
(2012) invokedzonalflows as a possibility to explainthe
tivity f (see Section 2.2.2).
submillimeter and millimeter-sized particles observed in η
protoplanetary disks. They used artificial static density Thegasdensityisevolvedwiththecontinuityequation
bumps introduced as sinusoidal density perturbations
∂ρ ∂ρ
with different amplitudes (e.g., A = 0.1 and A = 0.3) ∂t +(u·∇)ρ+u(y0)∂y =−ρ∇·u+fD(ρ) , (3)
and different wavelengths (L = 0.3...3H). They found
thata30%densityperturbation(withL=1H)isneces- where the last term on the right-hand side describes
sarytostopthedriftofthedustgrains. Thepresentwork mass diffusion (see Section 2.2.3). We use an isother-
is the first three-dimensional MHD study that combines malequationofstateP =c2ρ, wherethe speedofsound
s
zonal flows and the reaction of dust particles on them. is c =HΩ; H is the gas pressure scale height.
s
Our paper is organized as follows. In Section 2 we
discuss the setup of the simulations in this paper. In 2.2. Dissipation
Section 3 we study the zonal flow properties and their
Maxwell and Reynolds stresses as well as the MRI re-
dependency on the physical box size. The behavior of
lease kinetic and magnetic energy at large scales. This
dust particles in zonal flows is described in Section 4.
energy cascades down to small scales. Since numerical
A discussion and conclusions follows in Section 5 and
simulations have a finite resolution, this small-scale en-
Section 6 provides a summary and an outlook.
ergyneedstobedissipated. Weusenumericaldissipation
inthe formofhyper-andshockviscosity(Section2.2.1),
2. SIMULATION SETUP
hyper-resistivity (Section 2.2.2), and hyper- and shock
We use the Pencil Code,1 a sixth-order spatial and diffusion (Section 2.2.3).
third-ordertemporal finite difference code, for our simu-
2.2.1. Viscosity
lations. WesimulatedthestandardidealMHDequations
in a local shearing box with vertical stratification. The Theviscositytermf inEquation(1)is expressedby
ν
simulation boxes are centered at an arbitrary distance r
to the star. The radial direction is denoted by x, the f =ν ∇6u+ S(3) ∇lnρ
azimuthaldirectionbyy,andthe verticaldirectionbyz. ν 3 ·
The Keplerian frequency is Ω. We include dust particle +νshh[∇∇ u(cid:16)+(∇ u)(∇(cid:17)ilnρ)]
· · ·
dynamics, without back-reactionto the gas and without +(∇ν )∇ u. (4)
sh
·
self-gravity.
We restricted our models to hyper- (ν ) and shock (ν )
3 sh
viscosity. Thus, the regularNavier-Stokesviscosityterm
2.1. Gas Dynamics
isneglected. Thethird-orderrate-of-straintensorS(3) is
The gas velocity u relative to the Keplerian shear is
defined by
evolved via the equation of motion
∂5u
S(3) = i . (5)
∂u ∂u ij ∂x5
+(u ∇)u+u(0) = j
∂t · y ∂y
The high-order Laplacian ∇6 in Equation (4) is ex-
1
2Ωuyxˆ Ωuxyˆ+Ω2zzˆ panded as ∇6 = ∂6/∂x6 +∂6/∂y6 +∂6/∂z6. Further-
− 2
more, the shock viscosity is expressed by
1 1
+ J B ∇P +f (u,ρ) . (1)
ρ × − ρ ν νsh =csh max[−∇·u]+ min(δx,δy,δz)2. (6)
On the left-hand side of the equation, the second and In the fashiono(cid:10)f von Neumann(cid:11)& Richtmyer (1950) it is
third terms are the advection terms by the perturbed proportional to positive2 flow convergence. We take the
1 Details on the Pencil Code and download informationcan be 2 Symbolized by the plus sign in Equation (6). We only apply
foundathttp://www.nordita.org/software/pencil-code/. shockviscositywherethevelocityflowisconverging.
Planetesimal Formation in Zonal Flows 3
maximum over five zones, and smoothed it to the sec- measures are in units of the pressure scale height H =
ond order. As suggested by von Neumann & Richtmyer c Ω−1. Density is stated inunits of the initial mid-plane
s
(1950), we setthe shockviscositycoefficientto c =1.0 gas density ρ . Magnetic field strength is measured in
sh 0
to dissipate energy in shocks at high z above the mid- unitsofc (µ ρ )−1. Energyandstressareinunitsofthe
s 0 0
plane of the disk. mean thermal pressure in the box P =c2 ρ .
h i sh i
Since our simulations are dimensionless, they can be
2.2.2. Resitivity placed at any distance r to the star. Only by defining a
The effects of resistivity are captured by the term globalpressure gradient ∂Pglobal/∂r, which balances the
Coriolis force in
f =η ∇6A, (7)
η 3 1∂P
global
=2Ω∆v, (11)
where η is the hyper-resistivity. ρ ∂r
3
we restrict our simulations to a specific distance to the
2.2.3. Diffusion
starwherethechosenpressuregradientapplies. The pa-
Mass diffusion is computed with rameter∆v =u(0) u isthedifferencetotheazimuthal
y y
−
fD =D3∇6ρ+Dsh∇2ρ+∇Dsh ∇ρ, (8) Keplerian velocity. We fix ∆v = 0.05cs (see also Sec-
· tion 2.6). Numerically, the global pressure gradient acts
where D3 is the hyper-diffusion parameter and Dsh is as an external force on gas and dust.
expanded as in Equation (6).
2.6. Initial Conditions
2.3. Dust Dynamics
The gas density is set to an isothermal hydrostatic
Dust particles are simulated as individual super- equilibrium ρ(z) = ρ exp( z2/2H2). We start with
0
particles i with position xi and velocity vi. Each super- random noise fluctuations −in the gas velocity with
particle position is evolved with δu = 10−3c . The azimuthal component of the
s
magnetic vector potential is initialized with A =
dx(i) y
dt =v(i)+u(y0)yˆ. (9) Ak0=cosk(k=xx4).c7o6sH(k−y1ya)ncdosA(kz=z)0w.0h4ecre(µthρro)u−g1h.out kx =
y z 0 s 0 0
The change of velocity for each particle is evolved Particles are released after the gas turbulence is satu-
through rated. Wemeasuredthistobeafter20T forthelargest
orb
runs. Forconvenience,weusedthesamesaturationtime
dv(i) =2Ωv(i)xˆ 1Ωv(i)yˆ Ω2zzˆ for all our simulations. Particles have a Stokes number
dt y − 2 x − ofSt=τfΩ=1,unlessotherwisestated. Theinitialpar-
1 ticle distribution is Gaussian in z and uniform in x and
v(i) u(x(i)) , (10)
y. The particle velocity is initialized with the station-
−τ −
f
h i ary solution (Nakagawa et al. 1986) for the radial and
where the first and second terms are due to the Coriolis azimuthal velocity
force. Thethirdtermcorrespondstothe verticalgravity
ofthestar. Particlesonlyfeelthegasdrag(thelastterm vx 2∆v
=
in Equation (10)) of nearby cells, but are not subjected c −τΩ+(τΩ)−1
s f f
topressureorLorentzforces. τfdenotesthefrictiontime, vy ∆v
a measure for the size of the particles. = . (12)
c −1+(τΩ)2
s f
2.4. Boundary Conditions We get ∆v from the solution of Equation (11)
For our simulations, we use shearing box boundary
2
conditions in radial (shear-periodic) and azimuthal (pe- ∆v 1 H ∂lnP
= . (13)
riodic) directions. In the vertical direction we also use c −2 r ∂lnr
s (cid:18) (cid:19)
periodicboundaryconditions. Althoughperiodicbound-
We initialized ∆v =0.05c for our simulations.
aryconditionsinverticaldirectionarenotphysical,these s
boundary conditions conserve the average flux of the
2.7. Simulation Parameters
magneticfield. Simulationswithoutflowboundaries(not
includedinthispaper)showednoconsiderablemassflux The parameter space covered by our simulations
across the vertical boundary and did not change the av- is summarized in Figure 1. The vertical extent
erage properties of the zonal flow. is always set to L = 2.64H.3 One simula-
z
tion set (A) covers the boxes with a squared base,
2.5. Dimensions i.e., radial and azimuthal extent are kept the same:
L = L = 1.32, 2.64, 5.28, 10.56, 21.12 H.
We use the dimensionless unit system c = Ω = µ = x y
s 0 { }
These are marked with blue boxes in Figure 1 and
ρ = 1. Velocity is measured in units of the local sound
0 are called runs S, M, L, XL, and XXL. The devia-
speedc . Gasvelocitiesarealwaysdenotedbyuwhereas
s
tion to the global density profile in the largest box
particlevelocitiesarealwaysdenotedbyv. Allvelocities
are differences to the Keplerian orbital velocity v =
K 3 L = 1.32 has been chosen as the basic box size, because
(0) (0)
i(n0,uunyits,0o)f,twhheelroecaulyorb=ita−l(t3im/2e)ΩTx. T=im2πeΩis−1m.eLaesnugrethd sLuxpe=rso1n.3ic2Kapepplreorxiaimnastheelayrmfloawrks(Jtohheantrsaennseittioanl.f2r0o0m9as)u.bsonic to
orb
4 Dittrich, Klahr, & Johansen
Table 1
RunParameters
SimulationSet Run Lx×Ly×Lz Nx×Ny×Nz ν3=η3=D3 nparticles St Shear ∆t
(1) (2) (3) (4) (5) (6) (7) (8) (9)
A S 1.32×1.32×2.64 36×36×72 4.0×10−10 62,500 1.0 FDA 121
A,B,C M 2.64×2.64×2.64 72×72×72 4.0×10−10 250,000 1.0 FDA 121
A,E L 5.28×5.28×2.64 144×144×72 4.0×10−10 1,000,000 1.0 FDA 121
A,E XL 10.56×10.56×2.64 288×288×72 4.0×10−10 4,000,000 1.0 FDA 121
A XXL 21.12×21.12×2.64 576×576×72 4.0×10−10 4,000,000 1.0 FDA 121
B x-S 1.32×2.64×2.64 36×72×72 4.0×10−10 125,000 1.0 FDA 121
B x-L 5.28×2.64×2.64 144×72×72 4.0×10−10 500,000 1.0 FDA 121
B x-XL 10.56×2.64×2.64 288×72×72 4.0×10−10 1,000,000 1.0 FDA 121
C y-S 2.64×1.32×2.64 72×36×72 4.0×10−10 125,000 1.0 FDA 121
C y-L 2.64×5.28×2.64 72×144×72 4.0×10−10 500,000 1.0 FDA 121
C y-XL 2.64×10.56×2.64 72×288×72 4.0×10−10 1,000,000 1.0 FDA 121
D LspecMR 5.28×5.28×2.64 144×144×72 4.0×10−10 1,200,000 0.01...100 FDA 121
D LspecHR 5.28×5.28×2.64 256×256×128 2.0×10−11 120,000,000 0.01...100 FDA 121
D LspecMRs 5.28×5.28×2.64 144×144×72 4.0×10−10 14,000,000 0.01...1.0 FDA 121
D MspecMRb 2.64×2.64×2.64 72×72×72 4.0×10−10 6,000,000 1.0...100 FDA 223
E L SAFI 5.28×5.28×2.64 144×144×72 4.0×10−10 1,000,000 1.0 SAFI 121
E XL SAFI 10.56×10.56×2.64 288×288×72 4.0×10−10 4,000,000 1.0 SAFI 121
Notes. Column1: simulationset. Column2: nameofrun. Column3: boxsizeinunitsofpressurescaleheights. Column4: grid
resolution. Column5: dissipationcoefficients. Column6: numberofparticlesinsimulation. Column7: Stokes numberSt=τfΩ. Column
8: shearadvection scheme. Column9: totalruntimeinorbitsTorb.
yellow in Figure 1). All simulations are stratified and
XXL
have dust particles with different couplings to the gas.
21.12 576×576×72
ThesimulationsdisplayedinFigure1haveparticleswith
N=4×106
p a Stokes number of St=1.
y−XL XL Details on run parameters of those and six more sim-
10.56 72×288×72 288×288×72 ulations are found in Table 1. The first set of simula-
N=106 N=4×106 tions (A,B, and C) in Table 1 are the simulations with
p p
medium resolution, i.e., 36 grid cells4 per 1.32 pressure
y−L L
H scale heights. Simulation set D was carried out to in-
/y 5.28 72×144×72 144×144×72 vestigate the behavior of different particle sizes in the
L N=5×105 N=106
p p presence of zonal flows. Run LspecMR is very much like
x−S M x−L x−XL runL,butwith12differentparticleStokesnumbers. The
2.64 36×72×72 72×72×72 144×72×72 288×72×72 run LspecHR has a resolutionof 64 grid cells per 1.32H.
N=125,000 N=250,000 N=5×105 N=106 Runs LspecMR and LspecHR have 12 different parti-
p p p p
cle species, with Stokes numbers of St = 0.01...100.0.
S y−S
The runs LspecMRs and MspecMRb have particles with
1.32 36×36×72 72×36×72 StokesnumbersofSt=0.01...1.0andSt=1.0...100.0,
N=62,500 N=125,000
p p respectively. These two simulations were carried out to
studyparticlebehaviorwithmoreparticlespergridcell5
1.32 2.64 5.28 10.56 21.12 atmediumresolution. Thecorrespondingsizesfordiffer-
L /H ent protoplanetarydisk models are found in Section 5.2.
x
Simulation set E is a comparison of runs L and
Figure 1. Parameterspaceofradialandazimuthalboxsizesthat
wassimulatedforthispaper. Everysimulationhasaverticalextent XL to the same runs (L SAFI and XL SAFI) with
of 2.64H. The first line in each box states the name of the run, the Shear Advection by Fourier Interpolation (SAFI)
the second line the number of grid cells used, and the third line scheme. Here,allvariablesq(x,y,z)aretransformedinto
givesthenumberofsimulatedsuper-particles. Moredetailsonall
Fourier space in the y-direction to get qˆ(x,k ,z). Then
simulationsarefoundinTable1. y
each Fourier mode is multiplied by exp[ik u(0)(x)δt] to
y y
(0)
shift by u (x)δt in real space and is inverse Fourier
can be quite severe at the inner and outer bound- y
transformed to real space. This method reduces the ad-
ary of the largest simulation. Thus, the results from
vectionerrortothestandardFiniteDifferenceAdvection
run XXL have to be treated with caution. Another
(FDA) scheme in the PencilCode (more details onFDA
set of simulations (B) varies the radial size of the box,
and SAFI are found in Johansen et al. 2009a).
L = 1.32, 2.64, 5.28, 10.56 H, with constant box
x
{ }
size in azimuthal direction, L = 2.64H. This set is
y
marked red in Figure 1 and includes runs x-S, M, x-L, 4 Wechose36gridcellsinsteadoftheusual32gridcells. That
choicewasdoneduetothearchitecture(12CPUspernode)ofthe
and x-XL. The third set of simulations (C) varies the
usedcluster,THEOintheMPGcomputingcenterinGarching.
azimuthal extent, Ly = 1.32, 2.64, 5.28, 10.56 H, 5 ∼ 9 and ∼ 16 particles per grid cell for runs LspecMRs
{ }
while the radial extent is kept constant, Lx = 2.64H. and MspecMRb compared to ∼ 0.8 particles per grid cell for
This set includes runs y-S, M, y-L, and y-XL (marked runLspecMR.
Planetesimal Formation in Zonal Flows 5
Table 2
TurbulenceProperties
Run h21u2xi h12u2yi h21u2zi h12Bx2i h12By2i h12Bz2i hρuxuyi h−BxByi α
(1) (2) (3) (4) (5) (6) (7) (8) (9) (10)
S 2.1×10−3 3.3×10−3 1.5×10−3 8.8×10−4 6.4×10−3 3.5×10−4 7.7×10−4 3.4×10−3 2.8×10−3
M 3.9×10−3 5.2×10−3 2.2×10−3 1.9×10−3 1.2×10−2 7.6×10−4 1.6×10−3 6.6×10−3 5.5×10−3
L 5.0×10−3 5.6×10−3 2.3×10−3 2.0×10−3 1.3×10−2 8.0×10−4 1.9×10−3 6.9×10−3 5.9×10−3
XL 5.2×10−3 5.0×10−3 2.1×10−3 1.7×10−3 1.1×10−2 6.7×10−4 1.8×10−3 6.1×10−3 5.2×10−3
XXL 5.1×10−3 4.6×10−3 2.0×10−3 1.6×10−3 1.0×10−2 6.2×10−4 1.7×10−3 5.7×10−3 4.9×10−3
x-S 4.0×10−3 5.2×10−3 2.4×10−3 2.0×10−3 1.3×10−2 8.4×10−4 1.7×10−3 7.0×10−3 5.8×10−3
x-L 3.8×10−3 5.2×10−3 2.1×10−3 1.7×10−3 1.1×10−2 6.9×10−4 1.5×10−3 6.1×10−3 5.1×10−3
x-XL 3.8×10−3 5.0×10−3 2.2×10−3 1.8×10−3 1.2×10−2 7.1×10−4 1.5×10−3 6.2×10−3 5.2×10−3
y-S 2.3×10−3 3.7×10−3 1.6×10−3 1.0×10−3 7.1×10−3 4.1×10−4 8.6×10−4 3.8×10−3 3.1×10−3
y-L 5.2×10−3 5.7×10−3 2.5×10−3 2.2×10−3 1.4×10−2 8.8×10−4 2.0×10−3 7.5×10−3 6.3×10−3
y-XL 5.4×10−3 5.0×10−3 2.2×10−3 1.8×10−3 1.2×10−2 7.2×10−4 1.9×10−3 6.4×10−3 5.6×10−3
LspecMR 4.8×10−3 5.3×10−3 2.2×10−3 1.9×10−3 1.2×10−2 7.5×10−4 1.8×10−3 6.6×10−3 5.6×10−3
LspecHR 3.0×10−3 4.2×10−3 1.4×10−3 1.3×10−3 7.5×10−3 5.4×10−4 1.0×10−3 4.2×10−3 3.5×10−3
LspecMRs 5.3×10−3 6.2×10−3 2.5×10−3 2.2×10−3 1.4×10−2 9.2×10−4 2.1×10−3 7.7×10−3 6.5×10−3
MspecMRb 4.3×10−3 5.7×10−3 2.5×10−3 2.1×10−3 1.4×10−2 8.7×10−4 1.8×10−3 7.2×10−3 6.0×10−3
L SAFI 5.1×10−3 5.8×10−3 2.4×10−3 2.1×10−3 1.3×10−2 8.4×10−4 2.0×10−3 7.2×10−3 6.1×10−3
XL SAFI 5.4×10−3 5.3×10−3 2.2×10−3 1.8×10−3 1.2×10−2 7.3×10−4 1.9×10−3 6.5×10−3 5.6×10−3
Notes. Column1: nameofrun. Columns2-4: kineticenergy. Columns5-7: magneticenergy. Column8: Reynoldsstress. Column9:
Maxwellstress. Column10: α-value,followingEquation(14). Stressesandenergieshavebeennormalizedtothemeanthermalpressure
inthebox,hPi=c2hρi.
s
Every simulation is run for 121 local orbits T = Simon et al. (2012):
orb
2πΩ−1, except for LspecMRb which runs for 223T in
orb
order to follow the evolution of the slowly settling large 2π va,j
Q = | |, (15)
particles. After 20Torb, when the initial conditions are j Ω∆xj
sufficiently forgotten and the turbulence saturated, the
particles are started. where the Alfv´en speed is defined as va,j 2 = Bj 2/ ρ .
| | h i h i
The notation x denotes volume averaging, x shows a
h i
time average. Sorathia et al. (2012) show that Q &
z
3. ZONAL FLOW PROPERTIES 10 15 for poorly resolved azimuthal quality factors
−
(Q 10) is required to resolve the MRI. Larger val-
TurbulencepropertiesaresummarizedinTable2. The y ∼
ues of the azimuthal quality factor (Q & 25) allow for
kinetic andmagneticenergyaswellas the Reynoldsand y
lowerverticalquality factors. The azimuthalcomponent
Maxwell stress almost doubles when increasing the box
of the magnetic field is very well resolved (Q & 25) for
sizefrom(1.32H)2 2.64H (runS)to(2.64H)3(runM). y
× all simulations, but runs S and y-S. The vertical com-
Further increasing the box size does not change the re-
ponent has values between 6 and 8. We thus conclude
sultingenergiesandstressesbymuch. Theradiallyshort
that all simulations, but runs S and y-S have sufficient
boxofrunx-S with1.32H (2.64H)2 hassimilarresults
× resolution for the MRI.
on these values. However, the azimuthally short box of
In Figure 2 a snapshot of the runs y-XL, XL, M, and
run y-S has turbulent energies and stresses comparable
x-XL are shown in scale, giving a real size comparison
to run S. These measurements show that the turbulence
of high- and low-pressure regions. The large-scale sinu-
parameters are saturated for boxes with an azimuthal
soidal form of the dominant mode is observable in these
extent of at least 2.64H. This confirms the results from
plots. Thehighermodesaremuchshorterlivedandseem
Fromang & Stone (2009) who found that the turbulence
to be non-axisymmetric density waves affected by the
properties do not change when the box size is increased
shear (Heinemann & Papaloizou 2009). The amplitudes
radially,iftheazimuthaldimensionislargeenough. The
ofthepressuredifferencesarehigherinazimuthallylarge
α-value (Shakura & Sunyaev 1973) in Column 10 in Ta-
boxes. Only the axisymmetric density waves are long-
ble 2 is calculated via
lived and strong enough to make up a significant contri-
bution to the pressure bump structure in an azimuthal
2( ρu u B B )
α= h x yi−h x yi , (14) aswellasatemporalaverageoversomelocalorbits. The
3 P
azimuthal average of Figure 2 is seen in Figure 3. Here,
h i
theaxisymmetricstructureisclearlyvisibleandastrong
where P = c2 ρ . The factor of 2/3 originates from correlation between the particle location and a positive
h i sh i
theshearparameterq = dlnΩ/dlnR. Weuseq =3/2, radialgradientofthe gasdensity isseen. Theblackdots
−
appropriate for a Keplerian disk. For further details see in Figures 2 and 3 show the radial and azimuthal posi-
Brandenburg et al.(1995,page748). TheMaxwellstress tion of every 100th particle. The particles are trapped
isaroundthreetimeshigherthantheReynoldsstressand by the axisymmetric pressure bumps. Also, the shapes
thus dominates the α-value. of spiral density waves (Heinemann & Papaloizou 2009)
In order to verify that our numerical resolution is suf- can be seen in the structures. The particle distribution
ficient, we examined the quality factor as described in with respect to the gas flow will be discussed in more
6 Dittrich, Klahr, & Johansen
y−XL XL t= 85.0T
orb
5.3 1.30
3.9 1.24
2.6 1.18
1.3
1.12
H 0.0
y/ 1.06
−1.3 =0)
1.00Σ(t
Σ /
−2.6
0.94
−3.9
0.88
−5.3
M x−XL
0.82
1.3
H 0.0 0.76
y/
−1.3 0.70
−1.3 0.0 1.3 −5.3 −3.9 −2.6 −1.3 0.0 1.3 2.6 3.9 5.3
x/H x/H
Figure 2. Collageoffourgassurfacedensityrepresentations oftherunsy-XL,XL,M,andx-XL.Eachsnapshotwastakenafter85Torb.
Theseplotsshowthatgasoverdensitiesaremostpronouncedinthelargestbox. Thenon-axisymmetricstructureshaveveryshortlifetimes;
less than a tenth of an orbit. The pressure bump structures are more visible when the density is averaged over the azimuthal direction
(Figure3). Theblack dots represent the position of every100th particle, integrated invertical direction. Theparticles aretrapped both
inaxisymmetricpressurebumpsandinspiraldensitywaves asdescribedinHeinemann&Papaloizou(2009).
detail in Section 4. metric pressure bumps get weaker too.
All of our simulations show signs of zonal flows. Figure5additionallyshowsspace–timeplotsoftheaz-
Strength and lifetime of the zonal flows and the asso- imuthal gas velocity and the radialgradients of gas den-
ciated pressure bumps differ very much with the physi- sity and azimuthal gas velocity of run XXL. In the right
cal box size. Space–time plots of all different simulation panels,thepositionofthehighestdustdensityareshown
sizes are shown in Figure 4 and the upper left panel of as dots. Particles get clearly slowed by the maxima of
Figure 5. The pressure bumps are generally more pro- the azimuthal gas velocity, i.e., the large-scale maxima
nounced in simulations with a larger radialextent. Sim- of the pressure gradients. The velocity has large-scale
ulationsetAstrictlyfollowsthisgeneraltrend. Pressure structures that are very similar to those of the density
bump features grow in strength and lifetime with the gradient, as expected in geostrophic balance. Thus, the
physical box size, always staying at the largest radial structureofthevelocitygradientcanbeapproximatedas
scale. This rule applies to all but the largest runs XL, the large-scale structure of the second derivative of the
x-XL, and XXL. There, instead of the formerly predom- gasdensity. It is shownin the lowerrightpanelthat the
inant k = 1 (ω = 2πk /L in Fourier mode sin(ω x)) particles get stopped at the minima of the radialderiva-
x x x x x
mode, the mode k = 2 (higher modes for run XXL) is tiveoftheazimuthalgasvelocity(andthusatminimaof
x
occupied by the pressure bumps. In simulation set B, the second derivative of the gas density) as analytically
the strength and size of the pressure bumps converges predicted (see, e.g., Klahr & Lin 2001).
for simulations with a radial extent larger than 5.28H. We calculated the correlation time of the pressure
Thelifetimeofthepressurebumpsevendecreasesforthe bumps and the zonal flows in the same way as it was
largestsimulationinthisset. SimulationsetC isqualita- calculated in Johansen et al. (2009a). We use the den-
tively different from the other simulationsets. Strength, sityρ,averagedoverazimuthalandverticaldirections,at
size, and lifetime of the pressure bumps seem to be in- agiventimet. Then,weaverageovereachpointinradial
versely proportional to the azimuthal extent of the box, direction the time it takes for the density at each point
when the verticalandradialbox sizes arekept constant. tochangebyavaluecorrespondingtothestandarddevi-
This effect was already seen in, e.g., Simon et al. (2012) ation of the gas density. These measurements are taken
andFlock et al.(2012). Bothgroupsshowthatthemag- for every local orbit. The measurements are averaged
netic field consists of two components: a local turbulent over the time between saturation of the turbulence and
component that is responsible for the zonal flows and a time when the correlation does not extend the corre-
a global azimuthal component. Since the total energy lation time to the final time of the simulation. Finally,
stays approximately constant, the local component gets the averagesare multiplied by two, in order to cover the
weaker and consequently zonal flows as well as axisym- full temporal extent of the correlated structures. The
Planetesimal Formation in Zonal Flows 7
y−XL XL t= 85.0T
orb
5.3 1.15
3.9 1.12
2.6 1.09
1.3
1.06
H 0.0
y/ 1.03
0)
=
−1.3 1.00Σ> / (ty
−2.6 Σ<
0.97
−3.9
0.94
−5.3
M x−XL
0.91
1.3
H 0.0 0.88
y/
−1.3 0.85
−1.3 0.0 1.3 −5.3 −3.9 −2.6 −1.3 0.0 1.3 2.6 3.9 5.3
x/H x/H
Figure 3. Surface density distribution of Figure 2, averaged in azimuthal direction and averaged over the mean surface density. This
revealsaxisymmetricpressurebumpsandvalleys. Particlesaretrappedontheinnersideofthedensitymaxima,atplaceswithapositive
density gradient to overcome the negative global pressuregradient. These pressurebumps are stable formany orbits (compare Figures 4
and5).
correlation times measured in this fashion are in good The runs L SAFI and XL SAFI were carried out to
agreement with the lifetime of the overdensities that is compare the turbulence and zonal flow parameters with
seeninFigures4and5. However,achangeofpositionof the runs L and XL. They were run to check that zonal
thestructures,asseeninrunXL(checkFigure4),isnot flowsarenoeffectsfromtheshearadvectionschemethat
accountedfor. Thus,correlationtimesaremorelikelyto was used in the Pencil Code. Comparing the values in
be underestimated than overestimated. Also, we cannot Tables 2 and 3 shows that there is little change in the
be entirely sure whether this behavior is really drift or measured properties of the zonal flows and the associ-
structure decay and reformation. ated pressure bumps. However, the computation time
The results of the correlation time determination are increasesifoneusestheSAFIscheme. Thus,thisscheme
showninTable3andintheupperpanelofFigure6. For was only used to confirm our results.
the diagonalsimulationset, (A), the correlationtime in-
creases with box size. It seems to saturate toward the
4. PARTICLE BEHAVIOR IN ZONAL FLOWS
largestbox size. The trendto longercorrelationtimes is
alsoevidentforsimulationsetB.Hereonlyrunx-XLhas Particleaccumulationsandplanetesimalformationcan
a shorter correlationtime than expected. This might be occurin clumps andfilaments ofthe overdensitiesin the
an effect of the strongly stretched simulation box. The dust. Inoursimulations,wedonotincludegravitational
correlationtime decreasesslightly withanincreasingaz- interaction between the particles. Thus, we only study
imuthal box size in simulation set C (not shown in the the passively developed overdensities of the dust to see
figure). The lower panel in Figure 6 shows a measure- whenandwhetheroverdensitiessufficientforthestream-
ment for the physical size of the zonal flow features. We ing instability can be reached. By not having explicit
Fourier-transformed the vertically and azimuthally av- feedback one can retroactively study the concentration
eraged gas density and azimuthal gas velocity for each factor for various initial dust-to-gas ratios. Simulations
time step and averaged the amplitudes of the first four including feedback will have to be done in future stud-
modesoverthetimeof20...120localorbits. Thelength ies. Figure 2 shows the position of every 100th parti-
was normalized for the size of the simulation box, to get cle in selected simulations. These plots clearly show the
the physical size of the modes by λ = L /k with the trend for particles to accumulate in the downstream of
x x x
wave number k . The turbulence is always strongest at high-pressure regions. Particles are pulled toward pres-
x
the largest modes for simulations with L . 5H. The sure gradientmaxima (Klahr & Lin 2001). In the upper
x
highest amplitude for both quantities in the largestsim- right panel in Figure 2, a snapshot of run XL after 85
ulationdomainisfoundbetween5H and7H (upto10H local orbits is shown. The particles clump up at posi-
for ρˆ). These measurements are also found in Table 3. tions just left of the maxima in the kx = 1 of the gas
density; these are the locations of positive zonal flows,
8 Dittrich, Klahr, & Johansen
Table 3
ZonalFlowProperties
Run ρrms |ρˆ(kx=1)| |ρˆ(kx=2)| |ρˆ(kx=3)| |uˆ˜y(kx=1)| |uˆ˜y(kx=2)| |uˆ˜y(kx=3)| τcorr
(1) (2) (3) (4) (5) (6) (7) (8) (9)
S 6.1×10−3 4.1×10−3 7.5×10−4 4.0×10−4 9.9×10−3 3.4×10−3 2.1×10−3 7.6
M 2.0×10−2 1.0×10−2 2.5×10−3 1.6×10−3 1.2×10−2 5.2×10−3 3.1×10−3 11.2
L 3.9×10−2 2.1×10−2 5.6×10−3 3.2×10−3 1.3×10−2 6.5×10−3 3.7×10−3 23.2
XL 4.3×10−2 1.5×10−2 1.1×10−2 5.0×10−3 4.6×10−3 6.6×10−3 4.3×10−3 43.2
XXL 4.0×10−2 5.0×10−3 7.9×10−3 7.9×10−3 7.8×10−4 2.4×10−3 3.5×10−3 47.3
x-S 1.6×10−2 3.9×10−3 1.7×10−3 1.2×10−3 8.0×10−3 3.2×10−3 2.0×10−3 4.4
x-L 3.3×10−2 2.3×10−2 7.3×10−3 2.7×10−3 1.4×10−2 8.8×10−3 4.7×10−3 37.6
x-XL 2.7×10−2 1.2×10−2 1.0×10−2 6.4×10−3 3.5×10−3 6.2×10−3 5.6×10−3 20.2
y-S 1.2×10−2 9.9×10−3 2.4×10−3 1.1×10−3 1.2×10−2 5.9×10−3 4.1×10−3 14.4
y-L 3.0×10−2 6.9×10−3 3.8×10−3 3.1×10−3 7.8×10−3 3.7×10−3 2.3×10−3 10.8
y-XL 3.6×10−2 5.6×10−3 4.1×10−3 3.3×10−3 5.7×10−3 2.6×10−3 1.6×10−3 10.3
LspecMR 3.7×10−2 1.8×10−2 5.5×10−3 3.1×10−3 1.1×10−2 6.1×10−3 3.6×10−3 21.8
LspecHR 4.2×10−2 9.8×10−3 3.2×10−3 2.0×10−3 5.8×10−3 3.6×10−3 2.0×10−3 23.4
LspecMRs 4.3×10−2 2.4×10−2 5.4×10−3 3.2×10−3 1.4×10−2 6.1×10−3 3.9×10−3 10.9
MspecMRb 1.9×10−2 8.6×10−3 2.6×10−3 1.6×10−3 1.0×10−2 5.3×10−3 3.1×10−3 26.4
LSAFI 3.9×10−2 2.1×10−2 5.3×10−3 3.2×10−3 1.2×10−2 6.1×10−3 3.9×10−3 25.6
XLSAFI 4.8×10−2 2.4×10−2 1.0×10−2 5.4×10−3 7.0×10−3 6.2×10−3 4.8×10−3 48.6
Notes. Column1: nameofrun. Column2: root-mean-squaredensityρrms=ph(ρ−ρ)2i. Columns3-5: Fourieramplitudeofradial
densitymodeskx=1...3,normalizedbymeandensityinthebox. Columns6-8: Fourieramplitudeofazimuthalvelocitymodes
kx=1...3withu˜y =uy−uy. Column9: correlationtime,inorbitsT =2πΩ−1,ofthelargestradialdensitymode.
i.e., regions where the azimuthal gas velocity is higher one can see a clear correlation between the zonal flow
than the pressure-supportedKeplerian flow. strengthandthegasdensityenhancements. Thebottom
panelinFigure7showstheevolutionoftheα-parameter,
4.1. Particles in Zonal Flows calculated as in Equation (14).
Themaximumofthedustoverdensitythatoccursdur-
IntheupperrightpanelofFigure5,theazimuthalgas
ing one simulation is plotted against the box size in the
velocity development of run XXL is shown, overplotted
upper panel of Figure 8. The general trend shows that
with the position of the most massive clump for each
radiallylargerboxeshavehigherparticleconcentrations.
time step. The azimuthal gas velocity coincides with
Anincreasedazimuthalextentdoesnothaveaneffecton
the derivative of the gas density, but it is much eas-
the particle concentrations. The most surprising result
ier to interpret. The speckled structure of the deriva-
is in runy-S. It showsa very highparticle concentration
tives comes due to the high power in the smaller scales.
that occurs early in the simulation (compare Figure 7).
However,the large-scale structure is still visible and the
This is most likely a stochastic coincidence. The lower
geostrophiccorrelationbetweenthestructuresofu (x,t)
y panelinFigure8showsaplotofthemaximumdustover-
and d/dx[ρ(x,t)] is directly observed. Since they have
density against the correlation time of the zonal flows.
the same large-scale structure, the particle position is
The error margin are calculated with the standard de-
mucheasierinterpretedattheazimuthalgasvelocityplot
viation of the temporal evolution of the two quantities.
than on the density gradient plot. Sometimes the radial
We see a clear trend that denser particle accumulations
displacement from one orbit to the next is too large to
develop with longer correlation times. The distribution
be explained by radial drift. That happens when an-
can be fitted by a power law. This gives an exponent of
other clump becomes more massive than the previous
dlogρ/dlogτ =0.38 0.05. The one point that does
one. These particlesaccumulate in regionswith highaz- corr ±
notoverlapwiththeerrormarginsofthefitisfromruny-
imuthalgasvelocities(seeupperrightpanelinFigure5).
S. If we take the maximum of the top panel in Figure 7
The only time when this is not true is at times from 80
after the two first maxima (i.e., after 45T ) and plot
to 100 local orbits. In this period, an inward-drifting orb
this value again in the parameter space of Figure 8, we
clump stayed coherent during the time of its drift. The
get the position marked with the blue square. It agrees
drift velocity of the most massive particle clump is indi-
well with the error margins of the fit.
rectly encryptedin this plot. Particlesare drifting much
In isothermal geostrophic balance, 2ρΩu = c2∂ρ/∂x,
slowerwhenthey aretrappedbyapressuregradient. As y s
the azimuthalgasvelocityfollowsthe radialdensity gra-
all particles drift inward this leads to accumulation of
dient. That this is true for large scales as shown in Fig-
particles in regions where the perturbed pressure gradi-
ure 9. The upper left panel shows the evolution of the
ent is positive.
azimuthally and vertically averaged azimuthal compo-
The maximal accumulation of particles for runs XL
nentofthe gasvelocity. Overplottedarethe locationsof
and y-S are plotted in the top panel in Figure 7.
themaximainthedustdensity. Intheupperrightpanel
The second panel shows the evolution of the quantity
the dust density evolution of the same run L is plotted.
u u 2,ameasureforthestrengthofthezonal
h y−h yiyzix In comparing the location and times of the maxima and
flows. The third panel in Figure 7 shows the evolu-
minima on these two plots, one clearly sees that max-
ption of the strength of the gas density enhancement as
ima in the dust density occur often at times and loca-
ρ ρ 2. Comparing the second and third panels,
h −h iyzix
p
Planetesimal Formation in Zonal Flows 9
120 120 120 120
S x−S y−S M
100 100 100 100
Torb 80 Torb 80 Torb 80 Torb 80
t/ 60 t/ 60 t/ 60 t/ 60
40 40 40 40
20 20 20 20
−0.64 0.00 0.64 −0.64 0.00 0.64 −1.30 −0.65 0.00 0.65 1.30 −1.30 −0.65 0.00 0.65 1.30
x/H x/H x/H x/H
120 0.15 120
y−L L
100 0.12 100
t/Torb 6800 00..0069ρ(t)> t/Torb 6800
40 < 40
/
20 0.03 ] 20
>
−1.30 −0.65 0.00 0.65 1.30 0.00 (t) −2.62 −1.97 −1.31 −0.66 0.00 0.66 1.31 1.97 2.62
x/H ρ x/H
<
120 y−XL −0.03− 120 x−L
100 −0.06) 100
t
Torb 80 −0.09ρ(x, Torb 80
t/ 60 [ t/ 60
−0.12
40 40
20 −0.15 20
−1.30 −0.65 0.00 0.65 1.30 −2.62 −1.97 −1.31 −0.66 0.00 0.66 1.31 1.97 2.62
x/H x/H
120
XL
100
Torb 80
t/ 60
40
20
−5.26 −4.60 −3.95 −3.29 −2.63 −1.97 −1.32 −0.66 0.00 0.66 1.32 1.97 2.63 3.29 3.95 4.60 5.26
x/H
120
x−XL
100
Torb 80
t/ 60
40
20
−5.26 −4.60 −3.95 −3.29 −2.63 −1.97 −1.32 −0.66 0.00 0.66 1.32 1.97 2.63 3.29 3.95 4.60 5.26
x/H
Figure 4. Evolution of the gas density perturbation of all runs from simulationsets A, B, and C. Run XXL is shown inFigure 5. The
density is averaged in vertical and azimuthal direction and plotted in radial direction over time. The lifetime, the size, as well as the
strength of the pressure bumps are clearly increasing with increasing box size in simulation set A, i.e., runs S, M, L, XL, and XXL. In
simulationset B, i.e., runs x-S,M, x-L,and x-XL,we have the same increase of lifetime, size, and strength of the pressurebumps. Only
fortheverylargesimulation(Lx=10.56H),thereisnoapparentdifferenceinpressurebumpsizeandstrengthtorunx-L.Forsimulation
setC,i.e.,runsy-S,M,y-L,andy-XL,thestrengthofthepressurebumpsisapparentlyconstantthroughoutthissetofsimulations. Even
thelifetimedecreases slightlywithincreasingboxsize.
tions where one finds maxima in the gas velocity. Two points. This is indicated by the color scale, showing the
attempts to quantify this observation are shown in the amountof points in eachof the boxes in the scatter plot
lower row of Figure 9. In the left panel, the particle space. There is a clear trend for high dust density con-
density and the azimuthal velocity from the two upper centrations to appear at high gas velocities. Without
panels are plotted against each other, regardless of po- radial drift particles would concentrate where u = 0,
y
sition and time. In the right panel, a snapshot of the i.e., between the sub- and super-Keplerian flow. Due
simulation (as in Figure 2) was taken at 85 local orbits, to the radial drift particles accumulate slightly down-
the time when the maximum dust density enhancement stream at the formed pressure bumps. Those happen to
occurs. Theparticledensity aswellasthe azimuthalgas beatthemaximaoftheazimuthalgasvelocity. Withthe
velocitywere integratedin verticaldirectionand plotted geostrophic balance, high velocities are also regions of a
againsteachother,regardlessoftheirradialorazimuthal high radial density gradient. These plots prove that the
position in the simulation, in this scatter plot. In order particlesinthesimulationsaretrappedbythelong-lived
to visualize high densities of points in these plots, we pressure gradients that occur due to stable zonal flows.
computed a two-dimensional histogram of the scattered If the dust-to-gas ratio increases to values larger than
10 Dittrich, Klahr, & Johansen
120 0.15 120 0.060
0.12 0.048
100 0.09 100 0.036
>
0.06 (t) 0.024
ρ
< cs
80 0.03 >]/ 80 0.012 v]/
Torb 0.00ρ(t) Torb 0.000∆)+
t/ 60 −0.03− < t/ 60 −0.012(x,t
) uy
−0.06x,t −0.024[
(
ρ
[
40 −0.09 40 −0.036
−0.12 −0.048
20 −0.15 20 −0.060
−10 −5 0 5 10 −10 −5 0 5 10
x/H x/H
120 0.10 120 0.10
0.08 0.08
100 0.06 100 0.06
0.04 ) 0.04
>
t/Torb 6800 −− 000.0..00220ρρ((x,t)/<(t) t/Torb 6800 −− 000.0..00220dx(u(x,t)/c)ys
dx d/
−0.04d/ −0.04
40 −0.06 40 −0.06
−0.08 −0.08
20 −0.10 20 −0.10
−10 −5 0 5 10 −10 −5 0 5 10
x/H x/H
Figure 5. ToppanelsshowtheevolutionofthegasdensityperturbationandtheazimuthalgasvelocityofrunXXL,whereasthebottom
rowshowstheradialderivativeofthesequantities. Thederivativesareveryspeckled,sincesmall-scalefluctuationsgivestrongeramplitudes
to the derivatives. However, the underlying large-scale structure is still visible. The azimuthal gas velocity follows the radial gas density
gradient,asexpected forageostrophicbalance. Hence, itispossibletointerprettheradialderivativeoftheazimuthal gasvelocityasthe
secondderivativeofthe gasdensity. Intheupperrightpanel,theblackdots representthepositionofthemostmassiveparticleclumpin
thesimulationateachtime. Itisclearlyshownthatparticlesgettrappedinregionsofpositivezonalflowdownstreamofpressurebumps.
unity, the streaming instability (Youdin & Goodman they would in a laminar disk. However, the size of the
2005; Johansen & Youdin 2007; Youdin & Johansen simulation has little effect on the actual drift velocity,
2007)istriggered. Thisincreasesthedustdensityfurther as shown in the lower panel of Figure 10. It shows a
on timescales shorter than an orbital period. To follow timeaverageofthe particledriftvelocityplottedagainst
the streaminginstabilitydevelopment,the back-reaction the box size. The uncertainties are too large to reveal
ofthedustparticlestothegasphasemustbeconsidered a trend. Thus, the reduction of the radial drift velocity
infuturenumericalsimulation. Thiseffectwasneglected apparently only depends on the amplitude of the zonal
in this set of simulations. Otherwise the initial dust-to- flow, but not on the correlation time. Looking at the
gas ratio would have been an additional free parameter largest run XXL, we can estimate that the radial drift
to be studied. gets reduced by about 28% (drop of the absolute value
from 0.05c to (0.036 0.003)c ).
s s
±
4.2. Radial Drift
Radial drift velocities of the particles in the simula-
4.3. Clustering
tions with different box sizes are shown in Figure 10.
The upper panel shows the measured and expected ra- The clustering degree of the particle distribution can
dial drift of two simulations (M and XL). They show be estimated with the distribution of the dust surface
that particles drift slower in turbulent simulations than density Σ (Pan et al. 2011). The initial distribution is
p