Turbulence In the Outer Regions of Protoplanetary Disks. I. Weak Accretion with No Vertical Magnetic Flux

1.4 MB·English

by  Jacob B. Simon

Draftversion October 17,2012
PreprinttypesetusingLATEXstyleemulateapjv.5/2/11
TURBULENCE IN THE OUTER REGIONS OF PROTOPLANETARY DISKS.
I. WEAK ACCRETION WITH NO VERTICAL MAGNETIC FLUX
Jacob B. Simon1, Xue-Ning Bai2,3,4, James M. Stone2, Philip J. Armitage1,5, and Kris Beckwith1,6
Draft versionOctober 17, 2012
ABSTRACT
Weuselocalnumericalsimulationstoinvestigatethestrengthandnatureofmagnetohydrodynamic
2
(MHD)turbulenceintheouterregionsofprotoplanetarydisks,whereambipolardiffusionisthedomi-
1
nantnon-idealMHDeffect. Thesimulationsincludeverticalstratificationandassumezeronetvertical
0
magnetic flux. We employ a super time-stepping technique to ameliorate the Courant restriction on
2
the diffusive time step. We find that in idealized stratified simulations, with a spatially constant am-

t bipolar Elsasser number Am, turbulence driven by the magnetorotational instability (MRI) behaves
c
in a similar manner as in prior unstratified calculations. Turbulence dies away for Am 1, and
O
≤
becomes progressively more vigorous as ambipolar diffusion is decreased. Near-ideal MHD behavior
5 is recovered for Am 103. In the intermediate regime (10 Am 103) ambipolar diffusion leads
≥ ≤ ≤
1 to substantial increases in both the period of the MRI dynamo cycle and the characteristic scales
of magnetic field structures. To quantify the impact of ambipolar physics on disk accretion, we run

] simulationsat30AUand100AUthatincludeaverticalAmprofilederivedfromfarultraviolet(FUV)
R
ionizeddiskmodels. ThesemodelsdevelopaverticallylayeredstructureanalogoustotheOhmicdead
S zone that is present at smaller radii. We find that, although the levels of surface turbulence can be
. strong(andconsistentwith constraintsonturbulent line widths atthese radii),the inferredaccretion
h
rates are at least an order of magnitude smaller than those observed in T Tauri stars. We speculate
p
that this discrepancy may be due to the neglect of vertical magnetic fields. If this is the case, then
-
o theMRI alonemayresultindisjointclassesofdiskevolution,withthosediskslackingnetfieldsbeing
r weakly viscous at large radii.
t
s Subject headings: accretion, accretion disks — (magnetohydrodynamics:) MHD — turbulence —
a protoplanetary disks
[

1
v 1. INTRODUCTION particles. Once planetesimals haveformed,gravitational
coupling to turbulent fluctuations in the disk may af-
4 The structure and evolution of protoplanetary disks
6 play a crucial role in the formation of stars and their fect their growth (Ida et al. 2008). Finally, the strength
1 planetary systems. Disk gas is observed to accrete onto andnatureofturbulencedetermineswhetherthe critical
4 the central star at rates that require some form of an- co-orbital contribution to the Type I migration torque
0. gular momentum transport substantially stronger than remains unsaturated (Paardekooper et al. 2011).
Ataminimum,turbulenceinprotoplanetarydiskswill
1 that provided by molecular viscosity. Turbulence has
begeneratedinregionswherethemagnetorationalinsta-
2 long been suggested as the source of enhanced trans-
bility (MRI; Balbus & Hawley 1998) operates. Indeed,
1 port (Shakura & Syunyaev 1973). This turbulence not
simulations of the non-linear evolution of the MRI un-
: only allows accretion, but can also play an important
v deridealmagnetohydrodynamic(MHD)conditionsyield
role in the formation and subsequent evolution of plan-
i sustainedturbulencethattransportsangularmomentum
X ets. At early times, turbulence could act to inhibit
outwardat rates in generalagreementwith observations
dust settling and largely determine the collisional ve-
r (Hartmann et al.1998). However,largeregionsofproto-
a locities that affect the balance between fragmentation
planetary disks are expected to have very low ionization
and coagulationof these particles (Ormel & Cuzzi 2007;
fractions (e.g., Ilgner & Nelson 2006), which in turn re-
Youdin & Lithwick 2007; Birnstiel et al. 2011). Per-
sult in three significant non-ideal MHD effects: Ohmic
sistent pressure maxima predicted by some turbulence
diffusion, ambipolar diffusion, and the Hall effect (see,
models (Barge & Sommeria 1995; Johansen et al. 2009;
e.g., Armitage 2011). The relative importance of these
Uribe et al. 2011; Simon et al. 2012) may act to concen-
effects depends primarily upon the density (as well as
trate particles, enhancing their coagulation into larger
magnetic field strength). Ohmic diffusion is efficient at
[email protected] high densities, and is thus most important in the inner
1JILA,UniversityofColoradoandNIST,440UCB,Boulder, regions of the disk (outside a small zone very close to
CO80309-0440
the star where thermal ionization of alkali metals pro-
2DepartmentofAstrophysicalSciences,PrincetonUniversity,
Princeton,NJ08544 vides sufficient ionization throughout the disk column).
3HubbleFellow Atlowgasdensities,suchasintheouterdisk,ambipolar
4Current address: Harvard-Smithsonian Center for Astro- diffusion becomes dominant, while at intermediate den-
physics,60GardenSt.,MS-51,Cambridge,MA02138
5Department of Astrophysical and Planetary Sciences, Uni- sities, the Hall term is important (e.g., Kunz & Balbus
versityofColorado,Boulder,CO80309 2004).
6Tech-XCorporation,5621ArapahoeAve.,SuiteA,Boulder, MRI physics and the phenomenological consequences
CO80303
2
of the non-ideal terms have been best-characterized in the collision frequency between neutrals and ions ex-
the case of Ohmic diffusion. The evolution in this ceeds the angular frequency. Simulations in the alter-
limit depends upon the Elsasser number, defined as nate regime,where the recombinationtimescale for elec-
Λ v2 /ηΩ, where v is the Alfv´en velocity in the trons is assumed to be very long, were conducted by
≡ A,z A,z
verticaldirection,η istheOhmicresistivity,andΩisthe Hawley & Stone (1998) using a two-fluid approach in
angularfrequencyof Keplerianrotation. ForΛ less than whichtheionsandneutralswereevolvedseparately,only
order unity, MRI turbulence is severely quenched (Jin interactingthroughcollisions. Theprimaryresultofthis
1996;Turner et al.2007),whileforlargervaluestheMRI work was that for a collision frequency less than 0.01Ω,
saturationlevelwilldependonthestrengthofOhmicdif- the ions and neutrals behave independently, but for a
fusion; stronger diffusion leads to lower turbulence lev- frequency larger than 100Ω, the gas behaves as though
els. Combining these results with chemical models for its fully ionized. When both the collision and orbital
disks motivates the dead zone model of disk accretion frequencies are comparable, the saturation level of the
(Gammie 1996). In this model, the disk is well-ionized turbulence is primarily controlled by the ion density.
only in its surface layersdue to non-thermal sources (X- The importance of ambipolar diffusion for outer pro-
rays, cosmic rays, and far ultraviolet (FUV) photons) toplanetary disks, the advent of well-resolved mm-wave
that penetrate the disk from the exterior down to some observations of this region, and advances in numerical
column depth. Closer to the mid-plane, the ionization techniques,allmotivatemoredetailedstudiesofhowam-
fractionislow,resultinginasmallΛ,andnoMRI-driven bipolar diffusion affects the saturated state of the MRI.
turbulence (e.g., Gammie 1996; Fleming & Stone 2003; Bai & Stone (2011)carriedoutshearingbox simulations
Turner et al. 2007; Turner & Sano 2008; Oishi & Low in the strong-coupling limit to determine how the MRI
2009). Thus, MRI-driven accretion occurs in active lay- saturation level correlates with dimensionless number,
ers only, leaving much of the disk mass near the mid- Am, defined as the frequency for the neutrals to collide
plane magnetically inactive. withtheionsdividedbytheorbitalfrequency(see 2.1).
§
Qualitative changes to the predicted disk structure Theyfoundthatforsimulationswithanettoroidalmag-
mayequallyresultfromambipolardiffusionandtheHall netic flux, but no vertical magnetic flux, the turbulence
effect, though for these terms the understanding of the dies away for Am 1, in line with previous studies.
≤
non-linear behavior is incomplete. The linear regime of For sustained turbulence runs, the saturated turbulent
the MRI in the presence of the Hall term has been ex- stresses increase with increasing Am, eventually asymp-
ploredbyWardle(1999),Balbus & Terquem(2001),and toting towards the ideal MHD level. However, in the
Wardle & Salmeron(2012). Aprimaryresultfromthese presence of a net vertical magnetic flux, turbulence can
studies is that the growth rate of the MRI is strongly alwaysbe sustainedevenfor Am<1, assumingthat the
affected by the sign of Ω B, i.e., how the vertical mag- background vertical magnetic flux is weak enough. For
netic field is aligned with· the angular velocity vector. low Am values, however, the resulting turbulence levels
The only study of the non-linear, turbulent state of the are fairly small.
MRI in the presence of the Hall term was carried out in Following the above lines of investigation, we study
Sano & Stone (2002a)andSano & Stone (2002b). Their the effect of ambipolar diffusion on the MRI in the
numericalsimulationsincludedbothOhmicdiffusionand outer region of protoplanetary disks by performing nu-
the Hall term, with the Ohmic contribution dominating merical simulations using a more realistic disk structure
significantlyoverthe Halleffect. Inthis regime,the Hall than attempted previously. We include vertical strat-
term does not strongly influence the saturated state of ification (absent in all prior work except for that of
the MRI. However, the regime in which the Hall term Brandenburg et al. (1995)), which has been shown in
dominates has yet to be explored through simulations some previous MRI calculations to lead to significant
(Wardle & Salmeron 2012). qualitative changes in the non-linear evolution. For ex-
Ambipolardiffusionarisesfromtheimperfectcoupling ample,Davis et al.(2010)showedthatintheidealMHD
between ionized species and neutrals. The linear anal- case, the turbulence properties converge with numer-
yses of Blaes & Balbus (1994), Kunz & Balbus (2004), ical resolution, while in unstratified simulations (and
and Desch (2004) showed that the growth of the MRI is for a vertical domain size of one scale height or less,
damped when the collision frequency between the neu- (Stone, private communication)), the stress level de-
trals and the ions is smaller than the orbital frequency. creases with resolution with no signs of convergence
This is intuitive; neutrals need to communicate with the (Fromang & Papaloizou 2007). Similarly, Simon et al.
ions faster than the timescale over which the MRI acts (2011b)showedthatinthepresenceofOhmicresistivity,
(i.e., the dynamical one) in order for the neutrals to feel vertical gravity can lead to large amplitude fluctuations
any MRI-like effect at all. in the stress levels, a behavior that is absent without
Initial two and three-dimensional simulations of non- vertical gravity.
linear MRI turbulence in the presence of ambipolar Wealsoaimtotranslateouridealizedunderstandingof
diffusion were carried out by Low et al. (1995) and the ambipolar-dominated MRI into predictions for tur-
Brandenburg et al. (1995), respectively. These authors bulenceandaccretioninthe outerregionsofprotoplane-
consideredthesingle-fluid,“strong-coupling”limit,valid tary disks. We run simulations that mimic realistic con-
when the recombination timescale is much shorter than ditionsintheouterdisk(whichreflectthechemistrycal-
the dynamical time. This limit is generally applicable culations in Bai (2011a) and the FUV ionization model
to protoplanetary disks (Bai & Stone 2011; Bai 2011a). of Perez-Becker& Chiang (2011b)). These simulations
Their results agreed with the expectations of linear the- will provide a quantitative measure of the turbulence,
ory (Blaes & Balbus 1994); the MRI only operates if structure, and evolution of the outer disk regions.
Our investigation is divided into two sets of studies,
3
∂ρv 1
using different magnetic configurations. The first, which + (ρvv BB)+ P + B2
wepursuehere(PaperI),assumeszeronetverticalmag- ∂t ∇· − ∇(cid:18) 2 (cid:19) (2)
netic flux7. Although it is not the most geometrically =2qρΩ2x ρΩ2z 2Ω ρv
favorable configuration for the ambipolar-dominated − − ×
MRI (Bai & Stone 2011), this is the most studied field ∂B (J B) B
(v B)= × × , (3)
configuration in the literature (e.g., Stone et al. 1996; ∂t −∇× × ∇×(cid:20) γρ ρ (cid:21)
i
Fleming & Stone 2003), and it allows us to make direct
comparisons between the ambipolar MRI and the ideal where ρ is the mass density, ρv is the momentum den-
MHD case (Davis et al. 2010), as well as the case with sity, B is the magnetic field, P is the gas pressure, and
Ohmicresistivity(Simon et al.2011b). Itisalsoentirely q is the shear parameter, defined as q = dlnΩ/dlnR.
−
conceivable (though, maybe not particularly likely) that We use q = 3/2, appropriate for a Keplerian disk. We
there will be some regions of protoplanetary disks that assume an isothermal equation of state P = ρc2, where
s
have negligible vertical magnetic fields; our results will c is the isothermal sound speed. From left to right, the
s
apply to such regions. The second set of studies will in- source terms in equation (2) correspond to radial tidal
cludeanon-zeroverticalnetmagneticflux,andwedefer forces (gravityand centrifugal), verticalgravity,and the
that problem to Paper II. Coriolis force. The source term in equation (3) is the
The structure of the paper is as follows. In Section 2, effect of ambipolardiffusion on the magnetic field evolu-
we describe our equations, the methods used to solve tion,whereρ istheiondensity,andγ isthecoefficientof
i
them and the initial conditions for our simulations. In momentum transfer for ion-neutral collisions. Note that
Section 3.1, we study the effect of Am (here assumed oursystemofunitshasthemagneticpermeabilityµ=1,
constant in space and time) on vertically stratified MRI and the current density is
turbulence. Then, in Section 3.2, we apply a realistic
protoplanetary disk model to allow for a spatially and J = B. (4)
temporally varying Am. Section 4 discusses the implica- ∇×
tions of our results for real protoplanetarydisks, and we Adopting this shearing box approximation allows for
wrap up with conclusions in Section 5. better resolution of small scales within the disk, at
the expense of excluding global effects (those of scale
R ). These scales could be physically important
0
2. METHOD (∼Sorathia et al. 2012; Simon et al. 2012). However, the
2.1. Numerical Method trade-offisworthwhileforourpurposes,becauseweneed
to study not only models where ambipolar diffusion is
In this study, we use Athena, a second-order accurate
dominant, but also situations where diffusion is only im-
Godunovflux-conservativecodeforsolvingtheequations
portant on small scales.
of MHD. Athena uses the dimensionally unsplit corner
The numerical integration of the shearing box equa-
transport upwind (CTU) method of Colella (1990) cou-
tions require additions to the Athena algorithm, the de-
pled with the third-order in space piecewise parabolic
tails of which can be found in Stone & Gardiner (2010)
method (PPM) of Colella & Woodward (1984) and a
and the Appendix of Simon et al. (2011b). Briefly,
constrained transport (CT; Evans & Hawley 1988) al-
Crank-Nicholsondifferencingisusedtoconserveepicyclic
gorithm for preserving the B = 0 constraint. We
∇ · motion exactly and orbital advection to subtract off the
use the HLLD Riemann solver to calculate the nu-
background shear flow (Stone & Gardiner 2010). The y
merical fluxes (Miyoshi & Kusano 2005; Mignone 2007).
boundary conditions are strictly periodic, whereas the
A detailed description of the base Athena algorithm
x boundaries are shearing periodic (Hawley et al. 1995).
and the results of various test problems are given in
The vertical boundaries are the outflow boundary con-
Gardiner & Stone (2005), Gardiner & Stone (2008), and
ditions described in Simon et al. (2011b). The electro-
Stone et al. (2008).
motive forces (EMFs) at radial boundaries are properly
Oursetupisspecializedandnecessarilymorecomplex
remapped which guarantees that the net vertical mag-
than the base algorithm. First, our simulations utilize
neticfluxisstrictlyconservedtomachineprecisionusing
the shearing box approximation, which is a model for
CT (Stone & Gardiner 2010).
a local, co-rotating disk patch whose size is small com-
Theintegrationoftheambipolardiffusiontermalsore-
pared to the radialdistance from the centralobject, R .
0 quires some modifications to the algorithm. Ambipolar
This allows the construction of a local Cartesian frame
diffusion is implemented in a first-order operator-split
(x,y,z) that is defined in terms of the disk’s cylindrical
manner as in Bai & Stone (2011); the ambipolar diffu-
co-ordinates (R,φ,z′) via x = (R R ), y = R φ, and
z =z′. Thelocalpatchco-rotatesw−ith0anangula0rveloc- sion term is integrated separately from the ideal MHD
integrator. Furthermore,asisevidentfromequation(3),
ity Ω corresponding to the orbital frequency at R , the
0 the ambipolardiffusion termcan be written as an EMF.
center of the box; see Hawley et al. (1995). Thus, the
Thus integrating it is done via the CT method to pre-
equations to solve are: serve B =0. The ambipolar diffusion EMFs are also
∇·
remappedattheradialboundariesinthesamewayasthe
∂ρ ideal MHD EMFs in order to conserve the vertical mag-
+ (ρv)=0, (1)
netic net flux. In addition, before even remapping these
∂t ∇·
ambipolar diffusion EMFs at the radial boundaries, we
7 Of course due to the dynamo action of the MRI(Davisetal. also remap the toroidal current densities Jy (located at
2010; Simonetal. 2011b), the net radial and toroidal fields are cell edges) so that the line integral Jydy along the in-
allowedtoevolveintime. ner and outer radial boundaries areRequal. We find that
4
this procedure is essential to avoid spurious numerical We note that as ν 0, ∆t N2∆t so that
STS diff
features at shearing-box boundaries.8 the STS approachis a→symptotically→N times faster than
Throughout this paper, the strength of ambipolar dif- the standard explicit approach. However, the value of ν
fusionwillcharacterizedbythe ambipolarElsassernum- needs to be properly chosen to achieve the optimal bal-
ber ance between performance and accuracy. In general, the
γρ
i STSmethodprovidesbetteraccuracyasN decreasesand
Am , (5)
≡ Ω ν increases, whereas large N and small ν lead to higher
efficiency. Here, we choose ν = 1/4N2 with a limit of
which corresponds to the number of times a neutral
N 12. At N = 12, one achieves an acceleration fac-
molecule collides with the ions in a dynamical time ≤
(Ω−1). Am can be rewritten as tor of about 9. It is also found that further increasing
N wouldnotsignificantly increasethe efficiency without
sacrificingaccuracy(Stone,privatecommunicationbased
v2
Am= A , (6) on a Princeton Junior Project done by Sara Wellons).
ηAΩ In our calculations, we first compute the ideal MHD
time step ∆t and the diffusion timestep ∆t . The
which is a form reminiscent of the Ohmic Elsasser num- MHD diff
number of super time steps N can be found from the
ber. As shown by equation (5), Am is independent of
condition G[N 1,1/4(N 1)2] < ∆t /∆t
the Alfv´en speed; this comes about because the ambipo- − − MHD diff ≤
G(N,1/4N2). If N 12, then we modify ∆t so
lar diffusivity, η is defined as diff
A that ∆t G(N,≤1/4N2)∆t . Otherwise, we fix
MHD diff
v2 N = 12, and≡set ∆t = ∆t G[12,1/(4 122)].
ηA ≡ γρAi. (7) In this way, we alwaMysHDhave ∆tdMiffHD = ∆tST×S. As
we use an operator-split algorithm for magnetic diffu-
Thisdiffusivityisresponsiblefordeterminingthediffu- sion, we integrate N STS substeps of the ambipolar dif-
∆sivxe2/tiηmAe. stSeipncien athCisoudriaffnutsilvimityiteids cparlocpuolarttiioonn;al∆ttodiffth∝e fsutespionwittehrm∆wtSiTthS.∆τWj itbhefoSrTeSe,vowlevinhgavoenereMpeHatDedtitmhee
Alfv´enspeedsquared,itcanbecomeverylargeintheup- test problems (i.e., linear wave damping test and C-
per disk regions, making the Courant limited time step type shock test) performed in Bai & Stone (2011) and
extremely small in some of our calculations. found essentially the same results for N up to 10.
STS
To circumvent this issue, we have implemented the Moreover,we repeated the unstratified MRI simulations
super time-stepping (STS) technique of Alexiades et al. with Am = 1 (runs Z5 and M5 in Bai & Stone (2011))
(1996)toaccelerateourcalculations. TheSTStechnique and with STS turned on. In these simulations N
STS
has already been successfully implemented and tested reaches 12, and the stress level we find is the same as
for studying ambipolar diffusion in multi-fluid codes by reported in Bai & Stone (2011). Combining our tests
O’Sullivan & Downes (2006) and O’Sullivan & Downes withthe successfultests ofO’Sullivan & Downes(2006),
(2007) and in a single-fluid code by Choi et al. (2009). O’Sullivan & Downes (2007) and Choi et al. (2009), we
Ourimplementationissimilartotheirs,aswebrieflyde- are confident that the STS technique implemented here
scribe below. is capable of achieving substantial speed-up while main-
The STS method divides a compound timestep ∆tSTS taining accuracy.
into N unequal substeps ∆τ (j =1,...,N) with
j
2.2. Am Profiles
N
∆t = ∆τ . (8) For most of our simulations,we fix Am to be constant
STS j
Xj=1 inorderto study the effect ofambipolardiffusionon the
non-linearsaturationoftheMRIinthepresenceofverti-
Bychoosing∆τ judiciously,stabilitycanbe maintained cal stratification. This is designed to be the next logical
j
evenwhentheaveragedtimestep∆t /N ismuchlarger step in extending the work of Bai & Stone (2011) where
STS
than the normal stable diffusion timestep ∆t . The vertical stratification was not included. We have con-
diff
optimized lengths for the substeps were found to be sidered Am=1,10,100,300,103 and 104. Although this
(Alexiades et al. 1996) prescription of a constant Am profile is highly simpli-
fied, it is a necessary, incremental step between the con-
−1
2j 1π stant Am models without vertical gravity (Bai & Stone
∆τ =∆t (ν 1)cos − +ν+1 , (9)
j diff(cid:20) − (cid:18) N 2(cid:19) (cid:21) 2011) and simulations that incorporate a more realistic
prescriptionfor ambipolardiffusion, whichwe also carry
where 0 < ν < 1 is a free parameter. The sum of the out (see below).
substeps gives TheresultsofBai(2011a)showthatthephysicalvalue
ofAmintheouterregionsofPPDsisoforderunity. Re-
N (1+√ν)2N (1 √ν)2N cently, Perez-Becker& Chiang (2011b) pointed out that
∆t =∆t − −
STS diff2√ν(cid:20)(1+√ν)2N +(1 √ν)2N(cid:21) (10) the surfacelayerofprotoplanetarydisks shouldbe much
−
better ionized due to far ultraviolet (FUV) photon ion-
G(N,ν)∆t .
diff
≡ ization from the central star; these photons almost fully
8 We note that strict magnetic flux conservation (remap of the ionize the carbon and sulfur to overcome the effects of
EMFs)wasnotenforcedinBai&Stone(2011),inwhichcasethis
recombination onto dust grains. Their results suggest a
aatdidointsioinnalverretmicaaplnoeftJmyawgansetnicotflnuexceinssaBrayi.&NeSvteorntehe(l2e0s1s1,)thseimvualrai-- largeionizationfraction(f 10−5)downtoasmallpen-
tionsweretiny(lessthan0.01%),whichdidnotaffecttheirresults. etration depth of Σ ∼0.01 0.1 g cm−2, relatively
FUV
∼ −
5
independent of disk radius. Such an ionization fraction
shouldsignificantlyreducethestrengthofambipolardif-
fusion (i.e., increase Am) in the disk surface layers.
Thus, in our second set of simulations, we include the
effectofFUVionizationatthedisksurfacelayerstogive
Am a more realistic spatial and temporal dependence.
Since we are not including Ohmic resistivity (Gammie
1996) in our calculations, these particular models are
only appropriate for the outer regions of protoplanetary
disks (e.g., beyond 10 AU). We adopt a minimum-
∼
mass solar nebular disk model with surface density of
Σ = 1700R−3/2g cm−2 (Weidenschilling 1977; Hayashi
AU
1981), where R is the disk radius measured in AU.
AU
We canexpressthe value of Am within the FUV ionized
layer as follows (Bai & Stone 2012)
Am 3.3 107 f ρ R−5/4 , (11)
FUV ≈ × (cid:18)10−5(cid:19)(cid:18)ρ (cid:19) AU
0,mid Figure 1. Vertical profile for Am at RAU = 30. The profile
corresponds to the initial time of Z30AU, which is orbit 22 from
where f is the ionization fraction and ρ0,mid is the mid- the Am = 105 run with the same domain size. The value of Am
planedensity. Forsimplicity,wefixf =10−5,ρ0,mid =1, has been averaged horizontally. The units of the x axis are the
and assume a penetration depth of Σ = 0.1g cm−2 verticalscaleheight,H. Theasterisksdenotethelocationsofgrid
(which is slightly different from thaFtUiVn Bai & Stone zones. AmtransitionsfromAm=1toAm=AmFUV (asdefined
in the text), and this transition is smoothed over roughly 7 grid
2012). We conduct two simulations that correspond to zonesusingtheerrorfunction.
radiallocations atR=30AU and R=100AU. Assum-
ingthedensityprofileisGaussian(seeEquation15),one the run at R = 30 at the initial time, referring here
AU
finds that the base of the FUV layer (at which the col- to when the run was restarted from C1e5 (see below).
umn density equals Σ ) is located at z = 1.7H for The asterisks denote the grid cell centers; as previously
FUV b
R = 30 AU and z = 1.1H for R = 100 AU (H is the mentioned,thetransitionregionisresolvedby 7zones.
b
∼
vertical scale height as defined in equation (16) below).
In our simulations,we set Am=1 for z <z <z as a 2.3. Simulations
b b
−
proxy based on the calculations of Bai (2011a), and use
Wehaverunsimulationswithseveraldomainsizesand
Equation(11)fortheionizedsurfacelayersthedisk. For Am profiles. All of the simulations with Am < 105 or
simplicity, we keep the value of z fixed throughout the
b withspatially andtemporallyvaryingAm areinitialized
calculation.
from the turbulent state of a “starter” calculation with
From these considerations, the value of Am changes thesamedomainsizebutwithAm=105(i.e.,reasonably
quite dramatically from Am=1 to Am 3 104 at the
close to ideal MHD).
≈ ×
base of the FUV layer. This very large transition is
Thesestartersimulationsareinitializedwithadensity
smoothed over roughly 7 grid zones so as to prevent a
corresponding to isothermal hydrostatic equilibrium.
discontinuous transition in Am. The smoothing func-
tionsweapplyarebaseduponthe errorfunction(ERF). z2
Thus, the complete profile of Am for these runs is given ρ(x,y,z)=ρ0exp(cid:18)−H2(cid:19), (15)
by
whereρ =1isthemid-planedensity,andH isthescale
0
height in the disk,
AmFUV z≥zb+∆z
Am≡111++ 1122AAmmFFUUVVSS+−((zz)) −−zbzz−bb−+n∆∆n∆zzz<<≤zzz<<≤z−bzz+bb−+∆nzn∆∆zz H = √Ω2cs. (16)
AmFUV z≤−zb−∆z (12) The isothermal sound speed, cs = 7.07× 10−4, corre-
where S+(z) and S−(z) are the smoothing functions de- sponding to an initial value for the mid-plane gas pres-
sure of P = 5 10−7. With Ω = 0.001, the value for
fined as the scale0height×is H = 1. A density floor of 10−4 is
z−0.9z applied to the physical domain as too small a density
S+(z)≡1+ERF b , (13) leadstoa largeAlfv´enspeedandaverysmalltime step.
(cid:18) ∆z (cid:19)
Furthermore,numericalerrorsmake it difficult to evolve
S−(z)≡1−ERF z+0.9zb , (14) regionsofverysmallplasmaβ (ratioofthermalpressure
(cid:18) ∆z (cid:19) to magnetic pressure).
Here, n = 8 and ∆z = 0.05H. These numbers were The initial magnetic field is purely toroidal and has a
chosen to give a reasonably resolved transition region constantβ =100throughoutthedomain(thus,By2 hasa
betweenAm=1andAm . Foravisualrepresentation Gaussianshapelikethedensity). Randomperturbations
FUV
of the rather complex equation (12), we plot in Fig. 1 areaddedtothedensityandvelocitycomponentstoseed
the vertical profile of Am (averaged over x and y) for the MRI.
6
The remaining simulations are restartedfrom orbit 50
(orbit22forthesimulationswithdomainsize8H 16H Table 1
× × ShearingBoxSimulations
8H9)oftheir correspondingdomainsizesimulationwith
Am = 105. They are listed in Table 1. The label for
Label AmbipolarDiffusion DomainSize α
each calculation describes whether the value of Am is (Lx×Ly×Lz)H
constant with height, labeled C, or varies according to
C1 Am=1,constant 4×8×8 decayed
equation (12), labeled Z. For the constant Am simula-
C10 Am=10,constant 4×8×8 0.0046
tions,thenumberfollowingtheCisthevalueofAm. For C10L Am=10,constant 8×16×8 0.0055
the spatially varyingAm calculations,the number after- C100S Am=100,constant 2×4×8 0.070
C100 Am=100,constant 4×8×8 0.0062
wards (along with the AU) describes the radial location
C300 Am=300,constant 4×8×8 0.024
of the shearing box in our protoplanetary disk model in C1000 Am=1000, constant 4×8×8 0.022
unitsofAU.AnS(L)followingtheAmvaluecorresponds C10000 Am=104,constant 4×8×8 0.025
to a domain size of 2H 4H 8H (8H 16H 8H), C1e5 Am=105,constant 4×8×8 0.038
which is smaller (larger)×than×the 4H ×8H 8×H size Z30AU AmatR=30AU 8×16×8 0.0016
× × Z100AU AmatR=100AU 8×16×8 0.0015
of most of the constant Am calculations. The “starter”
simulationforthe4H 8H 8H runsisalsoincludedin
× ×
the table, labeled C1e5. Finally, all of our calculations value could not be examined because even with STS,
are carried out at a resolution of 36 grid zones per H. the diffusion limited time step is very small for Am =
1; running it out further would be very computation-
3. RESULTS
ally expensive. However, these results strongly indicate
3.1. Constant Am Calculations that the MRI turbulence has completely decayed away,
consistent with the results of Bai & Stone (2011). This
We begin by applying some standard diagnostics to
behavior will play an important role in the variable Am
the set of calculations with constant values of Am. The
simulations of Section 3.2.
firstsuchdiagnosticisthedensity-weightedMaxwelland
We time-average this normalized stress and define the
Reynoldsstresses(seeequation(37)ofBalbus & Hawley
Shakura-Sunyaevα parameter,
1998), defined as
W = hρvxδvy−BxByi (17) α= WRφ (18)
Rφ ρ c2s
h i
where the angled brackets denote a volume averageover where the overbar denotes the time average, which is
the whole domain. Figure 2 shows the time evolution of done from orbit 100 onwards for most of the constant
this total stress for the 4H 8H 8H runs, normalized Am runs with Am > 1; from orbit 72 onwards for runs
by the square of the sound×speed.×The dashed line indi- C10L, Z30AU, and Z100AU; and from orbit 25 to 53
cates the averagedvalue (from orbit25to 53)ofthe run for the Am = 105 “starter” simulation. Fig. 3 displays
withAm=105fromwhichalloftheother4H 8H 8H α versus Am for these runs. The arrow on the Am =
runswererestarted. ThedifferentAmvaluesa×reden×oted 1 run indicates that that the stress level is continually
by the color. The runs with Am > 1 appear to adjust decreasing. The trend of α with Am can be compared
on a roughly 50 orbit timescale after which a statisti- with the unstratified simulations shown in Fig. 10 of
cal steady state follows. In general, the stress increases Bai & Stone(2011). Thesetrendsroughlyagree,though
withincreasingAm(decreasingdiffusion)forthese runs. the results from Bai & Stone (2011) show a monotonic
However,theAm=10andAm=100runshaveroughly increase of α with Am, whereas our results show that
thesamevaluesatlatetimes,asdotheAm=300,1000, different Am values can lead to very similar α values.
and 10000 runs. This difference may be attributable to different back-
The Am = 1 case deserves extra attention. From ground magnetic field strengths as the background field
Fig. 2, it would appear that the turbulence completely evolves via the usual MRI dynamo (e.g., Davis et al.
dies away. A closer examination of the stress histo- 2010; Simon et al. 2011b). To understand this further,
ries show that the Maxwell stress levels out to a small, first let us consider the space time diagrams of the
but positive value, while continuing to slowly decrease toroidal field By for C10, C100, C300, and C1000 as
with time. The Reynolds stress approaches an oscilla- shown in Fig. 4. In these diagrams, the field has been
tory behavior which occasionally brings it below zero. averaged over x and y and is plotted in the (t,z) plane.
Space-time plots of various quantities in this run indi- The most obvious feature from these diagrams is that
cate that the gas is no longer MRI turbulent. The rem- the period of the dynamo flipping of By changes with
nantMaxwellstress is the result ofa residuallargescale Am;asambipolardiffusionbecomesstronger,theperiod
B and B field near the mid-plane, and the Reynolds increases. In particular,for Am =10,the periodis 50
x y ∼
stress appears to arise from residual waves propagating orbits,andforAm=100(onlyconsideringtimespastthe
through the box. The longer term behavior of this Am initial 50 orbittransientperiod), the period is 15 20
∼ −
orbits. For Am 300, the period is 10 orbits as is
≥ ∼
9Wechooseadifferentrestarttimeforthesecalculationsbecause usually observed in stratified MRI simulations.
wedecidedmidwaythroughrunningoursimulationsthatalarger The most relevant feature here, however is that the
number of cores issignificantly moreefficient forthe variableAm background toroidal field strength is different for differ-
runs. Thus,wehadto redothe Am= 105 calculations, andorbit
ent values of Am. In Bai & Stone (2011), it was found
22 was chosen because the density-weighted stress at this time
was roughly equal to orbit 50 of the lower core version of this that with zero net verticalflux, the stress levelincreases
calculation. with increasing net toroidal flux (which is mostly con-
7
profileofthetotalstress(whichhasagain,beenaveraged
in time and for all x and y) in Fig. 5.
The stress profile reveals the same behavior as that
in Fig. 3; there is a general trend of increasing stress
with increasing Am. Furthermore, this increase occurs
uniformly across all z. However, C300 and C1000 have
roughlythesamestressprofiles,andC10peaksataround
thesamevalueasC100. Examiningtheβ forthesepar-
y
ticular simulations, we see that C300 has a lower value
(stronger field) than does C1000. Similarly, C10 has a
significantly smaller β than C100. These results con-
y
firmthatitisindeedthelargerbackgroundtoroidalfield
strengththatmakethe stresslevelsinrunC10approach
that in run C100, and the stress in run C300 approach
that in run C1000. We note, however, that C1000 and
C10000 have both the same β profiles and the same
y
stress profiles. This could indicate that for Am > 1000,
theturbulencelevelsareapproachingthatofidealMHD.
The slightly higher stress for C1e5 would then be ex-
Figure 2. Density-weightedvolumeaverageofthetotal(Maxwell
andReynolds)stress,normalizedbythesquareofthesoundspeed, plained by its lower βy.
versus time for the standard 4H ×8H ×8H simulations. The It remains unclear why the background field strength
magentalinecorrespondstoAm=1,bluetoAm=10,redtoAm varies in the way that it does. Could this also be con-
= 100, orange to Am = 300, black to Am = 1000, and green to
Am = 104. The horizontal dashed line corresponds to the time- trolled by the value of Am? This is not unreasonable
averaged (from orbit 25 to 53) stress value for Am = 105 from considering that ambipolar diffusion already affects the
whichthe other runswereinitialized. After aninitialtransient of period of the toroidal field flipping. The question of ex-
∼50orbits,thesimulationswithAm&10aresustained. Thereis actly how Am and the dynamo are related is very open.
ageneral trendofincreasingstresslevelwithincreasingAm. The
Unfortunately, exploring it in detail would take us too
Am=1casehasturbulencethatdecays awayrapidly.
far from our goals in this paper, and so, we leave it for
future work.
The final diagnostic we employ is the two-point au-
tocorrelation function first used in the context of MRI
simulations in Guan et al. (2009). We employ this diag-
nosticforsimilarreasonsasthoseinSimon et al.(2012);
wewishtodeterminethestructureoftheturbulentmag-
netic field and check that the domain sizes we use are
sufficiently large to properly capture important turbu-
lent scales. Thus, we define the autocorrelationfunction
of the ith component of the perturbed magnetic field as
δB (t,x)δB (t,x+∆x)d3x
ACF(δB (∆x))= i i , (20)
i R δB (t,x)2d3x
i
R
where δB is the value of B after subtracting off the
i i
horizontal mean field. In equation form,
δB (x,y,z) B (x,y,z) B (z), (21)
i i i xy
≡ −h i
Figure 3. Time-averaged total stress (i.e., α) as a function of andtheaveragedenotedby isthehorizontalaverage.
xy
Amforthestandard4H×8H×8Hsimulations. Thereisageneral Note that we have defined thhie ACF to be normalizedby
trendofincreasingstresslevelwithincreasingAm.
its maximum value (at ∆x = ∆y = ∆z = 0). The ACF
of the total turbulent magnetic field is then defined as
served in unstratified simulations). Therefore, in our
ACF(δB)= ACF(δB )+ ACF(δB )+ ACF(δB ). The
stratified simulations, two effects are expected to deter- x y z
overbar denotes a time average done from orbit 100 to
minethesaturatedstressvalues: thevalueofAmandthe
125 in all cases.
background toroidal field strength, set by the dynamo.
Fromthe figure,it appearsthatthe Am =100ACFis
To demonstrate this effect more robustly, we calculate a
roughly consistent between a domain size of 2H 4H
versionoftheplasmaβ forthebackgroundtoroidalfield, × ×
8H and4H 8H 8H,thoughthereisaslightdifference
× ×
in the size of the tilted centroid. However,as we will see
β 2 P / B 2, (19)
y ≡ h i h yi shortly, 2H 4H 8H is actually too small for Am =
× ×
where the overbars indicate a time average (from orbit 100. The standard box size, 4H 8H 8H, appears
× ×
100onwards)and the angledbracketsdenote anaverage to properly contain the ACF for Am = 100, but not as
over x and y. This quantity is representative of the am- well for Am = 10. The centroid of the Am = 10 case is
plitude of the oscillatingbackgroundtoroidalfield. β is larger and appears to have a longer tail that intersects
y
afunctionofzonly,andweplotitalongwiththevertical the toroidal boundary. Going to an even larger domain,
8
Figure 4. Space-time plot of By averaged over x andy forAm= 10 (upper left), Am= 100 (upper right), Am =300 (lower left), and
Am=1000(lowerright). The“butterfly”dynamoispresentinallsimulations,buttheperiodoftheBy flippingincreaseswithdecreasing
Am. Inparticular,theperiodis∼40−50orbitsforAm=10and∼15−20orbitsforAm=100. Fortheother twocases, theperiodis
∼10orbits,equal tothatinidealMHDcalculations.
8H 16H 8H,theACF(B)forAm=10looksmorewell Figure 8 shows the W evolution for these two domain
Rφ
× ×
contained, though the very end of the tail does appear sizes with Am = 10. The use of a smaller domain size
totouchthe toroidalboundary. Thiseffectisnotasdra- does not appear to make a difference for this value of
matic as in the smaller domain. Going to an even larger Am. Furthermore, the B space-time plot for the larger
y
domain and running Am = 10 is prohibitively expensive domain looks very similar to the smaller domain. Evi-
given our current computational resources. dently, we can get away with using a smaller domain for
Returning to the smallest domain, we note some odd Am = 10. However, these ACFs suggest that caution
features. Despite the reasonable ACF, an inspection of be used when running ambipolar diffusion shearing box
the stress history and the α value (see Table 1) show calculations.
this calculation to be quite different than Am = 100 at One final thing to note from ACF(B) is that as Am is
larger domain sizes. An examination of the space-time decreased, the tilted centroid component appears to be-
evolution of B , Fig. 7, brings the point home further, comemoreelongated(hencetheneedforlargerdomains)
y
as it indicates that the dynamo behavior is completely and less tilted with respect to the y axis.
shutoffforthisparticularrun. Thisisagaininconsistent Tosummarizethissection,wefindtheturbulentstress
with the larger domain Am = 100. Thus, we conclude level dependence on Am in vertically stratified simula-
that 2H 4H 8H is too small of a domain for Am = tions to be generally consistent with the results of un-
× ×
100 and will likely be too small for smaller values of Am stratified simulations (Bai & Stone 2011); α increases
as well. with Am, and for Am . 1, there is no turbulence. We
As indicated by these results, the standard box size of also find that as Am is decreased, larger domain sizes
4H 8H 8H is large enough for all of our simulations are needed in order to properly capture the turbulent
× ×
except for Am 10. It is computationally too expen- structures represented by the ACF. These runs serve as
≤
sive for us to run all of our simulations at the larger a baseline for interpreting the results from the next sec-
8H 16H 8H. So, we elect to use the smaller size tion, in which Am varies spatially and temporally based
× ×
as our standard, and we compare the evolution of the upon our model for surface layer ionization.
stress between the 4H 8H 8H and 8H 16H 8H
× × × ×
domains for Am = 10 to justify using a smaller domain
for comparison between Am = 10 and other Am values. 3.2. Variable Am Calculations
9
plots,thereisasignificantregioninzoverwhichtheMRI
isindeedactive. Theverticalstructureisconsistentwith
what we would expect from our ionization profile; there
is a significant “ambipolar dead zone” corresponding to
where Am = 1 and the higher Am regions correspond
to turbulent activity. We will calculate an actual mass
accretionrate atboth30AU and100AUin Section4.1.
Despite there being very little Maxwell stress in the
ambipolardeadzone,thereisstillReynoldsstressinthis
region. Thisstressislikelyproducedbytheactivelayers,
similartowhatisobservedinsimulationsthatincludean
Ohmic dead zone (Fleming & Stone 2003). It is possi-
ble,however,thatasinrunC1,the Reynoldsstresshere
is simply left over from the turbulent state from which
the run was restarted. To test this notion, we restarted
Z100AU at orbit 85 and set all perturbed velocity com-
ponentstozerothroughoutthedomain(theshearprofile
is of course maintained). We find that the active MRI
layers do indeed induce velocity fluctuations within the
ambipolardeadzone,whichleadsto apositiveReynolds
stress on average.
Finally, the horizontally averaged B space-time plot
y
shows interesting behavior. Within the ambipolar dead
zone, the toroidal field remains fixed in sign, though the
magnitude appears to be decreasing. Between 1 and 2
scale heights above and below the mid-plane, the usual
MRI dynamo reappears, with a B oscillation period of
y
10 orbits, identical to ideal MHD. Again, from the
∼
Maxwell stress plot, it is obvious that this region is tur-
bulent. For z >2H, however,there is a strong toroidal
| |
fieldthatremainsstationary. Thereareslightchangesin
Figure 5. Verticalprofilesofβy asdefinedbyequation(19)(top) magnitudeasthetoroidalfieldfromtheturbulentregion
and the total stress normalized by the square of the sound speed
propagates outwards.
(bottom). The quantities have been averaged over x and y and
over time from orbit 100 onwards as described in the text. The Examining the same diagrams for Z30AU reveals very
bluelinecorrespondstoAm=10,redtoAm=100,orangetoAm similarbehavior. ThereisathinlayerofstrongMaxwell
= 300, black to Am= 1000, green to Am= 104, and dashed line andReynoldsstresslocatedaround z 2H,alongwith
toAm=105(thetimeaverageforthisrunisdonefromorbit25to | |∼
apositiveReynoldsstressinducedintheambipolardead
theendofthecalculation). Thestressincreases withAmroughly
uniformlyatallheights. Thereisnoobvioustrendbetweenβyand zone by the active layers. In this calculation, however,
Am. the active region is much narrower in z when compared
to Z100AU. This is because the FUV photons do not
Wenowturnourattentiontothetwocalculationswith ionize the disk quite as deep as at R = 100 AU, and as
spatially and temporally varying Am (the “Z” simula- is usual, the MRI appears to be inactive for z & 2H.
| |
tions in Table 1). As stated previously, these simula- Thus, the MRI-active region is forced to be contained
tions adopt more realistic non-constant Am values and within a smaller z region.
directly address the questions such as “how vigorous is One thing to question is whether or not the phe-
outerdiskMRI-driventurbulence?”or“whatisthemass nomenology of Simon et al. (2011b) could play a role
accretion rate in the outer disk?” under the assumption here if we were to integrate these simulations further.
that the outer disk is not threaded by any net vertical In that paper, the authors found that the toroidal and
magnetic field. We run them all at the largest domain radial field in a resistivity dominated mid-plane region
size,8H 16H 8H,becausethereareregionsofAm 1 would evolve in time. Occasionally, the toroidal field
× × ≤
in these calculations. would grow to strong enough levels that the MRI would
Asbefore,webeginbyexaminingthedensity-weighted be temporarily reactivated, after which the Ohmic re-
stress normalized by the square of the sound speed, as sistivity would quench the turbulence again. In these
shown in Fig. 9. It is clear from the figure that the simulations, we do see that the toroidal field in the am-
stress levels are lowered dramatically compared to the bipolar dead zone changes in amplitude over time. Is it
ideal MHD case. Indeed, by calculating the average of possible, then, that the field could grow to large enough
the curve from orbit 72 onward, we find that the α val- values to eventually re-activate the MRI in that region?
ues are 10−3, one order of magnitude below what is Integratingoursimulationsoutformanyhundredsofor-
∼
expected from observations (Hartmann et al. 1998). bitsisprohibitivelyexpensive,andtherefore,wewillnot
While the turbulence levels initially decrease drasti- be able to observe any such variability in our stratified
cally, it has not completely died away. Consider the simulations. Instead, we have run an unstratified shear-
space-time diagrams in the (t,z) plane for various hor- ingboxofsize 8H 16H H withAm= 1anduniform
× ×
izontally averaged quantities of run Z100AU, as shown toroidal field of strength β = 10, chosen to determine if
in Fig. 10. From both the Maxwell and Reynolds stress a strong toroidal field MRI will be active with Am = 1.
10
Figure 6. Autocorrelationfunction(ACF)ofthemagneticfield,asdefinedbyequation(20),forsimulations(fromlefttoright)withAm
=100andsize2H×4H×8H,Am=100andsize4H×8H×8H,Am=10andsize4H×8H×8H,andAm=10andsize8H×16H×8H.
AsAmisdecreased, largerandlargerdomainsareneededtoproperlycontaintheACF.Furthermore,thetiltedcentroidfeaturebecomes
lesstiltedwithrespecttothey axisandmoreelongated asAmisdecreased.
Figure 7. Space-time plot of By averaged over x and y for the
Am=100runatdomainsize2H×4H×8H. Thereisaremnant
dynamobehaviorthatrapidlydiesawayasthesimulationadjusts
to the new value of Am. Eventually, the dynamo activity ceases
altogether, which is inconsistent with the larger domain counter-
part of Am = 100. This domain size appears to be too small to
properlycapturethedynamoatAm=100. Figure 8. Density-weightedvolumeaverageofthetotal(Maxwell
andReynolds)stress,normalizedbythesquareofthesoundspeed,
Westartthesimulationwitharelativelylargeamplitude versus time for the Am = 10 simulation at a domain size of
4H ×8H ×8H (solid line) and 8H ×16H ×8H (dotted line).
perturbation to the density and velocities, such that the
Note that the 8H ×16H ×8H run was restarted from orbit 22
initialconditionsshouldalreadybereasonablynonlinear. of a different “starter” simulation. For comparison purposes, we
As a control,we alsoranthis identical setupwith Am = translatedthesolutiontotherightby28orbits. Thecurvesshow
105. With Am = 1, we observeddecay of the initial per- anearlyidenticalevolution.
turbations and no development of any MRI turbulence,
whereas with Am = 105, the MRI becomes active and dead zone. Within this dead zone, the Reynolds stress
generates sustained turbulence for many orbits. Thus, dominatesoverthe Maxwellstress. Fromthe space-time
even in the presence of a strong field, Am = 1 is suffi- diagrams above, any nonzero Maxwell stress within this
cienttoquenchMRI-driventurbulence. We thereforedo region likely results from large scale correlations in the
notexpectthebehaviorobservedinSimon et al.(2011b) magnetic field rather than any sort of turbulence. Note
to occur in these simulations. that there are regions where the stress can go negative,
The vertical structure results can also be summarized and since the vertical axis is logarithmic, the curves are
in the time- and horizontally-averaged quantities as a simply cut off where the values drop below zero.
function of z, as is shown in Fig. 11. From the top two Thebottom twopanelsshowthe variousenergies(i.e.,
panels, it is obvious that there is a double-peak struc- thermal, kinetic, and magnetic) as a function of height.
ture to the stress; no doubt a result of the ambipolar The thermal pressure dominates over the vast majority