Diffusion coefficient of a passive contaminant in a local MHD model of a turbulent accretion disc

0.24 MB·English

by  Augusto Carballido

Mon.Not.R.Astron.Soc.000,000–000 (0000) Printed2February2008 (MNLATEXstylefilev2.2)
Diffusion coefficient of a passive contaminant in a local
MHD model of a turbulent accretion disc
Augusto Carballido1,2,⋆ James M. Stone2,† and James E. Pringle1‡
1 Institute of Astronomy, Universityof Cambridge, Madingley Road, Cambridge CB3 0HA, UK
5 2 Department of Astrophysical Sciences, Princeton University,Princeton NJ 08544, USA
0
0
2

n 2February2008
a
J

1 ABSTRACT
2 We calculate the radial diffusion coefficient for a passive contaminant in an accretion
disc which is turbulent due to the action of the magnetorotational instability. Nu-
1
merical MHD simulations are used to follow the evolution of a local patch of the disc
v
usingtheshearingboxformalism.Aseparatecontinuityequationforthemassfraction
4
of contaminant is integrated along with the MHD sytem, and radial profiles of this
7
fraction are obtained as a function of time. Solutions of a linear diffusion equation
4
are fitted to the numerical measured profiles of the contaminant, treating the diffu-
1
0 sion coefficient D as the fitting parameter. At early times, the value of D is found to
5 vary, however once the contaminant is spread over scales comparable to the box size,
0 it saturates at a steady value. The ratio of D to the transport coefficient of angular
/ momentum due to shear stress is small. If D canbe used as a proxy for the turbulent
h
p magnetic diffusivity, the effective magnetic Prandtl number Peff = ν/D (where ν is
- the coefficient of ”effective viscosity” associated with shear stress) would be large.
o
r Key words: accretion: accretion discs – magnetohydrodynamics
t
s
a
:
v
Xi 1 INTRODUCTION magnetohydrodynamicalturbulencedrivenbythemagneto-
rotational instability (MRI) (Balbus & Hawley 1991).
r In order to understand the morphology and chemical com-
a The transport of a passive scalar has been used as a
positionofprotoplanetaryaccretiondiscs,theradialmixing
criteriontostudydiffusiveprocessesinturbulentflows(e.g.
of vapour and solid material in accretion flows needs to be
Brandenburg, K¨apyl¨a & Mohammed 2004 and references
characterised(Morfill1983).Thedistributionofdustacross
therein). Research on the behaviour of passive tracers in
aturbulentdisc,asaconsequenceofbothadvectiveanddif-
turbulencehasshownthatthereisastrongvariabilityofdis-
fusive modes, will determine the mechanism which is most
sipation and mixing rates of both scalar and velocity fields,
favourablefortheassemblageofparticlesintoplanetesimals,
due to small scale intermittency. It has also been shown
meteoriticbodiesandeventuallyplanets.Forexample,there
that the central assumption of Kolmogorov scaling, that at
is evidence for radial upstream movement of nebular gas
smallscalesthescalarandvelocityfieldsareisotropicinthe
from the coexistence of calcium-aluminum inclusions and
limit of infiniteReynoldsand Pecl´et numbers,breaks down
chondrules in meteorites (Boss 1998). Studies by Takeuchi
(Warhaft 2000; Sreenivasan 1991). Although a fascinating
&Lin (2002) and byKeller&Gail (2004) indicatethat,for
verticallyisothermallaminardiscs,theremayexistanequa- subject, in this work we will not address many important
issues that belong to thevast field of passive scalars.
torialoutflowcapableoftransportingmaterialfromthehot,
innerregionstothecold,outerregions,therebyaffectingthe Clarke & Pringle (1988) analysed the evolution of the
composition of bodieslocated in thelatterzone.Themate- concentration of a trace contaminant in a steady, Keple-
rial is assumed to be well coupled to the gas, which is true rian disc. They point out a subtlety related to the identi-
for particles whose size is smaller than the mean free path fication of a diffusion coefficient with a viscosity for tur-
ofthegas molecules. However,thispicturemaybechanged bulent accretion. Despite the fact that both are associ-
if angular momentum transport in the disc is mediated by ated with velocity fluctuations in the turbulence, their ef-
fectsaredistinct:diffusiontendstoproduceauniformmix-
ing of material, whereas viscosity does not yield a uniform
distribution of angular momentum. They found that the
⋆ E-mail:[email protected]
† E-mail:[email protected] contaminant behaviour depends on the ratio of the diffu-
‡ E-mail:[email protected] sion coefficient to viscosity; only when this ratio is near
2
or greater than unity does diffusion take place apprecia- gation is similar, except turbulence in our models is driven
bly. Heyvaerts et al. (1996) have shown that in an accre- by a linear instability in the background flow, rather than
tiondiscinwhichtheeffectiveviscosityisduetoturbulence byanassumedforcing.Moreover,theorbitaldynamicsinan
generated by the Balbus-Hawley instability, the associated accretionflowstronglyaffectthevelocityandmagneticfield
magnetic Reynolds number is small, whereas the magnetic correlationsintheturbulence,whichmaydirectlyaffectthe
Prandtl number Pm = ν/η (where η and ν are the mag- turbulentdiffusion.
netic resistivity and the effective viscosity, respectively) is This paperis organised as follows. In Section 2 thenu-
not necessarily so. mericalmethodemployedisdescribed,includingthemethod
Inthiscontext,wewishtoinvestigatetherolethattur- bywhichtheturbulentdiffusioncoefficientiscalculated.Re-
bulent diffusion plays locally in the transport of a passive sults of a variety of simulations are described in Section 3,
tracer, within a small portion of a weakly magnetised ac- while a discussion and conclusions are presented in Section
cretion disc in which the MRI operates. We use numerical 4.
MHDsimulations which adopttheso-called “shearingbox”
model, discussed in detail by Hawley, Gammie & Balbus
(1995, hereafter HGB), a useful scheme for the study of lo-
2 METHOD
cal processes in accretion flows. We extend the numerical
scheme to explicitely follow the time evolution of a passive 2.1 Numerical Method
contaminant in the flow to determine the rate of turbulent
Weuseathree-dimensionalversionoftheZEUScode(Stone
mixing. We use the ideal MHD approximation to simplify
& Norman 1992a;b) in which the shearing box model of
theproblem;thisrequiresthemagneticfieldbewellcoupled
HGB is implemented. Briefly, a coordinate system is cen-
tothegas.Althoughthestructureofprotoplanetarydiscsis
teredatafiducialradiusr0 inasmallpatchofadisc,which
notwellknown,itislikelythatfarfromthecentralstar,the
has linear dimensions much less than r0. This new system
central regions of the disc are dense and cool enough that
corrotates at an angular velocity Ω0 = Ω(r0). The x,y,z
non-idealMHDeffectssuchasOhmicdissipation,ambipolar
coordinates correspond to cilindrical r,ϕ,z coordinates, re-
diffusion,andtheHalleffectwillsignificantlyaltertheMHD
spectively. The resulting MHD equations then incorporate
(Blaes&Balbus1994;Wardle1999;Sanoetal.2000; Stone
the effects of Coriolis forces and tidal stresses (see below).
et al. 2000; Balbus & Terquem 2001; Salmeron & Wardle
The numerical implementation of theshearing box involves
2003). Recent studies of the MRI in non-ideal MHD indi-
boundaryconditionssuchthatthecomputationaldomainis
cate turbulenceand transport my be strongly suppressed if
strictly periodic in the y,z directions at all times, whereas
the magnetic Reynolds number is close to one (Hawley &
in the x (radial) direction it is quasiperiodic: when a fluid
Stone1998; Sano& Miyama1999; Fleming, Stone,& Haw-
elementexitsatoneradialboundary,itreappearsattheop-
ley 2000; Sano & Stone 2002). It is possible that even if
posite radial boundary with a position and a velocity given
theMRIissuppressedinthecentralregionsofthedisc,the
bythe shearing motion.
surfacelayersmaybeionizedbynon-thermaleffectsandre-
Thepassivecontaminant isintroducedinthecodeasa
mainturbulent(Gammie1996;Glassgold etal.1997;Flem-
fractionofthegasdensity.InadditiontotheMHDequations
ing & Stone 2003). Of course, by adopting the ideal MHD
for theshearing box,
approximation, this study cannot be applied to regions of
thediscwheretheionisationfractionissolowthatitissta- ∂ρ
+ (ρv)=0 (1)
ble to the MRI: the turbulent transport of contaminants in ∂t ∇·
these regions will be essentially zero. In regions where the ∂v 1 B2 (B )B 2
MRI is suppressed but not eliminated, the saturation am- +v v= P + + ·∇ 2Ω v+3Ω xxˆ(2)
∂t ·∇ −ρ∇(cid:18) 8π(cid:19) 4πρ − ×
plitude of the turbulence will be reduced, thereby reducing
boththecoefficientofeffectiveviscosityν andtheturbulent ∂B
= (v B) (3)
coefficient D. Thus, theratio of thetwo (theeffective mag- ∂t ∇× ×
neticPrandtlnumberPeff)may beindependentof thelevel ∂ρǫ
of saturation of the instability. Studying diffusion in realis- + (ρǫv)+P v=0 (4)
∂t ∇· ∇·
tic,globalmodelsofprotoplaneterydiscsincludingnon-ideal
MHD and the vertical and radial dependence of the ionisa- P =ρǫ(γ 1) (5)
−
tionfraction isachallengingproblemandbeyondthescope
wherethesymbolshavetheirusualmeanings,thecodealso
of this study.
solves a continuity equation for thecontaminant fraction:
Previousworkhasinvestigateddiffusiveprocessesindif-
∂fρ
ferent magnetohydrodynamic contexts. Yousef et al. (2003) + (fρv)=0 (6)
performed forced turbulence simulations in which they ob- ∂t ∇·
tain a value of the turbulent magnetic Prandtl number of where f is a numberbetween 0 and 1.
aboutunity,regardless ofthevalueofthemicroscopic mag- AsinHGB,weusethestandardvaluesP =P0 =10−6
neticPrandtlnumber.Brandenburgetal.(2004)studynon- (initial pressure), H = 1 (box height), and Ω = 10−3. Our
local aspects of turbulent transport by adding a second- runsadoptthesameparametervaluesasthoseofmodelZ4
order time derivative to the equation for the concentration of HGB (in particular, we use a uniform strength vertical
of a passive contaminant, namely Fick’s law, which states field with plasma parameter β = 400), except that we use
the proportionality of the flux of passive scalar to the neg- aresolution of 88 184 88 grid zones, almost three times
× ×
ative concentration gradient. This extra term results in a greaterineachdimensionthanthatusedinZ4.Theboxsize
damped wave equation for the concentration. Our investi- is 1 2π 1.
× ×
Diffusioncoefficientofa passivecontaminantina localMHD modelofa turbulent accretiondisc 3
2.2 Simulations the fitting is performed independently for each time in the
evolution of thecontaminant.
HGB showed that turbulence develops after a few orbits
Itisworthpointingoutthatweareusingequation(7)as
in the evolution of the shearing box (one orbit=2π/Ω). We
afirst approximation, andit may not describethediffusion
first run the code for 16 orbits, to ensure that turbulence
process of thecontaminant completely.
has set in and is able to stir the contaminant efficiently.
At t=16 orbits the code is restarted and the contaminant
is introduced; we set f = 1 for 0.15 x 0.15 and
− ≤ ≤ 3 RESULTS
for all y and z, forming a slab inside the box. The system
is then allowed to evolve for several more orbits, and the
Figure 1 shows the evolution of the f = 0.5 isosurface of
rateof radial spreading of theslab is used todeterminethe
thecontaminantinRunA1at t=16,16.05 and16.1 orbits
turbulentdiffusion coefficient.
(thatisatthestartoftherun,and0.05and0.1orbitslater).
We expect statistical fluctuations in the evolution of f
At early times, wrinkling of the surfaces on small scales is
to be large given the chaotic nature of turbulence driven
evident,asradialadvectivetransportofthecontaminantby
by the MRI (Winters, Balbus, & Hawley 2003). Thus, to
turbulent eddies occurs. At later times, larger scale distor-
improve our measurement we have performed 22 runs be-
tionsappear,asthecontaminantisadvectedonlargerscales
ginningatdifferenttimes,reintroducingthecontaminantin
bylarger eddies.
its original distribution each time a run is restarted, and
Figure 2 plots the azimuthally and vertically averaged
evolvingtheflowfordifferentperiods.Weshallgroupthese
radial profile of the contaminant from t = 0.05 orbits to
simulations according to the length they were run, label-
t=2.75orbitsafterthecontaminantisintroduced.Theav-
ing each group with different letters, and each run within a
erages are taken from all the data available at each time.
groupbydifferentnumbers.RunA1wasevolvedforoneor-
Overthistimescale, steadydiffusion ofthecontaminantto-
bitbeginningatorbit16.RunsB1throughB8wereevolved
wardsauniformprofile(consistentwiththeradiallyperiodic
sequentially for 0.15 orbits each, starting at orbit 17 (so
boundaryconditions) is evident.
that run B1 provides a measurement of how the contami-
In Figure 3 we plot fits of the model curve, equation
nant mixes in the flow during the evolution from 17.00 to
(8),totheaverageradialprofilesofcontaminantfractionob-
17.15 orbits; run B2 from 17.15 to 17.30 orbits, etc.). Runs
tainedfromthesimulations,atthesamesixtimesasshown
C1throughC13 wereevolvedsequentiallyfor3orbitseach,
inFigure2.ThevalueofthediffusioncoefficientDthatpro-
starting at 18.2 orbits (so that run C1 is evolved from 18.2
videsthebestfitineachcaseisshownatthebottomofeach
to21.2 orbits, etc.). In thisway we areable to measure the
panel.Atotalof57fitswereperformed,correspondingto57
turbulentmixing in theflow over41.2 orbits of evolution.
different instants in the evolution of the contaminant; Fig-
From each run we obtain radial profiles of the az-
ure3plotssixthat arerepresentativeofthefittingprocess.
imuthally and vertically averaged contaminant fraction f
Theerror bars shown correspond to thestandard deviation
at intervals of 0.05 orbits. For each one of these snapshots,
of thesimulation data.
weaveragethecorrespondingradialcurvesoverthenumber
The values of D obtained by fitting equation (8) ev-
of runs performed, and hence obtain a single curve at each
ery 0.05 orbits of evolution are plotted as a function of the
time. For example, to obtain the curve that corresponds to
time after the contaminant is introduced in the turbulent
theradialprofileoff at 0.15orbits,wecanaverageoverall
flowinFigure4.Duringthefirst 0.75orbitsoftheevolu-
22 runs(sinceall were runfor at least that long). However, tion, D grows from approximate∼ly 1.2 10−2H2/orbit to
for the curve that corresponds to the radial profile of f at about 5 10−2H2/orbit. After about ×0.75 orbits, D be-
1.5orbits,weaverageoverthe13runsinthegrouplabelled ×
comes more nearly constant; its value does not fall be-
C. low 4 10−2H2/orbit, and does not exceed 5
10−2≈H2/or×bit, with an average value between 1 an≈d 3 o×r-
2.3 Calculation of diffusion coefficient bitsof (4.7 1.0) 10−2H2/orbit.Asdiscussed in thenext
± ×
section, thechangein behaviorat 0.75 orbitsappears tobe
In order to calculate a value of the diffusion coefficient as- relatedtothespatialextentofthecontaminant,andthesize
sociated with the turbulent flow, we have chosen to solve a ofthelargestturbulenteddiesinthebox.Thediffusiontime
diffusionequationforthecontaminantfractionintheradial tdiff H2/D, using the saturated value of D at late times,
direction with thesame initial condition asthat introduced is 21≡orbits.
in thecode: It is useful to compare the magnitude of D measured
∂f ∂2f from the simulations with the effective viscosity associated
=D , (7)
∂t ∂x2 with shear stress, using the Prandtl number, Peff ν/D.
≡
The valueof ν can becalculated through therelation
1 if x 0.15
f(x,0)=(cid:26) 0 if |x|≤>0.15 ν =αcsH = <<wxy >>H2Ω (9)
| | P
The solution of this initial valueproblem is
where cs is the sound speed, << wxy >> is the time- and
1 x+0.15 x 0.15 volume-averaged value of the r-ϕ component of the stress
f(x,t)= 2(cid:20)erf(cid:16) 2√Dt (cid:17)−erf(cid:16) 2−√Dt (cid:17)(cid:21) (8) tensor (magnetic plus kinetic), H is the disc scale height,
andΩistheangularfrequency.Thesecond equalityfollows
We then fit this solution to the curves obtained from from thefact that
tfihtetinsgimpualaratimonest,ert.reTaotintegstthwehdetiffhuersioDn icsoecffioncsietannttDinatsimthee, α= <<wxy >> (10)
P
4
Figure 1.Evolution ofthe f =0.5isosurfaceof contaminant inRun A,startingat t=16 orbits(left panel)when the contaminant is
introducedintothesaturated turbulence,t=16.05orbits(center panel),andt=16.1orbits(rightpanel).
and ZEUS code. The magnitude and convergence rate of these
errors have been studied in Stone & Norman (1992a). For
cs =HΩ (11)
thesecond-ordervanLeerschemeusedinthesecalculations,
Using a value for << wxy >> which is measured di- discontinuities are diffused until they span 10-15 zones in a
rectly from the simulations (including both the Reynolds few dynamical times, thereafter very little numerical diffu-
and Maxwell stress), we obtain ν 0.52H2/orbit. The sionoccurs.Figure2showstheprofilesreachnearlyuniform
≈
Prandtl numberis then profilesinlessthan3orbits,indicativeofamuchlargerdiffu-
sionratethanthatassociatedwithnumericalerrorsrevealed
ν in the advection tests in Stone & Norman (1992a). Thus,
Peff = D (12) while numerical diffusion may beimportant on small scales
0.52 and at early times, it is unlikely to determine the evolu-
= (13)
0.047 tionobservedinFigure2.Normallytheeffectsofnumerical
11 (14) diffusion would be assesed by repeating the simulations at
≈ different resolutions and looking for convergence. However,
increasingthenumericalresolutioncanalsochangetheprop-
erties of the turbulence driven by the MRI on small scales;
4 DISCUSSION
such resolution studies are less instructive in this case.
Using numerical MHD simulations of turbulence driven by
the MRI, we have measured the diffusion coefficient for a Figure 4 shows that the best fit value of the diffusion
passive contaminant mixed by the flow. When compared coefficientDfoundbyfittingtheanalyticsolution,equation
to the effective viscosity associated with angular momen- (8), to the numerically computed profiles increases nearly
tum transport, the diffusion coefficient is small, that is the linearly in time during the first 0.75 orbits of evolution,
Prandtl numberPeff =ν/D>1. andthereaftersaturatesatanearlyconstantvalueapproxi-
Of course, numerical errors are an important contri- mately four times theinitial value. Thus, theturbulentdif-
bution to the diffusion rate on scales comparable to the fusion coefficient seems to acquire a meaningful, approxi-
grid spacing. However, on larger (resolved) scales up to the matelyconstant valueonly afterthetracerhasbeenspread
box size, there is strong evidence that numerical effects do over large scales, comparable to the box size. The possi-
not determine the diffusion rate. The only numerical errors bility that the diffusion of contaminants in turbulent flows
in solving equation (6) arise in the advection step of the follows a non-Fickian diffusion law has been suggested pre-
Diffusioncoefficientofa passivecontaminantina localMHD modelofa turbulent accretiondisc 5
Figure 2.Evolutionoftheazimuthallyandverticallyaveragedradialprofileofcontaminantfraction.Thecorrespondingtimeafterthe
introductionofthecontaminant intotheturbulentflow(inorbits)isshownaboveeachframe.Thevaluesoff areshownateachradial
gridpointinthesimulation.
viously by Brandenburg et al. (2004). In our case, the vari- ourprimaryresultthatPeff =ν/D>1willnotchange.One
ation of D with time may be related to a variation in the should also be aware of the possibility that this transport
turbulent velocity fluctuations with scale. Because turbu- process follows anomalous diffusion (Bouchaud & Georges
lencedrivenbytheMRIfollowsanearlyKolmogorovscaling 1990).
(HGB), the largest velocity fluctuations are on the largest If the magnetic field acted as a passive scalar in the
scales. As the profile of the contaminant spreads, it sam- flow,andonecouldassociatethediffusioncoefficientDwith
ples larger scales, and larger amplitude eddies. Once the thecoefficientofresistivityη,themagneticPrandtlnumber
profile is spread across the largest scale in the simulation would bePm =ν/η 10, whereν is thetime- and volume-
≈
(the box size), no larger amplitude eddies are possible, and averaged shear stress (Reynoldsplus Maxwell) in thesimu-
D becomes approximately constant. At this stage we have lation. In studies of driven turbulence, Yousef et al. (2003)
not performed simulations with different box sizes; the ef- found that, for weak magnetic fields, Pm is close to unity.
fectofchangingtheboxsizemayalterthevalueofD,since There are a numberof reasons for this discrepancy. Firstly,
theamplitudeoftheMRI turbulenceis set bytheboxsize. themagneticfieldinshearingboxsimulationsoftheMRIis
However, thestress would also bemodified, and the overall clearlynotapassivescalar,thusDmaybeapoorproxyfor
effect would bebest examined via their ratio, Peff. the resistivity. Perhaps just as important, the momentum
The fact that D varies in time for t<0.75 orbits calls transport rate in the shearing box leads to a much larger
into question the applicability of using equation (7) to de- effective shear viscosity than in driven turbulence simula-
scribe the evolution of the contaminant at early times. A tions,becausethecontributionsfromtheMaxwellstressare
more accurate procedure for measuring D would be to use larger than the Reynolds stress (HGB). Through magnetic
the spatial distribution of f at t = 0.75 orbits as the ini- torques, momentum transport can occur indepent of mass
tial condition for equation (7), and solve for the resulting motionthatwouldleadtodiffusionofpassivecontaminants.
evolution, so that a solution is found only for times when The fact that Peff >> 1, so that angular momen-
D is known to be constant. Although this may result in a tum transport is much more efficient than diffusion of con-
slightlymodifiedvalueforDmeasuredfromthesimulations, taminants, has important consequences for the local radial
6
Figure 3.Fitsofthetheoretical curves[equation(8)] tothenumericallymeasuredaverageradialprofilesofthecontaminant shownin
Figure2.Thevalueofthediffusioncoefficientthatprovidesthebestfitisshownatthebottom ofeachpanel,inunitsofH2/orbit.
transportof smallparticlesembeddedinamagnetised disc. pitality, and would like to thank Gordon Ogilvie for fruit-
Global simulations will be required to determine the struc- ful discussions. JS is grateful for financial support from the
ture and amplitude of radial accretion flows driven by the RoyalSociety,theUniversityofCambridge, andNSFgrant
MRI, and therefore assess whether small particles can dif- AST-0413788 for portions of thiswork.
fusive radially outwards in the disc. One of the most im-
portant questions for global simulations of accretion discs
is to determine whether large scale poloidal magnetic fields
can be advected inwards by an accretion flow (leading to REFERENCES
strong vertical fields in the inner regions of the disc), or
BalbusS.A.,HawleyJ.F.1991,ApJ,376,214
whether they can diffuse outwards more rapidly than they
BalbusS.A.,TerquemC.2001,ApJ,552,235
arecarried in(Clarke&Pringle 1988). Ifthediffusion coef- BlaesO.M.,BalbusS.A.1994,ApJ,421,163
ficientDmeasuredherewererepresentativeofthevertically BossA.1998,Annu.Rev.EarthPlanet.Sci.,26,53
averaged resistive diffusion coefficient for poloidal magnetic BouchaudJ.P.,Georges A.1990,Phys.Rep.,195,127
fields, then Pm >> 1 would indicate fields would be ad- Brandenburg A., Ka¨pyla¨ P. J., Mohammed A. 2004, Physics of
vectedinward.However,sincepoloidalmagneticfieldsdrive Fluids,16,1020
theMRI,andsincereconnection(whichisanimportantsat- ClarkeC.J.,PringleJ.E.1988,MNRAS,235,365
urationmechanisminthenonlinearregimeoftheMRI)can FlemingT.,StoneJ.M.,HawleyJ.F.2000,ApJ,530,464
FlemingT.,StoneJ.M.,2003,ApJ,585,908
breakflux-freezing,itisunlikelyequation(6)canbeusedto
GammieC.F.1996,ApJ,457,355
describetheradialevolutionofpoloidalfields.Globalsimu-
GlassgoldA.E.,NajitaJ.,IgeaJ.1997, ApJ,480,344
lationsofthinKepleriandiscswhichcanaddresstheseissues
Hawley J. F., Gammie C. F., Balbus S. A. 1995, ApJ, 440, 742
are needed.
(HGB)
HawleyJ.F.,StoneJ.M..1998,ApJ,501,758
Acknowledgements: AC acknowledges support from HeyvaertsJ.,PriestE.R.,BardouA.1996,ApJ,473,403
CONACYT scholarship 167912, thanks the Department of KellerCh.,GailH.P.2004,A&A,415,1177
AstrophysicalSciencesatPrincetonUniversityfortheirhos- MorfillG.E.1983,Icarus,53,41
Diffusioncoefficientofa passivecontaminantina localMHD modelofa turbulent accretiondisc 7
Figure 4.Turbulentdiffusioncoefficient D,obtainedfromthecurve-fittingproceduredescribedinthetext, asafunctionoftime.
SalmeronR.,WardleM.,2003,MNRAS345,451.
SanoT.,MiyamaS.M.,1999,ApJ515,776
Sano T., Miyama S.M., Umebayashi T., Nakano T., 2000, ApJ
543,486
SanoT.,StoneJ.M.,2002, ApJ,570,314
SreenivasanK.R.,1991,Proc.R.Soc.London.Ser.A434,165
StoneJ.M.,NormanM.L.1992a, ApJS,80,753
StoneJ.M.,NormanM.L.1992b, ApJS,80,791
Stone J.M.,GammieC.F,Balbus S.A.,HawleyJ.F.,2000, Pro-
tostars and Planets IV, ed by V. Mannings, A. Boss, & S.
Russell,UAPress,Tucson,p589.
Takeuchi T.,LinD.N.C.2002,ApJ,581,1344
YousefT.A.,BrandenburgA.,Ru¨digerG.2003,A&A,411,321
WarhaftZ.,2000,Annu.Rev.FluidMech.,32,203
WardleM.,1999, MNRAS,307,849
WintersW.,BalbusS.A.,HawleyJ.F.,2003,MNRAS,340,519