Global simulations of protoplanetary disks with ohmic resistivity and ambipolar diffusion

2.9 MB·English

by  Oliver Gressel

DRAFTVERSIONJANUARY23,2015
PreprinttypesetusingLATEXstyleemulateapjv.5/2/11
GLOBALSIMULATIONSOFPROTOPLANETARYDISKSWITHOHMICRESISTIVITYAND
AMBIPOLARDIFFUSION
OLIVERGRESSEL1,NEALJ.TURNER2,RICHARDP.NELSON3ANDCOLINP.MCNALLY1
1NielsBohrInternationalAcademy,TheNielsBohrInstitute,Blegdamsvej17,DK-2100,CopenhagenØ,Denmark
2JetPropulsionLaboratory,CaliforniaInstituteofTechnology,Pasadena,CA91109,USA
3AstronomyUnit,QueenMaryUniversityofLondon,MileEndRoad,LondonE14NS,UK
DraftversionJanuary23,2015
ABSTRACT
5 Protoplanetary disks are believed to accrete onto their central T Tauri star because of magnetic stresses.
1 RecentlypublishedshearingboxsimulationsindicatethatOhmicresistivity,ambipolardiffusionandtheHall
0 effectallplayimportantrolesindiskevolution. Inthepresenceofaverticalmagneticfield,thediskremains
2 laminarbetween1–5au,andamagnetocentrifugaldiskwindformsthatprovidesanimportantmechanismfor
removing angular momentum. Questions remain, however, about the establishment of a true physical wind
n
solution in the shearing box simulations because of the symmetries inherent in the local approximation. We
a
J presentglobalMHDsimulationsofprotoplanetarydisksthatincludeOhmicresistivityandambipolardiffusion,
wherethetime-dependentgas-phaseelectronandionfractionsarecomputedunderFUVandX-rayionization
2
withasimplifiedrecombinationchemistry. Ourresultsshowthatthediskremainslaminar,andthataphysical
2
windsolutionarisesnaturallyinglobaldiskmodels. Thewindissufficientlyefficienttoexplaintheobserved

