Global Evolution of an Accretion Disk with Net Vertical Field: Coronal Accretion, Flux Transport, and Disk Winds

8 MB·

by  Zhaohuan Zhu

Draft version January 18, 2017
PreprinttypesetusingLATEXstyleemulateapjv.5/2/11
GLOBAL EVOLUTION OF AN ACCRETION DISK WITH NET VERTICAL FIELD: CORONAL ACCRETION,
FLUX TRANSPORT, AND DISK WINDS
Zhaohuan Zhu1, AND James M. Stone2
Draft version January 18, 2017
ABSTRACT
We report results from new global ideal MHD simulations that study thin accretion disks (with
thermal scale height H/R=0.1 and 0.05) threaded by net vertical magnetic fields. Our computations
spanthreeordersofmagnitudeinradius,extendallthewaytothepoletopreventthelossofmagnetic
7
fields, and are evolved for more than one thousand orbits at the inner edge of the disk, corresponding
1
to more than one viscous time over the inner decade in radius. In order to properly resolve the
0
magnetorotational instability (MRI) at all radii, static mesh refinement is used. We find that: (1)
2
inwardaccretionoccursmostlyintheuppermagneticallydominatedregionsofthedisk,similartothe

n ”coronal accretion” flows found in previous GRMHD simulations. Rapid inflow in the upper layers
a combined with slow outflow at the equator creates strong Rφ and zφ stresses in the mean field; the
J
vertically integrated Rφ stress α ∼ 0.5−1 when the initial field has β = 103 at the midplane. (2)
0
7 A quasi-static global field geometry is established in which flux transport by inflows at the surface
1 is balanced by turbulent diffusion. The resulting field is strongly pinched inwards at the surface.

A steady-state advection-diffusion model, with turbulent magnetic Prandtl number of order unity,
]
P reproduces this geometry well. (3) Weak unsteady disk winds are launched beyond the disk corona
at z/R ∼ 1 with the Alfv´en radius R /R ∼ 3. The wind crosses all three critical points within our
E A 0
computational domain. Although the wind is episodic, we show the time averaged properties are well
.
h described by steady wind theory. We measure the torque exerted on the disk by the wind; even with
p strongfields,β =103 atthemidplaneinitially,only5%oftheangularmomentumtransportisdriven
0
- by torque from the wind, and the wind mass flux from the inner decade of radius is only ∼ 0.4% of
o
themassaccretionrate. Withweakerfields(largerβ initially)orthinnerdisks, thewindcontributes
r 0
t even less to accretion and mass loss. (4) Within the disk, most of the accretion is driven by the Rφ
s
stressfromtheMRIandglobalmagneticfields. Oursimulationsmayshedlightonthelargevaluesof
a
[ α inferred in cataclysmic variables, high near-infrared flux in Herbig stars, fast inflow in transitional
disks, outflow in FU Orionis systems, and the transport of chondrites in our solar system.

1
Subject headings: accretion, accretion disks - dynamo - diffusion - magnetohydrodynamics (MHD) -
v
instabilities - turbulence - protoplanetary disks -
7
2
6
1. INTRODUCTION dle & Koenigl 1993; K¨onigl et al. 2010; Salmeron et al.
4
2011) is also important in these systems.
0 Accretion and outflow have been observed in a wide
Outflow/wind/jet launching mechanisms require large
. rangeofastrophysicalsystems,rangingfromprotoplane-
1 scalepoloidalmagneticfields. Inthemagnetocentrifugal
tary disks (Hartmann 1998) to supermassive black holes
0 wind model, the magnetic field is anchored in the rotat-
(Begelman et al. 1984). Two mechanisms are likely to
7 ingdisk. Ifthepoloidalcomponentofthefieldmakesan
drive accretion in most systems: MHD turbulence or a
1 angle of less than 60o with the disk surface, the centrifu-
magnetized disk wind. Turbulence can lead to net R−φ
: gal force can break the potential barrier and accelerates
v stresswithinthediskthattransportsangularmomentum
matter outwards, leading to outflow. The stress at the
i radially, while a disk wind can lead to a z−φ stress at
X wind base launches the outflow and, at the same time
the disk surface which carries angular momentum away
r vertically. Itisthoughtthemainmechanismdrivingtur- torques the disk, leading to accretion below the wind
a region. Such a picture has been confirmed by numerical
bulence is the magnetorotational instability (MRI, Bal-
simulationswithprescribedpoloidalmagneticfields(e.g.
bus & Hawley 1991, 1998), while in a Newtonian poten-
the review by Pudritz et al. 2007).
tial the main mechanism which produces a wind is the
Netverticalmagneticfieldsalsoplayanimportantrole
magnetocentrifugal effect associated with vertical fields
inMRIturbulence. Itisknownthattheturbulentstress
(Blandford & Payne 1982). Since many astrophysical
increases as the net field increases (Hawley et al. 1995),
disks are poorly ionized (e.g. protoplanetary disks), the
thus net field promotes both outflow and disk accretion.
effects of non-ideal MHD on MRI (e.g. the review by
Toaddresstherelativeimportanceofturbulenceanda
Turner et al. 2014b) and disk winds (Konigl 1989; War-
diskwindindrivingaccretion,weneedtorelyonnumeri-
calsimulations. However,simultaneouslyresolvingsmall
[email protected],[email protected]
scale turbulence and capturing the large scale disk wind
1Department of Physics and Astronomy, University of
Nevada, Las Vegas, 4505 South Maryland Parkway, Las Vegas, is challenging. In order to resolve the MRI, most previ-
NV89154,USA ous simulations with net vertical flux only study a small
2Department of Astrophysical Sciences, 4 Ivy Lane, Peyton patch of the disk using the shearing box approximation
Hall,PrincetonUniversity,Princeton,NJ08544,USA
2 Zhu & Stone
(Suzuki & Inutsuka 2009; Bai & Stone 2013; Fromang & Norman 1994) which observed that magnetic flux is
et al. 2013). However, the resulting wind in such simu- mainly transported in the corona of the disk and mag-
lations can flow in either radial direction (Bai & Stone netic fields are pinched within the corona. Thus, the
2013; Lesur et al. 2014), and the outflow rate drops dra- above 1-D estimates may not apply in more complex
maticallywhenatallerboxhasbeenused(Fromangetal. multi-dimensional field geometries. Since the radial dis-
2013). Globalsimulationsofturbulencewithinanaccre- tribution of global magnetic fields also determines the
tion disk with net vertical fields have been carried out collimation of disk winds (Anderson et al. 2005; Pudritz
recently (Suzuki & Inutsuka 2014; Gressel et al. 2015) et al. 2006), the processes which determine and main-
and strong outflows and high disk accretion rates have tain large scale magnetic fields in disks are essential for
been observed. However, these simulations have lim- sustaining outflow or even accretion.
ited radial range, only cover a wedge around the disk In §2, the theoretical framework on describing turbu-
midplane z ∈ [−0.5R,0.5R], and have not been evolved lenceanddiskwindispresented. Ournumericalmethod
long enough for a true steady-state to emerge. On the isintroducedin§3. Theresultsarepresentedin§4. After
otherhand,global(usuallytwo-dimensional)simulations a short discussion in §5, we conclude in §6.
ofdiskwindsthatdonotcapturetheMRIinthediskin-
teriorhavebeenreported(Stone&Norman1994;Ouyed 2. THEORETICALFRAMEWORK
& Pudritz 1997a; Krasnopolsky et al. 1999, 2003; Kato
Which is more important for disk accretion: turbu-
etal.2002;Porth&Fendt2010;Ramsey&Clarke2011).
lence or outflow? If we average the angular momentum
In such work, either only the wind region has been stud-
equation in the azimuthal direction and integrate it in
iedorexplicitresistivityhasbeenassumedinthediskre-
the vertical direction, we can derive
gion as a ”sub-grid” of the MRI (e.g. Casse & Keppens
2002, 2004; Fendt & Cˇemelji´c 2002; Zanni et al. 2007; ∂(cid:82)(cid:104)ρv (cid:105)dz 1 ∂ (cid:18) (cid:90) (cid:19)
φ =− R2 ((cid:104)ρv δv (cid:105)−(cid:104)B B (cid:105))dz
Tzeferacos et al. 2009, 2013). In this work, we carry out ∂t R2∂R R φ R φ
globalthree-dimensionalsimulationswhichcaptureboth
tmheeshM-rRefiInaenmdendtisaknwdisnpdecsieallf-pcoonlasrisbteonutnlydabryy cuosnindgitiboontsh. − 2π1R2∂Rv∂kRM˙acc −((cid:104)ρvzvφ(cid:105)−(cid:104)BzBφ(cid:105))(cid:12)(cid:12)(cid:12)(cid:12)zmax
An important question that can only be addressed in zmin
(2)
global simulations is the net rate of transport of vertical
magnetic field. Net vertical fields are both advected by or
the large-scale flow, and diffuse due to turbulence. In ∂(cid:82)(cid:104)ρδv (cid:105)dz 1 ∂ (cid:18) (cid:90) (cid:19)
the phenomenological model by van Ballegooijen (1989) φ =− R2 ((cid:104)ρv δv (cid:105)−(cid:104)B B (cid:105))dz
diskaccretionandfielddiffusionismodeledusingturbu- ∂t R2∂R R φ R φ
lseenntsivtiivsceolysitoynatnhderreeslaisttiivveitaym, hpoliwtuevdeerotfhbeortehs.ulAtsdecpleeanrdlys − 2Mπ˙aRcc2∂∂RRvk −((cid:104)ρvzδvφ(cid:105)−(cid:104)BzBφ(cid:105))(cid:12)(cid:12)(cid:12)(cid:12)zmax
summarizedinGuan&Gammie(2009), theevolutionof zmin
the poloidal fields is governed by (3)
η for the perturbed quantities (δv = v −v ). The sym-
∂ A =−v B − ∂ B +η∂ B φ φ k
t φ R z R z R R z bol (cid:104)(cid:105) denotes that the quantity has been averaged over
1 (cid:20)1 (cid:21) the φ direction 3, and M˙ = 2πR(cid:82) ρv dz. We refer
∂ A =v ∂ (RA )+η∂2A +η∂ ∂ (RA ) . acc R
t φ RR R φ z φ R R R φ to (cid:104)ρvRδvφ(cid:105) and −(cid:104)BRBφ(cid:105) as the radial Reynolds and
(1) Maxwellstressrespectively,todistinguishthemfromthe
vertical stresses T . If we normalize the stresses with
zφ
if only radial motion is considered. In Equation 1, A pressure, we can identify different contributions to the
φ
is the azimuthal component of the vector potential. On total stress α as
therighthandside,thefirsttermisadvectionbytheac-
cretion flow. Since vR ∼ ν/R in the viscous model, this αRφ,Rey =(cid:104)ρvRδvφ(cid:105)/(cid:104)p(cid:105) and αRφ,Max =−(cid:104)BRBφ(cid:105)/(cid:104)p(cid:105)
term is roughly νBz/R. The second term is the vertical αzφ,Rey =(cid:104)ρvzδvφ(cid:105)/(cid:104)p(cid:105) and αzφ,Max =−(cid:104)BzBφ(cid:105)/(cid:104)p(cid:105).
diffusion of radial field and is roughly ηB /z. If we as-
R
sume the large scale field enters the disk at an angle of Stresses and the α parameters in spherical-polar coordi-
∼45o,asinthediskwindmodel,wehaveB /z ∼B /H. nate can be defined in similar ways. Since the vertically
R z
Thethirdtermisradialdiffusion,roughlyηB /R. Thus, integratedR−φstressdeterminesthediskaccretionrate
z
the first and second terms dominate, and they balance as in Equation 3, we can define the vertically integrated
each other when ν/η ∼ R/H or the turbulent magnetic α parameter as
Prandtl number (Pr) ∼ R/H. However, local shear-
(cid:82)
T dz
ing box simulations suggest that Pr∼1 (Guan & Gam- Rφ
α = . (4)
mie 2009; Lesur & Longaretti 2009; Fromang & Stone int Σc2
s
2009), implying that large-scale fields will diffuse out-
wards faster than the inward advection (Lubow et al. where TRφ is the sum of both radial Reynolds and
1994). To maintain large scale magnetic fields, either Maxwell stress. If we choose vk so that (cid:104)ρδvφ(cid:105)=0 and
turbulentdiffusionissignificantlyreduced(Spruit&Uz- assume that the magnetic stress dominates at the disk
densky 2005; Rothstein & Lovelace 2008) or inward ac-
cretion is increased. An important insight is provided 3 vk hasbeenassumedtobeconstantalongz. Withoutthisas-
by Beckwith et al. (2009) (see also simulations by Stone sumption,therewillbeanadditionaltermrelatedtoM˙loss∂vk/∂z.
Global MHD Simulations 3
Fig. 1.—Densityiso-surface(leftside)andpoloidalsliceofdensityshowingthenestedgrid(rightside)fromasnapshotattime42T0 in
ourfiducialmodelwithinitialβ0 =1000. StaticmeshrefinementisusedtocaptureMRIturbulenceatthediskmidplane,andtoextend
thecalculationtothepole.
surface, Equation 3 can be written as enables mesh refinement, allowing global numerical sim-
M˙acc =−∂Rv2π/∂R(cid:32)∂∂R(cid:0)R2αRφ,intΣc2s(cid:1)−R2(cid:104)BzBφ(cid:105)(cid:12)(cid:12)(cid:12)(cid:12)zmax(cid:33)ugcyelaolitmniodenrtisrciascplsaaonnundricnsepgtheaermrlaicsragilen-prcoauldraivraillciornaoenradgriecn.oaFotreudsr)itnhaaertreemscoa(rreee.g,fu.tlhilnye
k zmin implemented so that angular momentum is conserved to
(5)
machineprecision,acrucialfeaturetoenabletheangular
Equation 5 suggests that disk accretion is due to both
momentum budget analysis as presented in §4.1.
T within the disk and T exerted at the disk surface.
Rφ zφ Since a disk wind flows radially, we adopt a spherical-
Normally, the R−φ stress is from disk turbulence, and
polar coordinate system (r, θ, φ) for our simulations,
thez−φstressatthesurfaceisfromamagnetocentrifu-
which should minimize the effects of the domain bound-
gal disk wind. If both the R−φ and z−φ stresses have
ary on the wind properties (Ustyugova et al. 1999). We
similarvalues,thesecondtermontherightwillbelarger
implement a special polar boundary in the θ direction
thanthefirsttermbyafactorofR/z whichcanbequite
allowing our simulation domain to extend all the way to
large for a thin disk . On the other hand, with vigor-
the pole, which is different from previous similar simu-
ousturbulencetheinternalstressmaybelargerthanthe
lations where a hole was carved out close to the pole.
surfacestress,thusitisunclearifdiskaccretionisdriven
This boundary condition prevents the loss of magnetic
by turbulence or wind when both processes are present.
fields at the previously carved-out pole, which is crucial
3. METHOD tostudytransportofmagneticfieldsindisks. Thedetails
on the implementation are given in the Appendix.
Wesolvethemagnetohydrodynamic(MHD)equations
Although we adopt spherical-polar coordinates for the
in the ideal MHD limit using Athena++ (Stone et al.
simulations, we transform some quantities to cylindrical
2017, in preparation). Athena++ is a newly developed
coordinates to study physical processes in disks. In this
gridbasedcodeusingahigher-orderGodunovschemefor
paper, weuse(R, φ, z)todenotepositionsincylindrical
MHD and the constrained transport (CT) to conserve
coordinatesand(r, θ, φ)todenotepositionsinspherical
the divergence-free property for magnetic fields. Com-
polar coordinates. In both coordinate systems φ repre-
pared with its predecessor Athena (Gardiner & Stone
sents the azimuthal direction (the direction of disk rota-
2005, 2008; Stone et al. 2008), Athena++ is highly op-
tion). Our simulation domain extends from r = 0.1
timized for speed and uses a flexible grid structure that min
4 Zhu & Stone
to r = 100 with θ from 0 to π. φ extends from 0 to We choose p = −2.25 and q = −1/2 in our simula-
max
2π except for the thin disk case for which θ is from 0 to tion so that the disk surface density Σ ∝ R−1. Then
π/4. ρ (R ,z = 0) is 1. The time unit in the simulation is
0 0
The initial density profile at the disk midplane is T = 2π/Ω(R = 1). The grid is uniformly spaced in
0
(cid:18) R (cid:19)p ln(r), θ, φ with 136×64×128 grid cells in the domain
ρ (R,z =0)=ρ (R ,z =0) , (6) of [ln(0.1), ln(100)]×[0, π ]×[0, 2π] at the root level.
0 0 0 R0 Thus, ∆r/r equals 0.052 at the root level. We use an
open boundary condition in the radial direction, the po-
whilethetemperatureisassumedtobeconstantoncylin-
lar boundary condition in the θ direction, and periodic
ders
(cid:18) R (cid:19)q boundary conditions in the φ direction.
T(R,z)=T(R ) . (7) We have carried out three simulations, one with
0 R0 (H/R)R=R0 =0.1andinitialplasmaβ=1000atthedisk
midplane (the fiducial case), one with (H/R) =0.1
Hydrostatic equilibrium in the R − z plane requires R=R0
that (e.g. Nelson et al. 2013) and β = 104 (the weaker field case), and one with
(H/R) = 0.05 and β = 103 (the thin disk case).
(cid:20)GM (cid:18) 1 1(cid:19)(cid:21) For theR=siRm0ulations with (H/R) = 0.1, three lev-
ρ0(R,z)=ρ0(R,z =0)exp c2 √R2+z2 − R , els of refinement have been adoRp=teRd0 towards the disk
s
(8) midplane, and ∆r/r at the midplane is 0.0065. Since
and H/R=0.1 at R = 1, the disk scale height is resolved by
15 grids at the finest level. For the thin disk case, 4 lev-
(cid:34) (cid:18) c (cid:19)2 qR (cid:35)1/2 els of mesh refinement have been used so that the disk
vφ(R,z)=vK (p+q) vφ,sK +1+q− √R2+z2 , stacatiloenhaeligchotsti,sφstiinll rtehseoltvheidnbdyisk15cagsreidos.nlTyoesxatveendcsomfrpoum-
(9) 0 to π/4. The (H/R) = 0.1 and β = 1000 case is
(cid:112) R=R0 0
where c = p/ρ is the isothermal sound speed, v = ourfiducialmodel,sincetheMRIisbetterresolvedatthe
s K
(cid:112) midplaneandthesimulationrunslonger. Wewillmainly
GM /R,andH =c /Ω . Alocalisothermalequation
∗ s K
present results for our fiducial case, but other cases will
ofstateisassumedduringthesimulation4. UsingEqua-
be discussed briefly in §4.4.
tions 8 and 9, the density and velocity will become infi-
In ideal MHD, the wavelength for the fastest growing
niteatthepole. Toavoidthis,weuseR=max(R,r )
min linear MRI mode satisfies
in the above equations.
(cid:112)
We adopt a density floor which varies with position λ=2π 16/15|v |/Ω=9.17β−1/2H (13)
A
ρfl =ρρρffflll,,,000(cid:16)(cid:16)(cid:16)rRRRRmR000i(cid:17)(cid:17)npp(cid:17)(cid:0)(cid:0)pzz(cid:0)1122z(cid:1)(cid:1)12i(cid:16)f(cid:1)5if−Rif2<Rr−Rrr≥mrmi<nrinmra(cid:17)imnn(cid:16)idn4rarmn<rimdni3−rnrR>+3r1m(cid:17)in tcf(ohHerlaalsttwhtaleheteyfiRmdeu=tocsai1tal..lu1Tnc9asht9seae5b)da.lenendWwsia3tivty×heflle1βon00og−=rt7h1p0faio0srra0rm,tehsλeoetleivosrethd∼ρefbrl0y,0c.24ai9s.s5e1Hsg0.r−sido6
min min
(10)
We choose this density floor so that it becomes small at
thediskatmosphereandbecomeslargeclosetotheinner
radial boundary.
Initial magnetic field is assumed to be vertical. To
maintain ∇·B = 0, we use the vector potential A to
initialize the magnetic fields. We set
 (cid:16) (cid:17)m
 21 ×R×B0 rmR0in if R≤rmin Fig. 2.— ThequalityfactorforBθ andBφinourfiducialmodel
A = att=42T0.
φ
 RBm0 Rmm++21 + BR0rmmmRi+n2(12 − m1+2) if R>rmin To demonstrate that our fiducial case has reached the
0 0
(11) necessary resolution to capture MRI, we plot the az-
where m=(p+q)/2. Thus, the vertical magnetic fields imuthally averaged quality factor ((cid:104)Q (cid:105) and (cid:104)Q (cid:105)) for
θ φ
are B and B in Figure 2. The quality factor is defined
 (cid:16) (cid:17)m θ φ
B0 rmR0in if R≤rmin agsrotwhiengnummobdeer(oNfogbrliedectelalsl.t2h0a1t0)r,esolve the fastest MRI
B = (12)
z
B0(cid:16)RR0(cid:17)m if R>rmin. Q =Qλθ =/(λrsθi/n(θr∆∆φθ)) ((1145))
φ φ
With this setup, the plasma β at the disk midplane be- (cid:112)
where λ = 2π 16/15|v |/Ω and λ =
yond r is a constant. θ A,θ K φ
min (cid:112)
2π 16/15|v |/Ω . The quantities v and v
A,φ K A,θ A,φ
4 To achieve this, we have used the adiabatic equation of state are the Alfven velocity calculated using the θ and
withγ=1.4,butinstantaneouscoolingisappliedateachtimestep. φ components of the magnetic field, and Ω is the
K
Global MHD Simulations 5
Keplerian angular velocity at the disk midplane. To get 4.1. Accretion Structure
converged results for MRI turbulence, Sorathia et al.
Although the whole disk accretes inwards, the accre-
(2012) have shown that if Q ≈ 10, Q needs to be
φ z tion flow has a complicated structure within the disk as
(cid:38) 10−15, and if Q (cid:38) 25, while Q can be smaller (∼
φ z showninFigure4. Thegreencurvesinthefigurearethe
6 in their Figure 8). As shown in Figure 2, the quality
velocity field lines and the white curves are where the
factor decreases towards the disk midplane despite the
plasma(cid:104)β(cid:105)=1. Themoststrikingfeatureisthatthedisk
fact that it has been boosted by a factor of 8 towards
accretes through the disk atmosphere where the disk is
the midplane due to mesh refinement. Q is larger than
φ magnetically dominated (β < 1) as shown in the right-
10 in the whole disk (>50 at R = 1), and Q is around
θ mostpanel. Suchsurfaceaccretionhasbeenseenasearly
10 at R=1. Thus, our simulation should have adequate
as Stone & Norman (1994) and recently been studied in
resolution to capture MRI turbulence except at the
GRMHD simulations by Beckwith et al. (2009) where it
midplane of the inner region (R<0.3).
was termed ”coronal accretion”.
The fast coronal accretion at the disk surface carries
4. RESULTS themagneticfieldsinwards,thuspinchingmagneticfields
After running for t=45.6 T , which is equivalent to atthedisksurface. Themagneticfieldlinesareshownin
0
1442 orbits at the inner boundary, our fiducial disk Figure 5. The field lines change from the vertical direc-
reaches a steady state within R∼3, i.e. the inner factor tion at the disk midplane to the horizontal direction at
of ∼30 in radius, as evident in Figure 3. A steady state the disk surface, and then they reverse their directions
throughout the whole disk cannot be established with higher above. The net fields along the height at R = 1
our simulation setup since we do not supply new mate- areshowninFigure6. Inthiswork,wedefinethecorona
rial at the outer boundary. As shown in the α region (the yellow shaded region) as the region which is
total,mid
panel of Figure 3, the MRI saturates at later times at magnetically dominated (β < 1) and has a negative ra-
largerradii. Bytheendofoursimulation, theinnerdisk dialvelocity(v <1). Awayfromthemidplane,B first
r R
reaches a steady state and the disk mass accretion rate becomesnegativeandthenpositiveduetotheradialdrag
is almost a constant within R∼3, shown in the M˙ panel fromthesurfaceaccretion. Mostimportantly,theresult-
of Figure 3. inghorizontalfieldlinesarestretchedazimuthallydueto
Once a quasi-steady state is reached at t=42 T , the Keplerian shear, similar to the growth of azimuthal
0
the surface density profile follows ∼ 0.12R−0.6, and the fields in the linear phase of MRI. At z > 0, the nega-
midplane αrφ profile follows ∼ 0.1R0.5. The midplane tiveBR closetothemidplaneisstretchedtothepositive
stress then follows ∝R−1.85. The integrated α follows Bφ and the positive BR at the upper corona region is
rφ
stretched to the negative B (also shown in the color
∼ 0.5R−0.4, so that the vertically integrated stress fol- φ
mapofFigure5). SuchnetB iscrucialformaintaining
lows ∝ R−1.5 and the disk accretion rate is a constant φ
coronal accretion and launching disk wind. Net B is
with radii. With these values and Equation 5, we can z
alwayspositivebutitgetsamplifiedattheupperbound-
deriveM˙ =-0.0038,whichisconsistentwiththedirectM˙
ary of the corona region since the field lines are pinched
measurement in the upper right panel of Figure 3. We
there. Such net field geometry exerts strong stresses or
can see that the vertically integrated α is much larger
torques to the disk. Since B and B have opposite
thanthemidplaneαbyafactorof∼10. Thisisbecause φ R
signs as shown above, the Rφ magnetic stress, -B B ,
the stress is almost uniform vertically, or even increases R φ
is positive from the midplane to the wind region. The
towardsthedisksurface,whilethedensitydropssharply
z−φ stress -B B has similar shapes as -B since B is
towardsthedisksurface. Ifthestressisperfectlyuniform z φ φ z
always positive. At z > 0, the stress changes from posi-
vertically to some disk heights, the vertically integrated
tive values in the lower corona region to negative values
Rφstress should follow R−0.85 and the midplane α and
at the wind base. These magnetic stresses are shown in
the vertically integrated α should have the same slope.
Figure 7. In the corona region, both Rφ and zφ stresses
But the stress panel in Figure 3 suggests that the ver-
are mostly due to the mean field as shown by the blue
tically integrated Rφ stress is much steeper with R−1.5,
curves.
whichisbecausethestressislargerattheinnerdisk’sat-
Figures 6 and 7 show disk quantities along the disk
mosphere where the fields are pinched as shown in §5.1.
height (z) at R =1. The corona region is again marked
Thus, the midplane α and the vertically integrated α
as the yellow shaded area. We can see that the corona
have different slopes. Since α ∼ 0.5 at R=1, 42 T is
int 0 region has a supersonic inflow velocity, reaching ∼ 2 c
alreadylongerthantheviscoustimescale(R2Ω/(α c2)) s
int s and transports a significant amount of mass inwards.
atR=1,explainingwhythediskhasaconstantaccretion
The density reaches a plateau in this magnetically sup-
rate.
ported corona region. Figure 6 also suggests that the
While the disk is losing mass through accretion, it is
disk is sub-Keplerian in the corona region. Considering
alsolosingmagneticfield. Thenetverticalmagneticfield
the corona region is magnetically dominated, the sub-
decreaseswithtime,shownintheB panelofFigure
θ,mid Keplerian motion could be due to that magnetic fields
3. The plasma β calculated with the net vertical mag-
in the corona connect to the midplane at larger radii.
netic field increases from 103 to 105 at the inner disk
In other words, the midplane magnetically breaks the
(the lower center panel). On the other hand, the ratio
corona, as described in detail below. The <ρV > panel
betweenthenetverticalfieldandthedisksurfacedensity r
in figure 7 also reveals that, on average, the midplane
doesnotvarymuchovertime(theB/Σpanel). Although
transportsmassoutwardsfromt=35to42T . Thewind
itistemptingtointerpretthisasfluxfreezingduringthe 0
region which is beyond the corona region carries little
accretionprocess,wewillshowlaterthattheglobalfield
mass outwards since the density is low there (more dis-
transport is much more complicated.
6 Zhu & Stone
2 0.01
0
t=10
t=21
1
mid t=42 Σ 0.00
ρ10 0 t=0 g10 −1 ˙M
g o
Lo L −0.01
−1 R 0.6
−
−2 −2 −0.02
−1.0 −0.5 0.0 0.5 1.0 −1.0 −0.5 0.0 0.5 1.0 −1.0 −0.5 0.0 0.5 1.0
0.0 0.5
−2 R−1.85 mid −0.5 R0.5 0.0 R−0.4
gT10 −3 αtotal, −1.0 α10int −0.5
Lo 10 og
g L
−4 R 1.5 Lo −1.5 −1.0
−
−5 −2.0 −1.5
−1.0 −0.5 0.0 0.5 1.0 −1.0 −0.5 0.0 0.5 1.0 −1.0 −0.5 0.0 0.5 1.0
)
>
0.005 −2 0.00
P
< >
> 0.000 2/2 −3 <Σ
mid > >/ −0.02
<Bθ, −0.005 Bθ,mid −4 θ,mid
< −5 B
< −0.04
(
0
1
g
−0.010 o −6
−1.0 −0.5 0.0 0.5 1.0 L −1.0 −0.5 0.0 0.5 1.0 −1.0 −0.5 0.0 0.5 1.0
Log R Log R Log R
10 10 10
Fig. 3.— The disk midplane density, surface density, mass accretion rate (upper panels), stresses (the dashed curves are vertically
integratedRφstressesandthesolidcurvesarerφstressesatthemidplane),midplaneα,andverticallyintegratedα(middlepanels),and
netverticalmagneticfield,(cid:104)Bθ,mid(cid:105)2/2(cid:104)P(cid:105),(cid:104)Bθ,mid(cid:105)/(cid:104)Σ(cid:105)(lowerpanels)atdifferenttimes. αtotal andαint areallcalculatedwiththerφ
orRφstress. Inthepanelsatthebottomrow,quantitiesarealsoaveragedalongtheradialdirectionover1diskscaleheight.
cussion on mass loss due to the wind is in §4.2). We then write Equation 16 as
To understand how magnetic stresses can lead to the
coronal accretion, we analyze the angular momentum ∂(cid:104)ρδv (cid:105)
φ =m +m˙ +m +m˙ . (17)
budget due to these stresses in different regions. Fig- ∂t rφ r θφ θ
ure 4 suggests that the boundary between the corona
and wind region is around z∼1.5 R, and the boundary Thelefttermisthechangeoftheazimuthalmomentum.
betweenthemidplaneandcoronaregionisaroundz∼0.1 The first term on the right (m ) is the radial gradient
rφ
R.Thus,themidplane,corona,andwindregionsarebet- ofther−φstress. Atthediskmidplane, thistermisas-
terrepresentedbyθwedgesinthespherical-polarcoordi- sociatedwithdiskturbulence. However,withinthewind
nate instead of z slabs in cylindrical coordinates. Thus, region,thistermacceleratesthewindsincethemagnetic
we rewrite and analyze the angular momentum equation fields are aligned with the radial direction. The second
under the spherical-polar coordinate coordinate. After termontheright(m˙ )istheangularmomentumcarried
r
inserting the mass conservation equation and dividing by the radial accretion flow. It is the radial accretion
the equation by r, we get term. Whenitispositive,thediskaccretesinwards. The
third term (m ) is the θ gradient of the θ−φ stress. It
∂(cid:104)ρδv (cid:105) 1 ∂(r3(cid:104)T (cid:105)) (cid:104)ρv (cid:105)∂rv θφ
φ =− rφ − r k can also be considered as the torque between different
∂t r3 ∂r r ∂r layers in the disk. The forth term (m˙ ) is the angular
θ
1 ∂(sin2θ(cid:104)T (cid:105)) (cid:104)ρv (cid:105)∂(sinθv ) momentum loss due to the flow in the θ direction, such
− θφ − θ k as disk wind.
rsinθ2 ∂θ rsinθ ∂θ
(16) These different terms are integrated over θ wedges as
shown in Figure 8. The wedges are chosen in a way that
where they represent different disk regions at R = 1. At the
disk midplane (upper left panel), the m˙ term is nega-
T ≡ρv δv −B B r
rφ r φ r φ tive at R=1, suggesting that the flow is outwards. With
T ≡ρv δv −B B
θφ θ φ θ φ our setup, we can see that the m˙ term balances the
r
Global MHD Simulations 7
Fig. 4.— The azimuthally averaged density at different scales at t=42 T0. The green curves are the velocity field lines calculated with
azimuthally averaged velocities. The dashed curves label the Alfven surface while the dotted curves label the fast magnetosonic surface.
Clearly our simulation domain is beyond all these critical points. The white curves label where β = 1. The red curves label where the
azimuthally averaged density is only larger than the density floor by 10%, indicating the majority grids in the azimuthal direction have
reachedthedensityfloorthere. Onlythediskatmosphereattheveryinnerregionhasreachedthedensityfloor.
addition of two larger terms: m and m . At some the corona region. Since T is negative and positive at
rφ θφ θφ
other radii, the balance of two large terms can also lead the upper and bottom surfaces of the corona region, the
to positive m˙ (or inward accretion). It is also possi- integrated m term is negative. Thus, the θ−φ stress
r θφ
ble that with some other disk parameters, the midplane torques the disk surface to flow inwards and the disk
could have different flow directions. As shown in Fig- midplanetoflowoutwards,againlikemagneticbreaking.
ure 3, T follows r−1.85 at the disk midplane and it is At R=1, the coronal accretion rate is ∼10 times larger
rφ
almostaconstantwithheight(Figure9). Usingthedefi- than the midplane outflow rate.
nitionofm ,m is−1.15T /r whichisnegative. The For the disk as a whole, including both the midplane
rφ rφ rφ
negative m term means that the r−φ stress tries to and the corona, the large θφ stress within the disk can-
rφ
makethediskaccreteinwards. However,theθ−φstress not lead to the overall disk accretion. Only the θφ stress
term (m ) which is mainly due to net magnetic fields exerted at the upper and lower surface of the disk can
θφ
is positive, trying to make the disk flow outwards. The torque the disk, leading to accretion. The lower right
positive value of m can be understood from Figure 9, panel in Figure 8 shows that the total m term is neg-
θφ θφ
whereT becomeslargemovingawayfromthemidplane ative due to the θ−φ stress exerted by the wind to the
θφ
due to the azimuthal field stretch. Then, the m term disk surfaces. Thus, the wind torque contributes to the
θφ
(∼ −∂T /∂θ) is positive, driving disk outwards. We disk accretion. However, the θ−φ stress term is smaller
θφ
could also regard the m term as the magnetic break- than the r−φ stress term by a factor of 20 by examin-
θφ
ing. It torques the midplane to move outwards but the ing the lower right panel in more detail. Thus, 95% of
surface inwards. Eventually, at the disk midplane, the the disk accretion is due to the r−φ stress. Such stress
θ−φ stress term wins over the r−φ stress term at some is from MRI turbulence at the midplane and the global
radii and the disk flows outwards at R=1. magnetic fields in the corona (the T panel in Figure
R,φ
Inthecoronaregion,bothtermsduetor−φandθ−φ 7).
stressesarenegative,drivingthedisktoaccreteinwards. Tounderstandhowthecoronalaccretionisdetermined
We note that the m term that is vertically integrated by different components of stresses, we fit simple curves
θφ
with a weight of sinθ2 is basically the difference between for both Trφ and Tθφ stresses and use Equation 16 to
the θ−φ stresses at the upper and bottom surfaces of calculate the accretion rate at different heights. We fit
8 Zhu & Stone
Fig. 5.— Similar to Figure 4 but the color map is B and the green curves show the magnetic field lines. The color bars in the left,
φ
middle,andrightpanelarefrom-0.001to0.001,-0.01to0.01,and-0.1to0.1respectively.
T with the constant of 0.00032 and T with T stress and 5% from T stress exerted at the disk
rφ θφ rφ θφ
surface, which is consistent with the value by examining
(cid:40)0.00025×(θ−0.77) if θ <1.4
the components in the lower right panel of Figure 8.
Tθφ = −0.00096×(θ−π/2) if 1.4<θ <1.74 Sincebothturbulenceandnetmagneticfieldscancon-
0.00025×(θ−2.37) if 1.74<θ tribute to stresses, we would like to know their relative
importance. We separate the magnetic stress B B into
shown as the red curve in the middle panel of Figure 9. i j
((cid:104)B (cid:105) +δB )((cid:104)B (cid:105) +δB ) where (cid:104)B(cid:105) is the net field
Then with the assumptions that the disk is in a steady i φ i j φ j φ
state, m˙ term is negligible, and T ∝ r−1.85, we can that has been averaged over the φ direction. Thus the
θ rφ azimuthally averaged stress (cid:104)B B (cid:105) can be divided into
derive the accretion rate as i j φ
the stress due to net fields (cid:104)B (cid:105) (cid:104)B (cid:105) and the stress
i φ j φ
−0.00073−0.0005− due to turbulence (cid:104)δB δB (cid:105) . The stresses are shown
−00..0000017×3+ct0g.θ00×19(θ+−0.77) if θ <1.4 irn−thφesutrpepsesrispadnoemlsinoafteiFdigbujyrφeth1e0MwhRiIchturrebvuelaelnsctehaatt tthhee
ρv =
r 0.0038×ctgθ×(θ−π/2) if 1.4<θ <1.74 disk midplane. Previous work has established that in
−0.00073−0.0005− MRI-resolved models αRey is ∼1/4 of αMax, and the ra-
0.001×ctgθ×(θ−2.37) if 1.74<θ tio between the Maxwell stress and the magnetic pres-
sure α = T /P ∼ 0.45 (Hawley et al. 1995). In
shown as the red curve in the right panel of Figure mag Max b
our simulation, the MRI turbulent midplane also satis-
9. We can clearly see that the T stress torques the
θφ fiestheserelationships(Figure11). Figure10alsoshows
disk midplane to flow outwards and the disk surface to
that the r −φ stress is dominated by the net fields in
flow inwards. The integrated mass flux (with the weight
the corona region, and by the turbulence again at the
of sinθ2) at r = 1 is thus ∼ −0.00073(θ − θ )-
max min transition between the corona and wind regions. Since
2(T (θ )−T (θ )). Since the disk region extends
θφ max θφ min the r − φ stress determines the radial inflow, we con-
from θ =0.69 to θ =2.45, we can derive the inte-
min max clude that both turbulence and net fields contribute to
grated flux to be -0.00135.
the radial accretion. On the other hand, the net fields
If we only use the T stress, we can derive ρv as -
rφ r dominate the θ−φ stress in the wind and corona region
0.00073 and the vertically integrated disk mass flux as
andcontributessignificantlytothestressevenatthedisk
∼ -0.00073(θ −θ ). With θ and θ plugged
max min max min midplane. Since the θ−φ stress determines the internal
in, the mass flux due to T stress is -0.00128. From
rφ flow structure within the disk, we conclude that the net
this simple model, 95% of disk accretion is due to the
Global MHD Simulations 9
1.0 0.20 2.0
0.15
0.5 0.10 1.5
0.05
R z φ
V 0.0 V 0.00 V 1.0
−0.05
−0.5 −0.10 0.5
−0.15
−1.0 −0.20 0.0
−2 0 2 −3 −2 −1 0 1 2 3 −3 −2 −1 0 1 2 3
0.02
0.010 0.05
0.01
0.005
R z φ
B 0.00 B B
0.00
0.000
−0.01
−0.005
−0.05
−0.02 −0.010
−2 0 2 −3 −2 −1 0 1 2 3 −2 0 2
z z z
Fig. 6.— VariousvelocityandmagneticcomponentsalongthezdirectionatR=1. Thequantitieshavebeenaveragedbothazimuthally
andovertime(t=35to42T0 witha∆t=0.1T0 interval).
fields determine the vertically-sheared motion. 1982). This is satisfied in our simulation as shown in
Finally, we would like to know the total mass flux in Figure 4 and 5 where the dashed curves are the Alfven
differentregions. Thus,wecutthreewedgesinthesimu- surfacewhilethedottedcurvesarethefastmagnetosonic
lationdomainandmeasuretheradialmassfluxes,shown surface. Alfven surface is where the velocity in the
inFigure12. Wecanseethatthecoronalaccretiondom- poloidal plane V ≡ (cid:112)(cid:104)v (cid:105)2+(cid:104)v (cid:105)2 equals the poloidal
p r θ
inates. AtR=1,themidplaneoutflowrateis10%ofthe (cid:112)
Alfven velocity V ≡ ((cid:104)B (cid:105)2+(cid:104)B (cid:105)2)/(cid:104)ρ(cid:105). Fast
coronalaccretionrate. However,atsomeotherradii,the Ap r θ
magnetosonic surface is where V equals the fast mag-
disk flows inwards at the midplane. In any case, the net (cid:113) p
flow rate at the midplane is always much smaller than netosonic speed V ≡ 1/2(V 2+c 2)+1/2|V2−c2|
F A s A s
the inflow rate in the corona region. The wind region (cid:112)
where V ≡ ((cid:104)B (cid:105)2+(cid:104)B (cid:105)2+(cid:104)B (cid:105)2)/(cid:104)ρ(cid:105). In the
has a very small mass flow rate. At larger distances, our A r θ φ
disk’s atmosphere V is much larger than c , so that
cutforthewindregionincludespartofthediskregionso A s
V ∼V .
thatthemassflowratedoesnotreflectthetruewindloss F A
For a steady axisymmetric outflow, there are four
rate. A more detailed analysis on the outflow properties
conserved quantities along the magnetic and velocity
is presented in the next subsection.
field lines (Weber & Davis 1967; Blandford & Payne
4.2. Wind Region 1982). In a steady flow, the induction equation becomes
∇ × (v × B)=0. If v and B are separated into the
Although the torque by the wind at the disk surface
poloidal and the toroidal components (v = v +ΩRφˆ
only accounts for 5% of the disk accretion rate, the disk p
and B = B +B ), it has been shown that v and B
wind is launched at all the radii in the disk and it can p φ p p
are in the same direction, and using mass conservation
be directly probed by observations (Bjerkeli et al. 2016).
equation we have the first constant
Thus, we will study its properties in more detail.
To properly study the disk outflow/wind, the simula-
tion domain needs to be larger than the Alfven surface
ρv
and the fast magnetosonic surface (Blandford & Payne k = p , (18)
B
p
10 Zhu & Stone
−3
ρ
0
0 ρ
> −4
<ρ >/cs 0 otal
g10 −2 Vr Rφ,t
o < T
L −5
−2
−4
−6
−1 0 1 −1 0 1 −1 0 1
z z z
4 0.00010
0.001
0.00005
>
2
1<β− ρV>r 0.000 φ,total 0.00000
10 < Tz
g
o
0
L
−0.00005
−0.001
−2 −0.00010
−1 0 1 −1 0 1 −1 0 1
z z z
Fig. 7.—VariousquantitiesalongtheverticaldirectionatR=1. Thequantitiesareaveragedoverboththeazimuthaldirectionandtime.
The time average is taken from snapshots from t=35 to t=42 T0 with ∆t = 0.1T0 interval. The blue curves in the TRφ and Tzφ panels
aremagneticstressesthatarecalculatedusingthemeanfields,-BR×Bφ and-Bz×Bφ. Themeanfieldsarebothazimuthallyandtime
averagedbeforebeingusedtocalculatethestress. Theyellowshadedregionlabelsthecoronaregion. Notethefastinwardflowatthedisk
corona.
which is the mass load parameter, affecting the dynami- ways using the first three constants above
cal properties of the wind (Ouyed & Pudritz 1999). Us-
ing the same equations, the second constant can also be 1 B B B v
e= v2+Φ+h+ φ φ − φ φ (22)
derived 2 ρ k
1 RB ω
ω =Ω− kBφ . (19) = 2v2+Φ+h− kφ , (23)
ρR
Although our disk is locally isothermal which is not
Withtheangularmomentumequation,wehavethethird barotropic, wewillseethathismuchsmallerthanother
constant terms(whichmeansthewindis“cold”)andwecanstill
treat e as a constant.
l=R(v − Bφ) (20) Using the azimuthally averaged velocity structure, we
φ k derive velocity streamlines and the four conserved quan-
tities along the streamlines. We find that the conserved
which is the specific angular momentum of the wind. In quantities are not constants. For example, the values of
a barotropic fluid, we can use Bernoulli’s equation to k and ω have a factor of 2 peak at r ∼ 20 at t=42 T0.
derive the fourth constant Whenwecheckthesequantitiesovertime,wenoticethat
the peak in k originates from the disk surface and prop-
1 B2 B·v agatesoutwards. Thisisnotsurprisingsincetheoutflow
e= v2+Φ+h+ − (21) is launched from the turbulent disk. Previous simula-
2 ρ k
tions have also observed the episodic disk wind (Ouyed
& Pudritz 1997b).
(cid:82)
where h is dp/ρ. We can also write e in some other Thus, we average all primitive quantities over time,