accretionrates. Furthermore,theionizationfractionatintermediatediskheightsislargeenoughformagneto-
]
P rotational channel modes to grow and subsequently develop into belts of horizontal field. Depending on the
E ionizationfraction, thesecanremainquasi-global, orbreak-upintodiscreteislandsofcoherentfieldpolarity.
. ThediskmodelswepresenthereshowadramaticdeparturefromourearliermodelsincludingOhmicresistivity
h
only. ItwillbeimportanttoexaminehowtheHalleffectmodifiestheevolution,andtoexploretheinfluence
p
thishasontheobservationalappearanceofsuchsystems,andonplanetformationandmigration.
-
o Subjectheadings: accretion,accretiondisks–MHD–methods: numerical–protoplanetarydisks
r
t
s
a 1. INTRODUCTION as sodium and potassium provides strong coupling between
[ the gas and magnetic field (Umebayashi & Nakano 1988),
Understanding the complex dynamical evolution of proto-
non-ideal MHD effects such as Ohmic resistivity, ambipolar
1 planetary disks (PPDs) is of key interest both for building a
diffusion(AD)andHalldriftareexpectedtobeimportantdue
v comprehensivetheoryofplanetformation,andasameansof
to the low ionization levels of the gas (e.g. Blaes & Balbus
1 explainingtheobservationallyinferredpropertiesoftheseob-
1994;Wardle1999;Sano, Miyama, Umebayashi, &Nakano
3 jects(seeTurneretal.2014,forarecentreview). Forexam-
2000; Balbus & Terquem 2001). External sources of ioniza-
4 ple, PPDs are known to accrete gas onto their host stars at a
5 typicalrateof10−8±1M yr−1 (Gullbringetal.1998;Hart- tionsuchasstellarX-rays,FUVphotonsandgalacticcosmic
0 mannetal.1998)andhav(cid:12)elifetimesintherange∼3–10Myr rays are expected to ionize the disk surface layers, provid-
1. (Sicilia-Aguilar et al. 2006). More recent observations have inggoodcouplingthere–eventhoughrecentresultssuggest
thatcosmicraysmaybehighlyattenuatedbythestellarwind
0 indicatedthepresenceofturbulenceinthediskssurrounding
(Cleevesetal.2013). Nevertheless, deepwithinthediskthe
5 theyoungstarsTWHydraeandHD163296,basedonanaly-
evolution should be controlled by the non-ideal effects. The
1 sisoftheirmolecularlineemission(Hughesetal.2011). Fur-
recognition of this led Gammie (1996) to propose what has
: ther evidence for turbulent broadening comes, for instance,
v now become the traditional dead-zone model in which the
fromtheCOrovibrationalbandheadshapemeasuredbyNa-
i disksurfacelayersaccretebysustainingMRIturbulence,with
X jitaetal.(2009)inV1331Cygni.DuringtheTTauri(classII)
theshieldedinteriormaintaininganinertandmagneticallyde-
phase,self-gravityisexpectedtoprovideanegligiblecontri-
r coupled“deadzone”, wheretheMRIisquenchedbytheac-
a bution to angular momentum transport (due to the low disk-
tionofOhmicresistivityonly.Thepotentialimportanceofthe
to-star mass ratio), and disk accretion is instead believed to
othernon-idealeffectshaslongbeenrecognized(e.g.Sano&
bedrivenbymagneticeffects. Twopossiblesourcesofangu-
Stone2002;Salmeron&Wardle2003),butitisonlyrecently
lar momentum transport are through a magnetocentrifugally
that nonlinear shearing box simulations including ambipolar
driven wind (Blandford & Payne 1982; Pudritz & Norman
diffusion and the Hall term have been performed in relevant
1983), or through the operation of the magnetorotational in-
parameterregimes,leadingtoamodifiedpictureofhowdisks
stability(MRI,Balbus&Hawley1991),whosenonlinearout-
accrete that deviates significantly from the traditional dead
comeintheideal-MHDlimitismagnetohydrodynamicturbu-
zoneone.
lence (Hawley & Balbus 1991; Hawley, Gammie, & Balbus
Generally speaking, it is expected that in disk regions be-
1995).
tween 1–5au, Ohmic resistivity will be dominant near the
Exceptfor theinnermost regionsofPPDs, wherethetem-
peratureT ∼>1000K,andtheionizationofalkalimetalssuch midplane, the Hall effect will be dominant at intermediate
disk heights and ambipolar diffusion will be most important
inlowdensityregionshigherupinthedisk(Salmeron&War-
[email protected] (OG); [email protected] (NJT);
[email protected](RPN);[email protected](CPMN) dle2008). Thehighestaltitudesurfacelayersareexpectedto
2
beionizedbystellarFUVphotons,andassuchwillevolvein that are threaded by vertical magnetic fields and which in-
the ideal MHD limit (Perez-Becker & Chiang 2011). Shear- cludethiscombinationofnonidealMHDeffects,andassuch
ingboxsimulationspresentedbyBai&Stone(2013)thatin- are most comparable with the shearing box simulations pre-
clude Ohmic resistivity and ambipolar diffusion for models sentedbyBai&Stone(2013),hereafterBS13. Afundamen-
computed at 1au demonstrate that AD quenches MRI tur- talquestionraisedbytheshearingboxsimulationsiswhether
bulence. They conclude that disks will be completely lam- or not a physical wind solution with the correct field geom-
inar between 1–5au with angular momentum transport oc- etry emerges spontaneously in global disk simulations in a
curringbecauseamagnetocentrifugalwindislaunchedfrom waythatisnotgenerallyobservedinshearingboxesbecause
near the disk surface. For reasonable assumptions about the oftheirradialsymmetry. Probablythemostimportantresult
magneticfieldstrengthandgeometry,accretionratesinagree- containedinthispaperisthatphysicalwindsolutionsdoin-
ment with the observations are obtained. The inclusion of deed arise spontaneously in our global simulations, demon-
the Hall term in the presence of a vertical magnetic field in- stratingtherobustnessofmanyofthefeaturesobtainedinthe
troducesdichotomousbehavior,arisingfromthefactthatthe shearingboxsimulationsofBS13.
coupleddynamicsdependsonthesignofΩ·B,whereΩand This paper is organized as follows: In Section 2 we de-
Bsignifythediskangularmomentumvectorandtheambient scribe the numerical method, the disk model, as well as the
magnetic field direction, respectively. Shearing box simula- ionization chemistry. Our results are presented in Section 3,
tions presented by Lesur et al. (2014) and Bai (2014) show where we mainly focus on the emerging wind solutions and
that when Ω · B > 0, the Hall effect leads to the amplifi- thedynamicalstabilityandevolutionofformingcurrentlay-
cationofthehorizontalmagneticfieldswithinapproximately ers. In Section 3.7, we will moreover report an instability
twoscaleheightsofthediskmidplane,andthegenerationof thatisdrivenbythesharptransitionintotheFUVdominated
alargescaleMaxwellstressthroughmagneticbraking. Inad- layer, and in Section 3.8 we assess whether secondary insta-
dition, the magnetocentrifugally driven wind is also seen to bilitiescandrivenon-axisymmetricevolution. Wesummarize
strengthen. When Ω · B < 0, however, magnetic stresses our results in Section 4 and discuss implications for planet
and winds are seen to weaken relative to the opposite case, formationtheoryinSection5.
and relative to the evolution observed when Hall effects are
neglected.QuitehowthisHalldichotomymapsontoobserva-
tionsofPPDsremainsanopenquestion.
The internal dynamics of protoplanetary disks in the re-
gion 1–5au are clearly of great significance for planet for-
mation. The growth and settling of grains depends on the
level of turbulence, with a laminar disk or one displaying 2. METHODS
weak turbulence providing the most favorable conditions – We present magnetohydrodynamic (MHD) simulations of
althoughitshouldbenotedthatsometurbulentmixingisre- protoplanetary accretion disks employing 2D axisymmetric
quired to maintain a population of small grains in the sur- and3Duniformly-spacedspherical-polarmeshes. Inthefol-
faceregionsofPPDssothattheirspectralenergydistribution lowing, the coordinates (r,θ,φ) denote spherical radius, co-
can be reproduced by radiative transfer models (Dullemond latitude and azimuth, respectively. We moreover ignore ex-
& Dominik 2005). The collisional growth of planetesimals plicit factors of the magnetic permeability µ in our nota-
0
hasbeenshowntobeaffectedstronglybythelevelofturbu- tion. Simulations were run using the single-fluid NIRVANA-
lence (Gressel et al. 2011, 2012), and the migration of low IIIcode,whichimplementsasecond-order-accurateGodunov
mass planets is also highly sensitive to the presence or oth- scheme (Ziegler 2004). The code is formulated on orthogo-
erwise of turbulence, with significant stress levels being re- nal curvilinear meshes (Ziegler 2011) and employs the con-
quired to unsaturate corotation torques (Paardekooper et al. strainedtransportmethod(Evans&Hawley1988)tomaintain
2011;Baruteauetal.2011). thedivergence-freeconstraintofthemagneticfield. Asanal-
Thispaperisthelatestinaseriesinwhichweareexploring terationtothepubliclyavailabledistributionofthecode, we
the influence of magnetic fields and non ideal MHD effects hereadopttheupwindreconstructiontechniqueofGardiner&
ontheformationofplanets,withtheeventualgoalofproduc- Stone(2008)toobtaintheedge-centeredelectricfieldcompo-
ingcomprehensivemodelsofPPDsthatwillbeusedtostudy nents needed for the magnetic field update (also see Skinner
planet building and evolution. In earlier work (Nelson & &Ostriker2010).Thismodificationbecamenecessarytotake
Gressel2010;Gresseletal.2011,2012)wehaveusedacom- advantageofthemoreaccurateHLLDapproximateRiemann
bination of global and shearing box simulations to examine solver of Miyoshi & Kusano (2005), which offers advanced
thedynamicalevolutionofdustgrains,bouldersandplanetes- numericalaccuracyforflowsinthesubsonicregime.
imalsinturbulentdisks, withandwithoutdeadzones. More Our numerical model is based on solving the standard re-
recently, wehavestudiedtheinfluenceofmagneticfieldson sistiveMHDequationsincludingOhmicresistivityaswellas
gap formation and gas accretion onto a forming giant planet anelectromotiveforceresultingfromthemutualcollisionof
usingglobalsimulationsthatalsoincludedOhmicresistivity ionsandneutrals. Giventhetypicalnumberdensitieswithin
(Gressel et al. 2013). In this paper, we include both Ohmic PPDs,theapplicabilityofthestrong-couplingapproximation
resistivityandambipolardiffusioninglobaldisksimulations, iswarrantedbydetailedchemicalmodeling(Bai2011). Ac-
andfollowthedynamicalevolutionoftheresultingdiskmod- cordinglythegasdynamicscanbemodeledbyasingle-fluid
els as a precursor to examining how ambipolar diffusion af- representing the motion of the neutral species, and the ad-
fects gas accretion onto a giant planet. Of particular interest ditional term can simply be evolved in a time-explicit fash-
isthenatureoftheaccretionflowinthedisk,thenatureofany ion (Choi et al. 2009). The total-energy formulation of the
magnetocentrifugalwindthatislaunched,andhowthesevary NIRVANA-III codeisexpressedinconservedvariablesρ,ρv,
as a function of small changes in model parameters. These ande≡(cid:15)+ρv2/2+B2/2,andifwedefinethetotalpressure,
arethefirstquasi-globalsimulationstobeconductedofPPDs p(cid:63),asthesumofthegasandmagneticpressures,wecanwrite
3
ourequationsystemas look-up table, that has been derived self-consistently from a
detailed chemical model accounting for grain-charging and
∂ ρ+∇·(ρv)=0,
t gas-phasechemistry. Theionizationratesthatenterthisequi-
∂ (ρv)+∇·[ρvv+p(cid:63)I−BB]=−ρ∇Φ, libriumchemistryaregovernedbyionizingradiationentering
t
thedisk. Thecolumndensitiesthatgoverntheattenuationof
∂ e+∇·[(e+p(cid:63))v−(v·B)B]=∇·[(E +E )×B] the external ionizing radiation are also computed on-the-fly,
t O AD
−ρ(∇Φ)·v+Γ, thatis,theyreflecttheinstantaneousstateofthedisc. Inprin-
ciple,thisallowsforaself-limitingoftheemergingwindvia
∂tB−∇×(v×B+EO+EAD)=0, (1) shieldingofionizingradiation(Bans&Ko¨nigl2012).
Inadditiontosomeminormodificationstothegraincharg-
wheretheelectromotiveforcesstemmingfromtheOhmicand
ingprescription(asdetailedbelow),wehaveimprovedthere-
ambipolardiffusiontermsaregivenby
alismofourionizationmodelcomparedtoourpreviouswork.
E ≡−η (∇×B), and (2) The main alterations concern the inclusion of two additional
O O
radiation sources, an un-scattered direct X-ray component,
(cid:104) (cid:105)
E ≡η (∇×B)×bˆ ×bˆ, (3) andhardUVphotons. Thesehaveincommonashallowpen-
AD AD
etrationdepthbutcomparativelyhighfluxlevels,thusmainly
(withbˆ ≡ B/|B|theunitvectoralongB),respectively. The affectingtheionizationlevelofthedisk’ssurfacelayers.
gravitational potential Φ(r) ≡ −GM /r is a simple point-
(cid:12) 2.1.1. FUVionizationlayer
mass potential of a solar-mass star. We moreover obtain the
gaspressureasp = (γ −1)(cid:15),whereγ = 7/5isappropriate Adopting a very simple prescription based on the recent
foradiatomicmoleculargas. ThesourcetermΓ,mimicking study by Perez-Becker & Chiang (2011), we have added the
radiativeheatingandcooling,isincludedtorelaxthespecific effect of an FUV ionization layer based on an assumed ion-
thermalenergydensity(cid:15)toitsinitialradialtemperaturepro- izationfractionoff =2×10−5(representingfullionization
file. Therelaxationisdoneonafractionofthelocaldynam- ofthegas-phasecarbonandsulfuratomssusceptibletolosing
ical time scale, which is short enough to avoid strong devia- electronswhen struckbyFUV photons), andacollision rate
tionsfromtheequilibriummodelbutlongenoughtosuppress coefficientof2×10−9m3s−1. WefurtherassumethatFUV
verticalcorrugationofthediskduetotheGoldreich-Schubert- photons penetrate to a gas column density of 0.03g cm−2.
Frickeinstability(seeNelson,Gressel,&Umurhan2013,for Due to their lower amplitude, we ignore any scattered, dif-
adetailedstudyofthisphenomenon). Weremarkthat,while fuseorambientFUVilluminationofthediskandonlyevalu-
thisinstabilityisphysicalinnature,itsappearancemaybeex- atesight-linespointingdirectlytowardsthecentralstar. Note
aggeratedinastrictlyisothermalsimulation,andmorerealis- that this deviates from the local simulations of BS13, where
ticmodelingincludingradiativetransfershouldbeemployed the vertical gas column was used for attenuation of the inci-
tostudyitsrelevance. dentflux. Ourtreatmentismotivatedbytheassumptionthat
Because we use a time-explicit method, large peak val- a large fraction of the FUV photons originate directly from
uesinthedissipationcoefficientsηO andηAD imposesevere thecentralstar(seeupperrightpaneloffigure9inBethell&
constraints on the numerically permissible integration time Bergin 2011, who study the forward scattering of FUV pho-
step. Weaddressthisbyusingastate-of-the-artsecondorder tonsindetail).
super-time-steppingschemerecentlyproposedbyMeyeretal.
(2012). Incomparisonwithconventionalfirst-ordermethods, 2.1.2. IlluminationbyX-rays
this Runge-Kutta-Legendre scheme is free of adjustable pa-
A similar treatment is adopted for the unscattered X-ray
rametersandhasbeenprovensuperiorintermsofbothaccu-
component, for which we use the fit to the Monte-Carlo
racy and robustness. To maintain the second-order accuracy
radiative-transferresultsofIgea&Glassgold(1999)givenby
ofourtimeintegration, weuseStrang-splittingforthediffu-
eq.(21)inBai&Goodman(2009). Deviatingfromourpre-
siveterms.
viouslocalsimulations(cf. eq(1)inGresseletal.2011),our
Inourpreviouslocalsimulations(cf. appendixB1inGres-
newX-rayilluminationis
sel et al. 2012), we have found that applying a constant cap
onηO canqualitativelychangethewayinwhichthetopand (cid:20) (cid:18)Σ (cid:19)α (cid:18)Σ (cid:19)α (cid:21)
bottom disk layers are coupled to each other. For a typical ζ =10−15s−1 exp− A +exp− B r−2
XR Σ Σ
windconfiguration,thehorizontalfieldcomponentsgenerally sc sc
change sign at the midplane. We imagine that in this situa- (cid:18)Σ (cid:19)β
tion a clipped ηO would affect the amount of field diffusing +6×10−12s−1 exp − C r−2 (4)
into the midplane region and conversely the strength of the Σab
current sheets above and below the magnetically decoupled
where Σ (Σ ) is the gas column to the top (bottom) of
layer. We therefore choose not to apply a cap on η . How- A B
O the disk surface, and Σ is the column density along radial
ever,welimitη suchthatΛ β ≡2Rm ≥0.1(also C
AD AD p AD rays towards the star. The radial column contains a contri-
seeBS13),whereβ ≡ 2p/B2 istheplasmaparameter,and
p bution from the inner disk (which is not part of the com-
Rm ≡c H/η istheequivalentofamagneticReynolds
AD s AD putational domain) and reaches in to an inner truncation ra-
number for AD. This we find greatly reduces the numerical
dius of r = 5R . Adopting values from Bai & Goodman
costwithoutnoticeablyalteringtheobtainedsolution. (2009), thtreshape(cid:12)coefficientsareΣ = 7×1023cm−2, and
sc
α = 0.65forscatteredX-rays, andΣ = 1.5×1021cm−2,
2.1. Externalionizationanddiskchemistry andβ = 0.4,forabsorptionofthedireacbtcomponent,respec-
The Ohmic and ambipolar diffusivities are evaluated dy- tively. Bothcontributionsadditionallydecaywiththesquared
namicallyforeachgrid-cellandarebasedonaprecomputed radius to account for the decrease in X-ray luminosity. To
4
floor acts to mimic a disk with realistic radiation thermody-
namics,wheretheirradiationheatingoftheupperdisklayers
naturallycreatesanincreasedpressurescaleheightandhence
a much shallower density profile. Because η (and hence
AD
Λ )dependsonthemagnetizationoftheplasma,theactual
AD
densityprofileinthediskhalohasaninfluenceonhowwell
thegasiscoupledtothemagneticfieldlines,whichinturnaf-
fectstheefficiencyofthemagneto-centrifugalwindlaunching
mechanism. This caveat should be kept in mind when inter-
pretingthemassloadingofthewind,whichweexpecttobea
functionofthethermodynamicsinrealprotoplanetarydisks.
2.2. Globaldiskmodel
We now describe the underlying protostellar disk model,
which is largely identical to the one used in Gressel et al.
(2013). The equilibrium disk model is based on a locally-
isothermal temperature, T, decreasing with cylindrical ra-
dius as T(R) = T (R/R )q. For q = −1, such a profile
0 0
leads to a constant opening angle throughout the disk as it
iscommonlyusedforglobalnumericalsimulationsofPPDs,
Figure1. ProfilesoftheinitialElsassernumbers(atradiusr = 1au)for and which provides us with the reference disk model here.
ohmicresistivity(lightorange), andambipolardiffusion(red). Solidlines
To study what effect disk flaring has on the absorption of
representthefiducialcaseofadust-to-gasmassratioof10−3,whereasdash-
dottedlinesshowprofilesforavalueof10−4.Wealsoplottheinitialprofile direct FUV photons, we also produce a moderately flaring
oftheplasmaβparameter(solidblackline). disk with q = −3/4. For this value, the resulting power-
law index for the local opening angle, h(R) ≡ H(R)/R, is
accountfortypicalmedianvaluesof1030ergs−1 inobserva- ψ ≡(q+1)/2=0.125,whichisinagreementwithobserva-
tionalconstraintsforthisvaluebasedonsub-mmobservations
tional data of luminosities in young star clusters such as in
(Andrews et al. 2009), and from multi-wavelength imaging
the Orion Nebula (Garmire et al. 2000), we enhance the in-
cident X-ray flux by a factor of 5× compared to the above (Pinte et al. 2008).1 The radial power-law exponents for the
gas surface density, Σ(R), are −1/2 for our fiducial model,
statedvalues. Inviewoffuturework,weenvisageadditional
and−3/8fortheflaringdiskmodel,respectively.
improvementstotheX-rayionizationmodelbyadoptingthe
Ourtemperaturedistributioniscomplementedbyapower-
new results of Ercolano & Glassgold (2013), who consider
law for the midplane density, ρ (R) = ρ (R/R )p, with
the enhanced ionizing effects of X-rays due to the presence mid 0 0
p=−3/2. Theequilibriumdisksolutionisgivenby
ofheavierelements(assumingsolarabundance). Thesealter-
ationshavebeenfoundtoleadtoenhancedionizingfluxesat (cid:18) R (cid:19)p (cid:18)GM (cid:20)1 1(cid:21)(cid:19)
intermediatecolumndensitiescomparedtotheoriginalresults ρ(r)=ρ exp (cid:12) − , and (5)
0 R c2 r R
ofIgea&Glassgold(1999). 0 s
2.1.3. Modificationstothediskchemistry (cid:34) (cid:18)H(cid:19)2 qR(cid:35)12
Ω(r)=Ω (R) (p+q) +(1+q)− , (6)
OurchemicalnetworkisbasedontheoneusedbyIlgner& K R r
Nelson(2006),labeledmodel4,withgrainsurfacereactions
and a simplified gas-phase reaction set involving one repre- which can be derived from the requirement of hydrostatic/-
sentative metal and one molecular ion. Small changes in re- dynamic force balance in the vertical and radial direc-
cent years include correcting the electron sticking probabili- tions, and where the Keplerian angular velocity ΩK(R) =
(cid:112)
ties for thegrain charge (Bai 2011)and consistently treating GM R−3/2. In the above equations, the square of the
(cid:12)
themolecularionsuchasHCO+ (withamassof29protons) isothermal sound speed results from the prescribed tempera-
sincethisisthelong-livedspeciesinaseriesofseveralreac- tureprofileasc2 = c2 (R/R )q,withaparameterc corre-
s s0 0 s0
tions. Furtherdetailsonthereactionnetworkanddiffusivities spondingtoavalueofH(R)≡c Ω−1 =0.05R.
s K
canbefoundinsection2.2. ofLandryetal.(2013)aswellas Guidedbypreviousresultsoflocal3Dsimulationsthatex-
section4.2ofMohantyetal.(2013). hibit quasi-axisymmetric structure, we mainly focus on 2D-
TheresultinginitialionizationprofileexpressedinEls√asser axisymmetricsimulations. Ourcomputationaldomaincovers
numbersΛO/AD ≡vA2z(ΩηO/AD)−1 withvAz ≡Bz/ ρis r ∈ [0.5,5.5]au. Inthelatitudinaldirection, thegridcovers
plottedinFig.1ataradiusof1au,togetherwiththeheight- eightpressurescaleheightsoneachsideofthediskmidplane,
dependentplasmaβparameter.Theimplicationsofthisfigure that is, θ ∈ [π/2 − ϑ,π/2 + ϑ], with ϑ = 8H/r. In the
foroursimulationsarediscussedinSect3.1. Theregionwith case of the flaring disk, the coverage of scale heights varies
|z|∼>5H,wheretheplasmaparameterbecomesconstantwith from 8H at r = 1au to 6.5H at r = 5au. The grid reso-
heightiscausedbyapplyingalowerlimitonthegasdensity lution is N ×N ×N = 512×192×64, which means
r θ φ
toavoidexcessiveAlfve´nspeedsinthediskhalo. Following thatweareabletoresolverelevantMRIwavelengths. Inthe
BS13, we chosethislimit at 10−6 (inour globalmodel, this axisymmetricmodelsweuseN ×N = 1024×384cells,
r θ
valueiswithrespecttotheinitialmidplanevalueatanygiven corresponding to ≥ 24 grid points per pressure scale height,
cylindrical radius). We remark that the equilibrium density
profile is artificially steep in our isothermal disk model with 1Atleastwhenassumingwell-mixeddustgrains,sincebothobservations
temperatureconstantoncylinders. Inthisrespect,thedensity aresensingthedustcontentofthediskratherthanthegasdensity.
5
matching the resolution of the three-dimensional box simu- component to zero and allowing non-zero parallel field), we
lations of BS13. For the two non-axisymmetric simulations, hereusepseudovacuumconditions,converselyenforcingthe
theazimuthalextentisrestrictedtoφ∈[0,ϕ],withϕ=π/2 perturbedparallelfieldtovanishatthesurfaceandonlyallow-
(aquarterwedge),orπ/4tolimitthecomputationalexpense. ingaperpendicularperturbedfield. Notethatweexcludethe
We furthermore note that in the ideal-MHD case, where tur- initialnet-verticalfieldfromtheproceduresuchthatonlyde-
bulencedevelops, theazimuthaldomainsizehasbeenfound viationsfromthebackgroundfieldaresubjecttothenormal-
to have an effect on the final saturated state(Beckwith et al. field condition. The vertical flux threading the disk is thus
2011;Flocketal.2012;Sorathiaetal.2012),atleastforthe preserved. Since the disk’s upper layers are magnetically
caseofzeronetflux. dominated, and theLorentz forceacts toalign theflow lines
Becauseourionizationmodeldependsonrealphysicalval- withthemagneticfield,lettingthefieldlinescrossthedomain
uesofthegascolumndensity,wehavetochoseameaningful boundaryisofcourseaprerequisiteforlaunchinganoutflow.
normalizationfactor ρ , which wefix accordingto asurface While we realize that forcing the radial and azimuthal com-
0
density of Σ = 150gcm−2 at the location R = 5au, com- ponentsoftheperturbedfieldtovanishattheboundarymay
patiblewiththetypicalminimum-masssolarnebula(MMNS unduly restrict the magnetic topology of the emerging wind
Hayashi 1981). In physical units, the disk temperature is solution, we postpone the study of less-restrictive but more
T = 540K at a radius of R = 1au, and T = 108K at cumbersomeboundariestofuturework.
5au, resembling typical expected conditions within the pro- Unlikeinaradially-periodiclocalshearing-boxsimulation,
tosolardisk.Whileourvaluesforρ,T andΣaresimilartothe our global model is critically affected by the inner radial
MMSNvaluesat5au,ouradoptionofdifferentvaluesforthe boundary condition. In a real protoplanetary disk, we can
radialpower-lawindicesmeansthatthesevaluesaredifferent expectthatalkalimetalswillbethermallyionizedinthisre-
at1aucomparedtothoseadoptedbyBS13. gion, leading to the development of the MRI on timescales
Our model disk is initially threaded by a weak vertical shortcomparedwiththeorbitaltimescaleattheinnerdomain
magnetic field B (R) corresponding to a midplane β ≡ boundaryofourmodel. Thisposesthequestionhowtobest
z p0
2p/B2 = 105 independent of radius for our fiducial disk. interface the MRI-turbulent inner disk with the magnetically
To achieve this, the vertical magnetic flux has to vary as a decoupled midplane of the outer disk that we model. It is
power-law in radius, taking into account the radial distribu- likely that MRI-generated fields can efficiently leak into the
tion of the midplane gas pressure that itself depends on the outerdiskviamagneticdiffusion. Inafirstattempttoaccount
densityandtemperature. Inourdiskmodel,thegaspressure, for the MRI-active inner disk, we gradually taper the diffu-
p(R) decreases outward, implying that B (R) falls off with sivity coefficients to zero within the inner 0.5au of our disk
z
radius,too. Whilesuchaconfigurationisgenerallyexpected model.Studyingtheinneredgeofthedeadzonewillhowever
whenaccountingforinwardadvectionandoutwarddiffusion require dedicated simulations (similar to the ones performed
oftheembeddedverticalflux(Okuzumietal.2014),ourpar- byDzyurkevichetal.2010)includingthistransitionregion.
ticular choice of keeping the relative field strength constant
3. RESULTS
is of course arbitrary. To preserve the solenoidal constraint
to machine accuracy, we compute the poloidal field compo- Themainaimofourpaperistoestablishthelaminarwind
nents from a suitably defined vector potential A (r,θ). The solutions that BS13 previously found in local shearing box
φ
initialdiskmodelisperturbedwithrandomwhite-noisefluc- simulations in the context of global disk simulations. In the
tuationsinthethreevelocitycomponents. Thermsamplitude interest of economic use of computational resources and to
oftheperturbationischosenatonepercentofthelocalsound guarantee adequate numerical resolution of our global mod-
speed. Wefurthermoreaddawhite-noiseperturbationinthe els,weprimarilyfocuson2Daxisymmetricsimulations,but
magnetic field with an rms amplitude of a few µG, which is havealsorunthree-dimensionalsimulationstocheckfornon-
on the sub-percent level compared with the net-vertical field axisymmetric solutions. Since all our runs include a net-
oftypicallyafewten mG. vertical flux, axisymmetry is warranted to obtain a reading
on the development of the MRI. In cases where the solution
2.3. Boundaryconditions
proves laminar, axisymmetry is likely to produce a reason-
To complete the description of our numerical setup, we ably accurate picture. We address the possibility of non-
need to specify boundary conditions (BCs). For the fluid axisymmetric secondary instabilities using a limited set of
variables,weemploystandard“outflow”conditionsatthein- three-dimensionalcalculations,describedinSection3.8.
ner and outer radial domain edges, r and r . This type The simulations conducted for this work are listed in Ta-
in out
is a combination of a zero-gradient condition in the case of ble 1, where we summarize key model parameters. We as-
nˆ·v(r ,θ)>0 (wherenˆ denotestheoutward-pointingnor- sume a fiducial dust-to-gas mass ratio of 10−3, that is, a de-
in
malvector),andreflectiveBCsintheoppositecase,thuspre- pletionbyafactoroftencomparedtointerstellarabundances.
venting inflow of material from outside the domain. At the StartingfromtheclassiclayeredPPD(model‘O-b6’)includ-
upperandlowerboundaries(thatis,intheθdirection),wefur- ingonlyOhmicdiffusivity,wecomparethisstandard“dead-
thermorebalancetheghostzonevaluessuchthatahydrostatic zone”diskwithacorrespondingmodel,‘OA-b6’,additionally
equilibriumisobtained. Thishelpstocontrolartificialjumps includingtheeffectsofambipolardiffusion. Ourfiducialref-
ofthefluidvariablesadjacenttothedomainboundaryasthey erencemodelis‘OA-b5’,withamidplaneplasmaparameter
arefrequentlyencounteredwithfinite-volumeschemes. of105. Inthepresenceofcombinedambipolardiffusionand
In contrast to our earlier global simulations of layered ac- Ohmicresistivity, weobservethewaningoftheMRI,which
cretion disks subject to Ohmic resistivity and containing an isinsteadreplacedbyalaminarwindsolution. Unlikeinge-
embedded gas-giant planet (Gressel et al. 2013), we here ometricallyconstrainedshearingboxsimulations(BS13,Bai
makeadifferentchoicefortheverticalmagnetic-fieldbound- 2013),wenaturallyobtainafieldconfigurationwithfieldlines
arycondition. WhereaswepreviouslyappliedmagneticBCs bendingoutwardonboth“hemispheres”ofthedisk. Notably,
oftheperfectconductortype(thatis,forcingthenormalfield thistopologyproducesthincurrentlayers,whichhaveprevi-
6
causestheambipolardiffusiontoincrease.Theirdiskthenre-
Table1
laxestocompletelylaminarstateinwhichangularmomentum
Summaryofsimulationruns.
transportisdrivenbyamagneto-centrifugalwind. Thelarger
Runlabel Ohm AD βp0 d/g q Resol. valueofΛAD inourmodelsmayallowMRIturbulencetobe
O-b6 ◦ − 106 10−3 -1 1024×384 sustained in these regions, or may instead lead to quenching
OA-b5 ◦ ◦ 105 10−3 -1 1024×384 oftheMRIafterithasreachedamorenonlinearstageofde-
OA-b6 ◦ ◦ 106 10−3 -1 1024×384 velopment. In this paper, we define our fiducial model to be
OA-b7 ◦ ◦ 107 10−3 -1 1024×384 one in which the dust-to-gas mass fraction is 10−3, and the
OA-b5-d4 ◦ ◦ 105 10−4 -1 1024×384 Elsassernumbersforthiscaseareshownbythesolidlinesin
OA-b5-flr ◦ ◦ 105 10−3 -3/4 1024×384 Fig. 1. The larger dust concentration reduces the gas-phase
OA-b5-flr-nx ◦ ◦ 105 10−3 -3/4 512×192×128 electron fraction and Elsasser numbers, and the local max-
OA-b5-nx ◦ ◦ 105 10−3 -1 512×192×64 imum in Λ now peaks moderately above unity (but still
AD
attainingalargerpeakvaluethanthefiducialmodelinBS13).
Modelsarelabeledaccordingtotheincludedmicro-physicaleffects(firsttwo
letters)andthestrengthofthenet-verticalmagneticfield(expressedinterms
ofthemidplanevalueβp0,prefixedwiththeletter‘b’).Furtherlabelsreferto 3.2. Traditionaldeadzonemodel
thedust-to-gasmassratio(‘d/g’,prefixedwiththeletter‘d’),thepower-law
Webegindiscussionofthesimulationresultsbyconsider-
index,q,oftheradialtemperatureprofile(‘flr’for“flaring”),andwhetherthe
azimuthaldimensionisincluded(‘nx’for“non-axisymmetric”). ing the 2D axisymmetric run O-b6, for which Ohmic resis-
tivity is included and ambipolar diffusion is neglected. To
ouslybeendiscussedbyBS13. Thestabilityandevolutionof introduceoursetup,andgivethereaderanimpressionofthe
thesecurrentlayerswillbeonefocusofourpaper. general appearance of our PPD model, Fig. 2 visualizes the
We perform additional analysis on the influence of further magneticfieldandpoloidalflowvectorsfromthemodel. The
key input parameters. With the global disk model allowing color coding of the azimuthal field strength shows the MRI-
direct illumination from the central star, it becomes possible turbulentsurfacelayerswithtangledpoloidalfieldlinessuper-
toaddressthequestionofhowthedisk’sionizationisaffected imposed in white. In the outer disk, where the MRI has not
by disk flaring. This is studied in model ‘OA-b5-flr’, where fully developed yet, the parity of the horizontal field is odd
we use a power-law index of q = −3/4 instead of q = −1 withrespecttoreflectionatthez =0line,consistentwiththe
toobtainamoderatelyflaringdisksurfaceasissupportedby winding-upofverticalfield,B ,bytheweakdifferentialrota-
z
observations (Pinte et al. 2008; Andrews et al. 2009). The tion∂ Ωofourdiskmodel.Whiletheinnerdomainboundary
z
role of dust depletion, driven by processes such as coagula- also shows an odd field symmetry, intermediate sections of
tion into larger grains and/or sedimentation, is considered in thediskshowamixedparity,reflectingthepresenceofMRI
model ‘OA-b5-d4’ with a reduced dust-to-gas mass fraction modeswithbothevenandoddverticalwavenumbers. Quite
of 10−4 (as used in BS13) compared with our fiducial value remarkably,theMRIchannelsextendoverseveralAUinthe
of10−3. Forthesakeofbrevity,wehererefrainfromvarying radial direction, although we remark that this may be overly
any of the many other input parameters of our model, as for pronounced in the axisymmetric case, where the growth of
instance,theX-rayorCRintensities,ortheFUVpenetration parasitic modes is likely to be artificially restricted (Pessah
depth,aswellasparametersaffectingthediskthermodynam- & Chan 2008). Because we use a net-vertical field, the lin-
ics. ear MRI modes appear more pronounced than in the other-
wisecomparable3DsimulationsofDzyurkevichetal.(2010)
3.1. Preliminaries withoutanetBz field.
Thepoloidalvelocityfieldplottedasblackarrowsismostly
Theinputstoourmodelsareverysimilartothoseincluded
disorderedintheturbulentregionsbutalsoshowssomelevel
in BS13, so it is instructive to compare the Elsasser number
ofcoherenceintheupperlayers,whereawindtopologycan
profilesinourmodelwiththeirfiducialmodelasameansof
beseen. Whilethewindisintermittentratherthansteadyand
understanding the similarities and differences between their
laminar, ageneraltendencyforoutflowisseen. Thisiscon-
results and ours. The fiducial model presented in BS13 is
sistentwithsimilarobservationsofdiskwindsbeinglaunched
computedat1au,andhasadust-to-gasmassfractionof10−4,
fromfully-MRI-activeaccretiondisks(see,e.g.,Suzuki&In-
sowecancomparethiswiththedot-dashedlinesinFigure1.
utsuka2009;Fromangetal.2013). Wequantifythemassloss
Theprofilesshownthere,andinfigure1ofBS13aresimilar,
with β decreasing below unity at disk heights |z| (cid:38) 4.5H, viatheverticaldisksurfacesbyevaluating
P
pturdevese.nItinngthMeirRiIntiutirablucloenncdeitiforonmanddevoeulrosp,iΛngOaitntchreesaesehsigmhoanltoi-- M˙wind ≡ 2π (cid:90) ro ρvθr sinθdr (cid:12)(cid:12)(cid:12)θ2 (7)
tonically from the midplane upward. The two initial states (cid:12)
r=ri θ=θ1
also have similar profiles for Λ , displaying local maxima
AD
at intermediate disk heights (z ∼ ±2.5H). The differences attheupperandlowerdisksurface,θ = θ1,θ2,respectively.
inthesurfacedensitiesintheBS13modelsandoursat1au, For model O-b6, we find a net mass loss rate of M˙ =
wind
combinedwithourinclusionofadirectX-raycomponentin 1.47±0.37×10−8M yr−1(alsocf. Table2below).
(cid:12)
additiontothescatteredcomponentmeansthatΛ ∼ 10at We observe the formation of field arcs and material flow-
AD
thislocalmaximuminourmodel,whereasitonlyrisesmod- ing downward along field lines towards the arc foot points.
estlyaboveunityinBS13.Wethereforeexpectthatourmodel This can be seen in the lower disk half around r = 2.8au,
withadust-to-gasmassfractionof10−4 shouldcontainnar- and r = 4.3au, for example, and illustrates the potential
rowregionsatintermediateheightsthatsupportthegrowthof roleofbuoyancyinstabilityinregulatingthedisk’smagneti-
MRIchannelmodes.ItisnoteworthythatBS13alsoobserved zation.Amplificationofthemagneticfieldthroughthegrowth
thedevelopmentoftheMRIatthebeginningoftheirsimula- of channel modes, and their break up through the action of
tions,butreportthatthesubsequentamplificationofthefield parasites, apparently leads to the formation of these locally
7
Figure2. VisualizationoftheglobaldiskstructureformodelO-b6withonlyOhmicresistivityandβp0 =106. Theazimuthalfieldsreachpeakstrengthsof
(cid:39)2.5Gintheactivedisklayers.Verticaloutflowisintermittent,butregionsofanorderedwindgeometrycanbeidentified(seeupperhalfaround3.5−4au).
uprisingfieldstructures. Asanalternativeexplanation,were-
mark that the dynamic evolution of such magnetic loops has
previouslybeenattributedtotheeffectofdifferentialrotation
ontothemagneticfootpoints(Romanovaetal.1998).
In contrast, within the magnetically decoupled midplane
layer, the magnetic field remains predominantly vertical, re-
taining the initial field configuration. Near the dead-zone
edge,theeffectofOhmicdiffusiongraduallydeclines. There,
the azimuthal magnetic field is allowed to change its value,
resulting in localized current layers. This is very similar to
thesituationobtainedinthedisksthatincludeambipolardif-
fusion,asdiscussedlater.
We illustrate the vertical disk structure in Fig. 3, where
we plot, in the upper panel, the total Reynolds and Maxwell
stressesrelativetotheinitialmidplanegaspressure,p . The
0
profilesshowthetypicalsignatureofalaminardeadzoneand
MRI-turbulent surface layers (Fleming & Stone 2003; Oishi
et al. 2007; Gressel et al. 2011, 2012), where the mechani-
cal and magnetic shear-stresses lead to transport of angular
momentum. Becauseoftherelativelyweaknet-verticalmag-
neticflux,thelevelofturbulenceismodest,andthevertically-
integrateddimensionlessMaxwellstressis(cid:39)8×10−5. With
the“viscous”massaccretionrate(i.e.,thateffectedbyinter-
nalstresses)approximatedby
Figure3. Verticalprofilesaveragedoverr ∈ [2,3]auandt = 25orbits
M˙ ≈ 3πΩ−1(cid:90) θ2 (cid:16)TReyn+TMaxw(cid:17) rdθ (cid:12)(cid:12)(cid:12) , (8) NfoergmatoivdeelvaOlu-bes6owfitthheOsthremsiscesre(usipsptievritpyanoenll)yaarendreaprmesiednptleadnebyβtphi0nl=ine1s.06.
visc Rφ Rφ (cid:12)
θ=θ1 (cid:12)r=rref
we estimate a value of M˙ (cid:39) 0.3 × 10−8M yr−1, at transportofangularmomentumalone. Thiscanbecompared
visc (cid:12)
r = 2.5au, if mass accretion were effected by turbulent to the actual mass accretion rate, which we can easily com-
ref
8
puteinourglobalmodelvia
(cid:12)
(cid:90) θ2 (cid:12)
M˙ ≡ 2π ρv r2sinθ dθ (cid:12) . (9)
accr r (cid:12)
θ=θ1 (cid:12)r=rref
Applying a radial average over r ∈ [2,3]au, and estimat-
ingfluctuationsoccurringwithinatimeintervalspanningt∈
[75,100]yr,wefindM˙ =(0.14±1.53)×10−8M yr−1
accr (cid:12)
(also see the last column of Table 2). Clearly, this diagnos-
ticturnsouttobesubjecttostrongfluctuationsformodelO-
b6–presumablyduetothepresenceofstrongradialmotion
within the MRI channel modes. Equation (9) will however
proveusefulwhenevaluatingtheradialdriftofmaterialinthe
localizedcurrentlayersseenintheAD-dominateddiskmod-
els. We note at this point that a simple energy conservation
argument prevents the steady-state mass loss rate (to infin-
ity)throughamagneto-centrifugalwindbeinglargerthanthe
accretionratethroughthedisk. ThelargevalueofM˙ re-
wind
portedaboverelativetoM˙ suggeststhatthewindlossrate
accr
is not yet converged. Indeed we note that BS13 report that
increasingtheverticalsizesoftheirshearingboxesleadstoa
reduction of the wind mass loss rates. A similar conclusion
isreachedbyFromangetal.(2013)forwindslaunchedfrom
Figure4. SameasFig.3,butformodelOA-b5includingbothOhmicresis-
turbulent disks. It appears that obtaining accurate and phys-
tivityandambipolardiffusion,andwithastrongernet-verticalfield.Notethe
ically meaningful estimates for the mass fluxes through the
differentscaleoftheabscissaintheupperpanel,reflectingthereducedlevel
upperboundariesofthediskwillrequiresimulationsthatare of“viscous”transport.Dottedlinesindicatethewindbaseatz=±zb.
trulyglobalintheverticaldirection,suchthatthefastmagne-
tosonicpointiscontainedwithinthesimulationdomainrather
ues of β . To test this, we have run simulations OA-b6, and
than outside of the boundary as it is at present for all of the p
models that we present in this paper (ensuring that the wind OA-b7,withβp0 =106,and107,respectively.
launching region is causally disconnected from the imposed
boundaryconditions). 3.3.1. Ambipolardiffusionwithweakverticalfields
In the lower panel of Fig. 3, the gas pressure, kinetic en- Inaccordancewiththeβ = ∞simulationsofBS13,we
p0
ergyandmagneticpressure,areshown. Withinthedead-zone findverylowlevelsofturbulenceinthesesimulations(cf.Ta-
layer,betweenz = ±4H,thediskissupportedbygaspres- ble 2). As a consequence, the radial mass transport via ac-
surealone,whichfollowsaGaussianprofile.Owingtothead- cretion (see M˙ values) is found to be negligible. At the
accr
ditionalmagneticpressuresupportwithintheactivelayer,the
same time, because of the weak vertical field, the magneto-
effective scale height increases roughly twofold there. Note
centrifugal wind is absent in model OA-b7, and only poorly
that unlike seen in the isothermal simulations of BS13, cf.
developed (and intermittent) in model OA-b6. We conclude
their fig. 3, the magnetic pressure does not significantly ex-
that the corresponding field amplitude of about 10mG (at
ceed the value of the gas pressure in large parts of the disk
1au) can be assumed as a minimal requirement for signif-
corona. Weattributethisdifferencetothedissipationheating
icant mass transport by magnetic effects. In the following,
(Hirose&Turner2011)presentinoursimulations,whichwe
weaccordinglyfocusonmodelsderivedfromourfiducialrun
note,however,maybeoverlypronouncedintheaxisymmet- OA-b5,withβ =105,thatis,B =31.5mG(4.2mG)at
p0 z0
ricmodels. Thekineticenergyonlybecomesimportantvery
1au(5au).
closetotheverticalboundariesofourmodel. Thisisalsore-
flectedinthepositionoftheAlfve´npoint,whichisrelatively
3.3.2. Thefiducialmodel
poorly constrained and which we infer as (7.60±0.45)H.
Values for these vertical points are listed in Table 2 for all InFigure4,weshowverticalprofilesoftheRφstresscom-
models. Wereporttime-andspace-averagedvaluesobtained ponents for our fiducial model, OA-b5, where owing to the
forr ∈[1,5]au,andwehaveappropriatelytakenintoaccount stronger net vertical field the MRI is ultimately suppressed
the flaring of the disk. Because of the turbulent character of byambipolardiffusionandalaminarwindsolutionisestab-
the upper disk layers, we cannot obtain a good estimate for lished instead. Dashed vertical lines indicate the base of the
thelocationofthebaseofthewindinmodelO-b6. wind,whichwedefineaccordingtoWardle&Ko¨nigl(1993)
as the vertical position, z = z , where the rotation becomes
b
3.3. Laminardiskmodelswithsurfacewinds super-Keplerian. Since we define our azimuthal velocity, vφ
as the deviation from the Keplerian background flow, this
We now discuss models that include both Ohmic resistiv-
can be verified in Fig. 4 by the fact that the Reynolds stress,
ity and ambipolar diffusion. A key finding of BS13 is that
theinclusionofambipolardiffusioninhibitslineargrowthof TRRφeyn ≡ρvRvφ,vanishesatthispoint. Thesameistruefor
the MRI in the intermediate disk layers, that is, the regions theverticaltensorcomponent,TReyn,whichmakesitconve-
zφ
thatarenotaffectedbyOhmicresistivity. Sincetheambipo-
nienttodefinethewindstress,
lardiffusioncoefficientdependsontheplasmaparameter,one
mightexpectthattheeffectofADislesssevereforhighval- Tzφ ≡ TzMφaxw(cid:12)(cid:12)z=+zb − TzMφaxw(cid:12)(cid:12)z=−zb (10)
9
at this point. For the Rφ component of the stress tensor, re-
sponsible for redistributing angular momentum radially, we
furthermore obtain average values within the disk body z ∈
[−z ,+z ],
b b
1 (cid:90) +zb(cid:16) (cid:17)
T ≡ TMaxw+TReyn dz, (11)
Rφ 2z Rφ Rφ
b −zb
which we list separately in Table 2, where values have been
normalizedtothemidplanepressure, p . Alreadyforamid-
0
planeplasmaparameteraslowas105,thewindstressexceeds
theRφ(accretion)stressbyafactoroffour. BS13alsoreport
thatthewindstressexceedstheaccretionstress;fortheirfidu-
cialrun,thediscrepancyisbiggerthananorderofmagnitude.
InthelowerpanelofFigure4,weseethatintheabsenceof
MRIturbulence,thehydrostaticbalanceextendsfurtherupin
the disk and the gas pressure remains the dominant stabiliz-
ingforceallthewayuptothebaseofthewind. Eveninthe
magneticallydominatedwindlayer,magneticpressuregradi-
entsremainmoderateandthescaleheightofthediskremains
roughlyconstantuptoz .
b
InFig.5,weshowazoom-inoftheinnerpartoftheglobal
field topology in our fiducial simulation OA-b5. The up-
per panel shows the resulting field configuration early on,
whenresistivelymodifiedMRIeigenmodesappearinginside
r = 1.5au are still clearly discernible. Localized channel
solutions form at z ≈ 3H near the interface between the
magnetically decoupled region, characterized by Λ (cid:28) 1
O
and purely vertical field lines with B ≈ 0, and the AD-
φ
dominated intermediate layer. In this small region, both El-
sassernumbersareaboveunity(cf.Fig.1),allowingforlinear
growthofMRIchannelmodes.Thistransitionlayerischarac-
terizedbyazigzag-shapedabruptchangeinthefieldlines(see
Fig. 8 below), associated with strong current sheets. At this
point,thelaminarwindsolutionisalreadywellestablishedin
theFUVionizedsurfacelayer,andhasspreadthroughoutthe
radialextentofthesimulationdomain.
While the disk wind itself has already reached a steady
state at this point in time, the strong field layers continue to
evolve. In the lower panel of Fig. 5, we accordingly show a
laterevolutionstage,wheretherelicMRIchannelshavemor-
phed into strong azimuthal field belts. These field belts can
be understood as “undead zones”, that is, a region that does
notsustainMRI,butwherediffusioncanbebalancedbythe
winding-upofpreexistingradialfieldviadifferentialrotation Figure5. Field topology of our fiducial simulation at different evolution
times. The azimuthal magnetic field (color) has been restricted to values
(Turner&Sano2008). Itappearsthattheinitialamplification
|Bφ| < 125mGforclarity;peakvaluesareafewhundredmG. Wealso
of the magnetic field through the growth of channel modes showprojectedmagneticfieldlines(white)andvelocityvectors(black).Ad-
causes theambipolar diffusion coefficient toincrease until it ditionallinesindicatetheposition,zb,ofthewindbase(dot-dash),andthe
radiallocationoftheprofilesplottedinFigs.6and7(dashedlines).
quenchesfurtherdevelopmentoftheMRI.Thisnotionissup-
ported by a reference simulation, where the diffusion coeffi-
cients η and η were held fixed at their initial value, and
O AD
where the MRI continues to produce quasi-turbulent fields supportedbyourpresentmodels,wherecurrentlayersappear
akintotheonesseeninFig.2fortheOhmic-onlycase. The to emerge from inner disk regions. What determines the ex-
reversalofthefielddirectionseenatr (cid:39) 1.5auisarelicof act shape of the surviving field configuration at late times is
the local nonlinear development of the MRI modes early on presentlyunclear. Wespeculatethattheparticularrealization
(as may be seen in the upper panel at r (cid:39) 0.75au). This seeninthelowerpanelofFig.5isnotnecessarilyauniqueso-
interfacebetweentheazimuthalfieldbeltsofoppositeparity lutiongiventhespecificdiskgeometryandionizationmodel,
movesradiallyoutwardataspeedof≈0.01auyr−1. butmaytosomeextenddependontheearlyevolutionhistory.
The emerging current layers are directly related to However, simulations that were initiated with different ran-
resistively-modified vertically-global MRI channel modes – domseeds,butwereotherwiseidentical, onlyshowedminor
seeLatteretal.(2010),andFig.8inSection3.5. WithMRI variations in the final field appearance. The configurations
channels potentially spanning large radial extents, as in the observed at late times are moreover found to remain quasi-
OhmicrunshowninFig.2,suchlayersmaybefeddiffusively stationary,atleastduringtheevolutiontimeof100yrthatwe
with magnetic field from the MRI-active inner disk. This is haveinvestigated.
10
Table2
Summaryofsimulationresults.
zb zA TRRφeyn TRMφaxw TzMφaxw M˙wind M˙accr
[H] [H] [10−6p0] [10−5p0] [10−5p0] [10−8M(cid:12)yr−1] [10−8M(cid:12)yr−1]
O-b6 — 7.60±0.45 6.87±14.4 7.44±0.95 — 1.58±0.59 0.04±1.91
OA-b5 5.23±0.07 7.31±0.17 3.63±0.19 2.22±0.06 9.82±0.08 0.79±0.01 0.43±0.01
OA-b6 7.22±0.48 6.94±0.37 −0.21±0.18 0.79±0.06 0.58±0.06 0.19±0.06 0.09±0.02
OA-b7 7.31±0.70 6.39±0.16 0.07±0.13 <0.01 0.22±0.01 0.03±0.01 0.00±0.03
OA-b5-d4 5.27±0.07 7.33±0.18 0.11±0.30 2.88±0.11 9.88±0.12 0.76±0.02 0.34±0.04
OA-b5-flr 4.81±0.03 6.90±0.31 0.26±0.21 1.78±0.02 14.3±0.02 1.57±0.01 0.64±0.02
OA-b5-flr-nx 4.78±0.03 7.50±0.30 2.28±9.24 1.87±0.04 13.0±0.04 1.58±0.01 0.64±0.03
OA-b5-nx 5.10±0.04 7.34±0.13 0.94±9.29 1.89±0.06 7.84±0.02 0.67±0.01 0.30±0.04
Theverticalpositionofthebaseofthewind,zb,andtheAlfve´npoint,zA,arefoundindependentontheradiallocationwhenmeasuredinlocalscale
heights,H. TheviscousaccretionstressesTRφareverticallyintegratedwithin|z| ≤ zb–notethedifferentunitsforReynoldsandMaxwellstresses.
Thewindstress,Tzφ,ismeasuredatz=±zb.Allstressesdependweaklyonradius;listedvaluesareatr=3au.
3.4. Typicalwindsolutions inner disk, probably because the vertical field lines become
too stiff (owing to B increasing with radius) to be suitable
PerhapsthemostnoteworthyresultfromrunOA-b5isthe z
forwindlaunchingintheouterdisk.
spontaneousemergenceoftheexpected“hourglass”geometry
Foraneutralsituationwith∂ B =0,westillobserveout-
for the magnetic field, allowing magnetic tension forces to R z
wardbendingofthefieldlines. Startingfromtheinnerradial
acceleratethewind.Fortheadoptedverticalfieldpolarity,the
domain, the outward propagating establishment of the wind
fieldmustbendradiallyoutwardneartheupperandlowerdisk
regionisfoundtostallatsomeradius,whereasthewindwas
surfaces, and the azimuthal field should point in the positive
quickly established throughout the entire domain in model
directioninthelowerhemisphereandreversedirectionabove
OA-b5. Asbefore(forthecaseofanoutwardmagneticpres-
themidplane,asobservedintheupperpanelofFigure5.
sure gradient), this is probably related to the field lines be-
comingtoostrongtosupportthewindlaunchingmechanism
3.4.1. Robustnessoftheemergingwindgeometry
at large radii. The wind indeed stalls further out in a run
The shearing box simulations of BS13 have reflectional with weaker overall net-vertical field. This poses the ques-
symmetryintheradialdirection, thatis, withrespecttomir- tionwhethertheverticalprofilesofΛ (z), andΛ (z)that
O AD
roringaboutx = 0, andhencecontainnoinformationabout we obtain from our ionization model put strong constraints
the location of the star. Instead of giving rise to a physical on permissible vertical field amplitudes as a requirement for
wind solution, their fiducial run produces a radial field that windlaunching.
possessesodd-symmetryaboutthemidplane,withaninward Apossiblereasonfortheoutwardbending,eveninthecase
pointingfieldinonehemisphereandanoutwardpointingone of neutral magnetic pressure forces, may be that the vertical
in the other. To test the robustness of the emerging solution shearpresentintheglobalmodels(becauseofthebaroclinic
in our setup, we initiate several instances of run OA-b5, us- conditions) naturally bends the azimuthal component of the
ingdifferentrandomnumberseedswhenperturbingtheinitial field in the correct direction required by the physical wind
velocityfield,andineachcasewerecoverthecorrectwindso- solution. We have checked that the outward bending of the
lutionwithoutwardbendingfieldlines. Ourinitialmodelhas vertical field lines is however equally seen in (strongly flar-
anet-verticalfieldwitharadialdependencesuchthatthemid- ing)globallyisothermalmodelsthatdonothaveanyvertical
planevalueofβp remainsconstant. Becausethegaspressure gradientsintherotationvelocityoranyradialgradientsinthe
itself decreases with radius, the corresponding Bz(R) (with vertical magnetic field. The outward pointing configuration
∂RBz < 0) results in an outward magnetic pressure force, isenergeticallyfavoredbecausesettingitupinvolvesdiluting
whichmoreoverhasastrongeraccelerationeffectinthelow- ratherthanconcentratingmagneticflux.
density upper disk layers. While a vertical flux distribution Oursimulationsdemonstrateinanycasethatglobalmodels
thathasthefluxcentrallyconcentratedcanbeexpectedwhen spontaneously develop the correct field geometry, although
accountingforinwardadvectionandoutwarddiffusionofflux we caution that this conclusion needs to be tested in future
(Okuzumi et al. 2014), our particular choice is of course ar- simulations that moreover adopt improved boundary condi-
bitrary. Radialmagneticgradientsarefurthermoreknownto tions for the magnetic field at the upper and lower (and po-
affecttheoutflow’smassloadinginaxisymmetriccalculations tentially the inner radial) surfaces of the simulation domain.
wherethediskisaboundarycondition(Pudritzetal.2006). Ultimately,theemergenceofthewindgeometrywillhaveto
In our simulations, the pressure gradient related to Bz(R) bestudiedinsimulationsthatdonotstartfromwell-motivated
may potentially affect the direction toward which the field (but nevertheless arbitrary) initial conditions, but do account
lines first bend from their starting configuration. We have for the formation of the PPD from a collapsing molecular
tested this by running reference models with different mag- cloud(Lietal.2014).
netic pressure gradients. If ∂RBz > 0, the initial condition Lastly,fortheparametersthatwehaveconsideredhere,we
has an unbalanced inward pressure force, and we do indeed do not find a self-limiting of the wind via shielding of FUV
find the field lines to spontaneously bend inward; this hap- photons (Bans & Ko¨nigl 2012). This is consistent with the
penssimultaneouslyintheupperandlowerhemisphereofthe T Tauri scenario of Panoglou et al. (2012), who performed
disk,suchthattheoverallsymmetrythatisobtainedisagain disk wind chemical modeling for various protostellar evolu-
even. The launching of the inward wind is restricted to the