Radial Transport and Meridional Circulation in Accretion Disks

0.67 MB·

by  Alexander A. Philippov

Draft version February 9, 2017
PreprinttypesetusingLATEXstyleemulateapjv.5/2/11
RADIAL TRANSPORT AND MERIDIONAL CIRCULATION IN ACCRETION DISKS
Alexander A. Philippov1,4 and Roman R. Rafikov2,3
Draft version February 9, 2017
ABSTRACT
Radialtransportofparticles,elementsandfluiddrivenbyinternalstressesinthree-dimensional(3D)
astrophysical accretion disks is an important phenomenon, potentially relevant for the outward dust
transportinprotoplanetarydisks,originoftherefractoryparticlesincomets,isotopicequilibrationin
theEarth-Moonsystem,etc. Togainbetterinsightintotheseprocesses,weexplorethedependenceof
7
meridionalcirculationin3Ddiskswithshearviscosityontheirthermalstratification,anddemonstrate
1
strong effect of the latter on the radial flow. Previous locally isothermal studies have normally found
0
a pattern of the radial outflow near the midplane, switching to inflow higher up. Here we show, both
2
analytically and numerically, that a flow, which is inward at all altitudes, is possible in disks with

b entropy and temperature steeply increasing with height. Such thermodynamic conditions may be
e typical in the optically thin, viscously heated accretion disks. Disks in which these conditions do not
F
hold should feature radial outflow near the midplane, as long as their internal stress is provided by

7 the shear viscosity. Our results can also be used for designing hydrodynamical disk simulations with
a prescribed pattern of the meridional circulation.

] Subject headings: accretion, accretion disks — protoplanetary disks — hydrodynamics
E
H
. 1. INTRODUCTION inward than is transported out by the midplane outflow,
h
resulting in a net (vertically integrated) inflow of mass
p Evolution of astrophysical accretion disks is believed
towardstheaccretingobject(Urpin1984;Jacquet2013).
o- tpoalobiezoduri&venLipnri1m9a9r5i)ly. bGyrotshseifreaitnuterernsaolfstthreisssepsro(cPeass- Thus,eventhough(cid:104)uR(cid:105)ρ =Σ−1(cid:82) ρuRdzisstillnegative,
r can be understood by treating the disk as a geometri- the radial outflow near the midplane actually transports
t
s cally thin structure and characterizing its properties us- mass out in 3D. This remarkable analytical prediction
a ingtheverticallyintegrated(orproperlyaveragedacross was later confirmed by numerical hydrodynamical cal-
[ thediskthickness)variablessuchasthesurfacedensityΣ culations of Siemiginowska (1988), Kley & Lin (1992),
2 (Shakura&Sunyaev1973;Lynden-Bell&Pringle1974). andRozyczkaetal.(1994),whichexplicitlyemployedα-
v Inthisapproachthesteadyaccretiondiskswiththeradi- prescription (Shakura & Sunyaev 1973) to describe the
2 ally constant mass accretion rate M˙ necessarily exhibit internal stress.
1 purely inward radial motion of gas, driven by the an- The meridional circulation, if indeed present in accre-
9 tion disks, should have very important implications. It
gular momentum redistribution in the disk due to its
1 may provide a natural way of coupling the inner disk to
internal stresses. In other words, the density-weighted
0 its outer parts: the material and information from the
vertical average of the radial velocity u is always nega-
1. tive, (cid:104)u (cid:105) =Σ−1(cid:82) ρu dz <0. R vicinity of an accreting object could potentially be di-
R ρ R rectlytransportedintheadvectivefashionfartheroutin
0 Situation becomes more complicated once one consid-
the disk by the near-midplane outflow.
7 ers the full three-dimensional disk structure, in partic-
Toprovidesomeillustrativeexamples,Takeuchi&Lin
1 ular, the profile of the radial velocity u as a func-
: tion of the vertical coordinate z. It was fiRrst shown by (2002) have suggested that meridional circulation can
v drive the outward transport of dust grains in protoplan-
Urpin (1984) that when the angular momentum trans-
Xi port in the disk is effected by the effective shear viscos- etary disks. More recently, the samples collected by
the Stardust mission from comet 81P/Wild 2 revealed
r ity (Shakura & Sunyaev 1973), the radial velocity of the
a gas is actually positive at z = 0, implying a radial ou- the presence of large (> 1µm) grains composed of high-
temperature minerals that appear to have formed in the
flow near the disk midplane. The magnitude of positive
inner regions of the Solar nebula (Brownlee et al. 2006;
u steadilydecreaseswiththeverticaldistancefromthe
R Zolensky et al. 2006). Their presence in comets that
midplane, and u ultimately changes sign at some al-
R form at tens of AU indicates that global outward radial
titude of order the disk scale height H. As a result, a
transportofdustparticlescouldhaveoperatedintheSo-
meridional circulation pattern sets in the poloidal plane
lar nebula (Ciesla 2007, 2009; Hughes & Armitage 2010;
ofthedisk. Radialgasinflowathighz carriesmoremass
Jacquet 2013).
MoonisthoughttohaveformedinacollisionofaMars-
1DepartmentofAstrophysicalSciences,PrincetonUniversity, size planetary embryo with the proto-Earth (Hartmann
IvyLane,Princeton,NJ08540 &Davis1975),primarilyoutofthematerialderivedfrom
2InstituteforAdvancedStudy,EinsteinDrive,Princeton,NJ
the impactor (Canup 2004). This should have resulted
08540
3Centre for Mathematical Sciences, Department of Applied inthedifferentcompositionsoftheMoonandtheEarth.
MathematicsandTheoreticalPhysics,UniversityofCambridge, However,theisotopicratiosofoxygenderivedfromlunar
WilberforceRoad,CambridgeCB30WA,UK
[email protected]
2
Apollo samples were found to be essentially identical5 to ally hold in accretion disks dominated by viscous dissi-
their terrestrial analogs (Wiechert et al. 2001). To ex- pation. Some earlier numerical studies (Siemiginowska
plainthispuzzlePahlevan&Stevenson(2007)suggested 1988; Kley & Lin 1992; Rozyczka et al. 1994) have at-
that equilibration of oxygen isotopes in the Earth-Moon tempted to use more refined treatments of the disk ther-
system resulted from the rapid radial mixing between modynamics. However, they were limited in resolution
theterrestrialmagmaoceanandtheproto-Lunardiskof and did not explore the full range of possible thermody-
vapor-melt debris produced by the impact. Pahlevan & namic regimes. Our present goal is to provide a thor-
Stevenson (2007) appealed to diffusive turbulent mixing ough analysis of the effects of thermal stratification on
as the mechanism of compositional equilibration inside the meridional circulation in the disk.
the disk. However, meridional circulation can also be This work is organized as follows. In §2 we present
an important contributor to this process by directly ad- general theoretical description of meridional circulation
vecting into the disk and mixing the material from the in accretion disks. We explore conditions, under which
magma ocean via the midplane outflow. the inflow occurs at the disk midplane in §3. In §4-5
Weakly magnetized accreting objects (cataclysmic we verify our theoretical results with 3D viscous numer-
variables,classicalnovaprogenitors,neutronstarsinlow- ical simulations. We discuss our findings in §6, focusing
mass X-ray binaries, etc.) are expected to develop a on applications to real systems in §6.3. Our results are
boundary layer between the inner edge of the accretion summarized in §7.
diskandthestellarsurface(Pophametal.1993;Popham
& Narayan 1995; Belyaev et al. 2012, 2013a,b). If the 2. THEORETICALCONSIDERATIONS
material from the near-surface layers of the accretor can We consider a geometrically thin, axisymmetric disk
be dredged up (evidence for this exists in e.g. classi- orbiting in a central potential of a point mass M ; we
(cid:63)
cal novae Truran & Livio 1986; Gehrz et al. 1998) into work in cylindrical coordinates (R,z). We are interested
the boundary layer by some internal processes, then the in the spatial structure of the radial velocity u (R,z),
R
near-midplane outflow could subsequently transport it which is non zero because of the internal stress T op-
ij
to the outer parts of the disk. Such elemental pollution erating in the disk. We are interested only in advective,
mightspreadthemetalswithunusualabundancesacross laminar motions of the fluid and do not consider turbu-
the disk, affecting its observational appearance (e.g. via lentdiffusion,whichisofteninvokedinstudiesofparticle
spectral signatures). transport in disks (Takeuchi & Lin 2002; Hughes & Ar-
Given the variety of situations, in which meridional mitage 2010). Thermal structure of the disk is specified
circulation may be important, it is natural to ask how in §2.1.
robust this phenomenon is. The majority of its ana- Equationsoftheradialandverticalbalance,describing
lytical and numerical investigations focused on disks, in the equilibrium disk structure, are
which angular momentum transport is accomplished by
1∂P ∂Φ
the conventional α-viscosity (Shakura & Sunyaev 1973). Ω2R− = , (1)
All of them invariably find u (R,0) > 0 implying the ρ∂R ∂R
R
near-midplane outflow. However, in many types of disks 1∂P ∂Φ
the transport is more likely to be effected by the mag- =− , (2)
ρ ∂z ∂z
netorotational instability (MRI, Balbus 2003), for which
the stress tensor T is anisotropic. The anisotropy of whereρandP arethegasdensityandpressure, Ωisthe
ij
Tij has been shown (Jacquet 2013) to play an important angular frequency, and Φ(R,z) = −GM(cid:63)(R2 +z2)−1/2
role for the meridional circulation, although the numer- isthecentralpotential(weneglectthediskself-gravity).
ical evidence on the issue is rather mixed at this point Intheseequationswealsoneglectedviscoustermsresult-
(Fromang et al. 2011; Flock et al. 2011; Suzuki & Inut- ing from the fluid motions in the meridional plane; we
suka 2014), see the discussion in §6. But even without comment on the significance of this simplification later.
magneticfields,globalsimulationsbyStoll&Kley(2014, Because of the variation of the central gravity with
2016)ofthepurelyhydrodynamicvertical shear instabil- height and radial pressure support in the disk the an-
ity(VSI;Urpin&Brandenburg1998;Urpin2003)exhibit gular frequency deviates from its Keplerian value Ω ≡
K
radial inflow at the midplane, changing to an outflow at (GM /R3)1/2:
(cid:63)
high z, contrary to the expectations for disks mediated
by the shear stress (Urpin 1984; Takeuchi & Lin 2002). (cid:18) 3 z2 1 P ∂lnP(cid:19)
Ω(R,z)≈Ω 1− + . (3)
Inthisstudyweaddressadifferentaspectoftheprob- K 4R2 2Ω2R2 ρ ∂lnR
lem and explore the effect of thermal properties of the K
disk on the meridional circulation, while still confining This approximate relation follows from our assumption
ourselvestotheshearviscosityasthesourceoftheangu- ofageometricallythindisk,forwhichH/R(cid:28)1, z/R(cid:28)
larmomentumtransport(whichallowsustoperformde- 1, where H(R,z) ≡ cs/ΩK is the local disk scaleheight,
tailed analytical investigation). Existing analytical and and c2(R,z) ≡ P/ρ = k T(R,z)/µ is the isothermal
s B
numerical work has typically assumed the disk to have sound speed. Note that the last term in parentheses is a
locally isothermal vertical structure (Urpin 1984; Fro- function of both6 R and z. The deviations from purely
mangetal.2011;Jacquet2013;Suzuki&Inutsuka2014). Keplerianrotationinequation(3)areatthe(z/R)2 (cid:28)1
This may be a reasonable assumption for the externally (second term in parentheses) and (H/R)2 (cid:28) 1 (third
irradiated protoplanetary disks, but should not gener- term) level.
5 Terrestrial O isotopic ratios are quite distinct from e.g. their 6 In particular, one can show with the aid of equation (2) that
martianvalues. Ωisindependentofz forabarotropicequationofstate.
3
Theφ-componentoftheequationofmotiondescribing (Takeuchi & Lin 2002; Fromang et al. 2011). The two
the angular momentum balance can be reduced (under coincideifρ (R)∝Rp,butthepowerlawdensityprofile
0
reasonable assumptions regarding the amplitudes of the isnotnecessaryforourresultstobevalid. Notethatthe
radial and vertical fluid velocities, Takeuchi & Lin 2002; expression(6)doesnotmakeanyassumptionsaboutthe
Fromang et al. 2011) to the following expression for the verticalorradialdependenceoftheviscosity,althoughit
radial velocity u of the disk fluid (Jacquet 2013): does rely on the stress model (5).
R
Thetwotermsinthefirstlineoftheequation(6)orig-
1(cid:18) ∂l (cid:19)−1 inate from T , while the last two terms are due to the
u (R,z)=− Rφ
R ρ ∂R vertical stress Tzφ. Since the midplane pressure P0 is a
(cid:20) (cid:21) decreasing function of R, it is clear that a vertical den-
× 1 ∂ (cid:0)R2T (cid:1)+ ∂ (RT ) . (4) sity profile steeply declining with z (i.e. high value of
R∂R Rφ ∂z zφ −∂2ρ/∂z2) should be favorable for suppressing the mid-
Here l ≡ ΩR2 is the specific angular momentum of the planeoutflowinthediskanddrivinggasinflow(uR <0)
at all heights.
disk fluid, while T and T represent the R−φ and
Rφ zφ Note that equations (1)-(2) neglect contributions from
z−φ components of the internal stress tensor. This ex-
pression is accurate at the (z/R)2,(H/R)2 (cid:28) 1 level; it viscous stresses arising from shear in the meridional
plane. Such terms were considered in studies of Kley
neglects certain terms proportional to the vertical veloc-
& Lin (1992), Kluzniak & Kita (2000), and Regev &
ityandpossibletimevariabilityoftheazimuthalvelocity
(Fromang et al. 2011). It demonstrates that the radial Gitelman (2002), resulting in the additional O(α2) rel-
fluid motion of the 3D disk is determined not only by ative contribution to the expression (3) for the angular
the horizontal but also by the vertical component of the frequency. Their inclusion would, in turn, lead to the
stress tensor. This turns out to be of crucial importance emergenceofanadditionalO(α3)contributionintheex-
for our work. On the other hand, multiplying equation pression (6) for uR. Given the expectation of α (cid:28) 1 in
(4) by ρ and integrating it over z, one finds that (cid:104)u (cid:105) real disks, the omission of the viscous terms should not
is determined by the T behavior only (Jacquet 201R3)ρ. present a problem7. This expectation is subsequently
Rφ
In real astrophysical disks internal stress can be pro- confirmed in §5 by the good agreement between our an-
videdbyavarietyofmechanisms,withtheMRI(Balbus alytical theory (which explicitly neglects viscous contri-
2003) being one of the important possibilities. However, butions to Ω(R,z)) and numerical results.
our present focus is still on standard shear viscosity pre-
2.1. Disk thermodynamics
scription, for which one has
So far our treatment was fully general. To make fur-
∂Ω ∂Ω therprogressweneedtolookinmoredetailintothedisk
T =−ρνR , T =−ρνR , (5)
Rφ ∂R zφ ∂z thermodynamics. We will assume a rather general equa-
tion of state (EOS) for the disk fluid in the polytropic
where ν is a kinematic viscosity.
form
Previous studies (Fromang et al. 2011; Jacquet 2013)
clearly demonstrated that the possible meridional out- P =esργ, s(R,Z)=c−1S(R,z), (8)
V
flow is most pronounced at the midplane of the disk, at
z = 0. For that reason, our primary goal will be to de- where γ is the polytropic index and s(R,z) is the scaled
termine the factors that affect the midplane value of the (dimensionless)gasentropyS(R,z);cV ≡(γ−1)−1kB/µ
radial velocity, u (R,0). To that effect, in our analysis is the specific heat capacity. Note that in this work
R
we will often take the limit z → 0 in evaluating differ- S(R,z) is specified explicitly, i.e. we do not attempt to
entexpressionsandaccountforthesymmetryofthedisk calculate it, e.g. from the energy equation by solving for
with respect to its midplane. The latter implies, in par- the radiation transfer in the disk.
ticular that the first derivatives of various fluid variables An often used assumption of a locally isothermal
with respect to z, e.g. ∂ρ/∂z, ∂Ω/∂z, etc., go to zero as disk structure (Takeuchi & Lin 2002; Fromang et al.
z →0. 2011; Jacquet 2013) implies that P/P0 = ρ/ρ0 =
In Appendix A we show that after some straightfor- exp(−z2/2H2), so that
0
ward manipulations the midplane value of the radial ve-
γ−1(cid:20) z (cid:21)2
locity can be written as siso(R,z)=s (R)+ , (9)
ν (cid:34)15 ∂ln(cid:0)ρ ν R3(cid:1) 0 2 H0(R)
uR(R,0)=R 2 −3 ∂0R0 where s0 = ln(P0/ργ0). A different limit of the locally
isentropic disks considered later in §3.2.1 has s(R,z) =
(cid:12) (cid:35) s(R).
P ∂lnP ∂2ρ(cid:12)
−δ − 0 0 × (cid:12) , (6) Using the ansatz (8) we show in Appendix A that the
ρ Ω2Kρ20 ∂lnR ∂z2(cid:12)(cid:12) lastterminsidethebracketsintheexpression(6)canbe
z→0
written as
where
(cid:32) (cid:12) (cid:33)
δ +δ P ∂2s(cid:12)
δ ≡∂lnρ /∂lnR, (7) ρ T 1+ (cid:12) , (10)
ρ 0 γ Ω2ρ∂z2(cid:12)
K (cid:12)
and subscript 0 implies the value of a particular variable z→0
7 Kluzniak & Kita (2000) and Regev & Gitelman (2002) find
at the midplane, e.g. ρ (R)≡ρ(R,0), P (R)≡P(R,0),
0 0 thatinclusionoftheviscoustermschangesmidplaneoutflowtoan
etc. Definition(7)impliesthatδρ iscloselyrelatedtothe inflowonlyforα(cid:38)0.7,althoughKley&Lin(1992)suggestalower
density power law exponent p used in previous studies valueofcriticalα≈0.06forthistransition.
4
where where we introduce the shorthand notation
(cid:12)
δ ≡∂lnT /∂lnR, (11) ∂2s(cid:12)
T 0 s(cid:48)(cid:48) ≡H2 (cid:12) (17)
zz 0∂z2(cid:12)
and T (R) ≡ T(R,0) is the midplane value of the gas (cid:12)
0 z→0
temperature. Definition of δ makes it closely related to
T for the (dimensionless) second derivative of the entropy
the temperature power law index q used in other studies
at the disk midplane, and
ofthemeridionalcirculationinaccretiondisks(Takeuchi
&Lin2002;Fromangetal.2011),whichassumeT0(R)∝ δF ≡∂lnFJ/∂lnR. (18)
Rq.
Nextweconnectthebehaviorofthesecondterminthe Once again, in equation (16) the first two terms arise
equation(6)totheglobalstructureofthedisk. Weintro- from TRφ, while the terms in the second line originate
duce the viscous angular momentum flux across a given from Tzφ.
radius F , which is equal to the viscous torque exerted Equation (16) represents the main analytical result of
J
bytheinnerdiskontheouterdisk. TheconceptofF is this work. It provides a connection between the ampli-
J
known to be very useful for describing viscous accretion tude (and direction) of the radial velocity at the disk
disks in a variety of situations, especially a steady state midplane and the vertical thermal stratification charac-
(Lynden-Bell & Pringle 1974; Rafikov 2013, 2016a). Ac- terized by s(cid:48)(cid:48) . This link has not been established in
zz
cording to the definition, F can be calculated by multi- previous studies of meridional circulation.
J
plyingT inequation(5)by2πR2 andintegratingover
Rφ
z. Using equation (5) this implies, to (H/R)2 accuracy, 3. CONDITIONSFORINFLOWATALLZ
that Equation (16) allows us to determine the conditions
underwhichthediskwillexhibitinflowatallz,including
∞ ∞
FJ =− (cid:90) 2πR3ρν∂∂RΩdz =3πl (cid:90) ρνdz, (12) tchrietemriiodnplfaonret.hSisetttoinbgeutRh(eRc,a0s)e:<0onefindsanecessary
γ
−∞ −∞ s(cid:48)(cid:48) >−
zz δ +δ
wherel≡Ω R2 andwehaveassumedaKeplerianrota- ρ T
K (cid:18) (cid:19)
tionprofile. Also, multiplyingequation(4)by2πRρand 9 3γ+2 1−γ
× −3δ + δ + δ , (19)
integrating over z one finds that the mass accretion rate 2 F 2γ T γ ρ
M˙ (defined to be positive for inflow) is related to F as
J
where we assumed that P (R) decreases with R so that
(Rafikov 2013) 0
δ +δ <0. When this inequality is fulfilled meridional
ρ T
∞ circulation is unable to convey information and mate-
(cid:90) ∂F
M˙ (R)=−2πR ρu dz = J, (13) rial from the inner disk to the outer disk — only inward
R ∂l propagation is allowed.
−∞ Equation (19) is the most general form of the inflow
criterion that does not make any assumptions about the
where the derivative is with respect to the specific angu-
lar momentum, so that ∂/∂l=2(Ω R)−1∂/∂R. valuesofδF,δρ,andδT —theycantakearbitraryvalues
K
allowing one to explore meridional circulation even in
We now make an additional assumption that the disk
evolving disks. Moreover, even though it does assume
structure obeys a certain similarity property, namely
a particular stress model given by the equation (5), it
that
makesonlyaweakassumptionabouttheactualbehavior
(cid:18) (cid:19)
z oftheviscosity(thatν dependsonz onlyincombination
S(R,z)=S R, , H (R)≡H(R,0),(14)
H0(R) 0 z/H0).
where H0 = cs(R,0)/ΩK is the local value of the scale- 3.1. Inflow criterion for a standard constant M˙ = disk
height at the disk midplane. This is not a very con-
In practice, one is often interested in accretion disks
straining assumption and its adoption should not limit
that have reached a steady state. One of the most pop-
the applicability of our subsequent results.
ular assumptions used in many studies is that of the
With such entropy behavior it follows from equation
radially constant mass accretion rate through the disk
(8) that ρ and P also depend on z only in combina-
tion z/H ; this makes it natural to expect that ν = (M˙ =const), with no torque applied at its center. In
0
ν(R,z/H ) as well. Then equation (12) implies that this case, integrating equation (13) one finds F = M˙ l
0 J
(Rafikov 2013), so that δ = 1/2. This transforms in-
F
F ∝ρ ν H l∝ρ ν T1/2R2, (15) equality (19) into the following inflow criterion:
J 0 0 0 0 0 0
(cid:18) (cid:19)
allowing us to tackle the second term in equation (6). γ 3γ+2 1−γ
s(cid:48)(cid:48) >− 3+ δ + δ . (20)
Plugging equations (10) and (15) into the expression zz δ +δ 2γ T γ ρ
ρ T
(6) we finally find that
NotethattheconstantM˙ assumptiondoesnotconstrain
(cid:20) (cid:18) (cid:19)
ν 9 δ
u (R,0)= −3 δ − T δρ or δT because of the freedom in choosing the radial
R R 2 F 2 variationoftheviscosityν . Thisconstraintisillustrated
0
(cid:21) in Figure 1 for two values of γ (5/3 and 7/5) and three
δ +δ
−δ + ρ T (1+s(cid:48)(cid:48) ) , (16) values of δ . One can see that steeper decay of ρ with
ρ γ zz ρ 0
5
Thisresults inyetanotherversionofthe inflowcriterion
for α disks, which depends only on the radial tempera-
ture profile (i.e. δ ):
T
6(2γ−1)+(6γ−1)δ
s(cid:48)(cid:48) > T. (22)
zz 6+δ
T
This constraint is also shown in Figure 1.
3.2. Effects of thermal structure of the disk
Next we assess how the different assumptions about
the vertical thermal structure of the constant M˙ disk
affect the possibility of the inflow at its midplane.
3.2.1. Locally isentropic disk
The simplest thermodynamic assumption is that of a
locally isentropicdisk, inwhichentropy(ands)doesnot
depend on z. Such disks naturally have temperature de-
creasing with height. The isentropic vertical stratifica-
tion may arise e.g. if the disk is convectively unstable.
Inthiscasethelefthandsideofequation(20)becomes
zero and the condition of inflow turns into a constraint
on the radial behavior of ρ (R) and T (R) in the form
0 0
6γ+(3γ+2)δ
δ > T. (23)
ρ 2(γ−1)
Forprotoplanetarydiskswithγ =7/2thisbecomesδ >
ρ
(42+31δ )/4. In hotter accretion disks with γ = 5/3
T
the constraint is δ >3(10+7δ )/4.
ρ T
It is clear that even if T (R) decays with radius as
0
rapidly as R−1 (i.e. δ = −1) the inflow at the disk
T
midplanewouldrequiremidplanedensityincreasing out-
wards. This is a pretty unusual arrangement, which
makes midplane inflow essentially impossible in locally
isentropic disks, including the disks which are convec-
tively unstable.
3.2.2. Locally isothermal disk
A popular assumption of the locally isothermal disk
Fig. 1.—Minimumvalueofs(cid:48)z(cid:48)z —thesecondderivativeinz/H0 structure (∂T(R,z)/∂z = 0) turns the left hand side of
of the scaled entropy s = S/cV — necessary for having radial the inflow criterion into γ−1, see equation (9). Then it
inflowatthediskmidplane(i.e. uR(R,0)<0occursaboveacor-
respondingcurve)inaconstantM˙ disk(δF =1/2),asafunction is easy to see from the inequality (20) that the inflow at
ofthemidplanetemperatureslopeδT (Eq. [11]). Twovaluesofγ all heights is possible only if δT < −6/5. This result is
are considered: (a) 7/5 and (b) 5/3. Black curve is a constraint completely independent of either δ or γ, thus it should
ρ
(22)thatappliesinadiskwithradiallyconstantα. Othercurves hold for an arbitrary radial profile of the midplane den-
correspondtoequation(20),whichdoesnotassumeα=const,for
sity ρ (R).
differentvaluesofthemidplanedensityslopeδρ (Eq. [7]). Dotted 0
linesarerelevantforthelocallyisothermaldisks(§3.2.2): horizon- Finding conditions in which a disk might have T0(R)
tallinemarksthevalueof∂2s/∂(z/H0)2=γ−1insuchadisk,see decayingwithradiusfasterthanR−6/5isnoteasy. Based
Eq. (9). The vertical line marks δT =−6/5 — in locally isother- on this we can conclude that the locally isothermal con-
δmTalbdelioskwstuhRis(Rva,l0u)e.< 0, i.e. midplane inflow, is possible only for stant M˙ disks are predetermined to naturally exhibit an
outflow at the midplane, in agreement with a number of
R (more negative δ ) makes it easier to achieve inflow previous studies (Urpin 1984; Jacquet 2013).
ρ
at the disk midplane, i.e. requires less extreme values of Summarizing the results of this and previous subsec-
s(cid:48)(cid:48) . tion (§3.2.1), we conclude that the inflow at z = 0
zz
Wecannowgoonestepfurtherandadoptaparticular requires the disk to have temperature increasing with
viscosity ansatz, namely the α-model of Shakura & Sun- height, starting from the midplane. Such disks natu-
yaev (1973), in which ν = αc2/Ω . If we additionally rally have entropy rising more steeply with z than in
s K
assume that the effective viscosity parameter α is inde- the isothermal case and may satisfy the conditions (20)
pendent of R, then locally ν0 ∝ RδT+3/2. Plugging this and (22) without making unrealistic assumptions about
and F =M˙ l into equation (15) one obtains that the behavior of ρ0(R) and T0(R).
J
3 3.3. Inflow criterion for F = const disk
δ =−3− δ . (21) J
ρ 2 T
6
Syunyaev & Shakura (1977) in the context of accretion
by the magnetized neutron stars in the propeller regime
(Illarionov & Sunyaev 1975). The global structure of
such disks is characterized by F (R) =const, implying
J
that δ =0.
F
Equation(16)demonstratesthatreachinginflowatthe
midplaneofsuchadiskismore difficultthaninthestan-
dard constant-M˙ disk. Indeed, the inflow criterion (20)
remains essentially the same, however the first term in-
side the parentheses changes from 3 to 9/2. As a result,
a higher value of s(cid:48)(cid:48) is needed to guarantee inflow at all
zz
heights (i.e. u (R,0) < 0) in the constant-F disk, as
R J
compared to the standard constant M˙ disk (§3.1). In-
equality (22) gets modified in a similar fashion.
Moreover,repeatingthecalculationsof§3.2.2onefinds
thattheinflowatz =0inalocallyisothermal,constant-
F diskwouldrequireδ <−9/5. Suchsteeplydeclining
J T
radial profiles of temperature are unlikely in real disks.
4. NUMERICALSETUP
In order to confirm the analytical result (16), we per-
formed viscous hydrodynamical simulations using new
Godunov code Athena++ code (Stone et. al. 2016,
in preparation). Compared with its predecessor Athena
(Gardiner&Stone2008),Athena++ishighlyoptimized
and uses flexible grid structures, significantly facilitat-
ing global numerical simulations. In this work we per-
form our simulations in spherical r,θ,φ coordinates, us-
ing uniform grid in logr and constraining ourselves to
small range in φ. We use the domain [1,10]×[π/2−
0.5,π/2+0.5]×[0,0.1] in these coordinates with numer-
ical resolution of 1024×512×16 cells. We verified that
simulations in the full φ range [0,2π] lead to same re-
sults. We also verified that our results are converged
with regard to numerical resolution.
For r and θ boundaries, we use the so-called ”do-
nothing” boundary condition, i.e. all fluid variables are
fixed at their initial values. In φ-direction, we impose a
periodic boundary condition. The Mach number of the
orbital flow at the inner radius is 10 in all our simula-
tions. The simulation time is always measured in units
of Ω−1, where Ω is the Keplerian angular velocity at
K K
the inner boundary.
In order for the disk structure to remain unchanged
Fig. 2.—Verticalprofilesofthedensity(a)andtemperature(b)
and maintain its initially prescribed entropy profile dur-
inaconstantM˙ diskwiththeradialtemperatureprofilesuchthat
ingthesimulation,weusetheopticallythincoolingfunc-
δT = −1 (see Eq. [11]) and two vertical entropy profiles charac-
terized by s(cid:48)(cid:48) =0.5 (blue) or 1.5 (green), see the definition (17). tion:
zz
Ainlgl ovnartiahbelevsalaureeonfosr(cid:48)m(cid:48) ,aloizuerdthbeyrmthoedirynmamidipclaannezavtazlu(e8s).caDnepreesnudlt- Λ=−ρT −T0, (24)
in either vertically dezczaying (for low s(cid:48)(cid:48) ) or rising (for high s(cid:48)(cid:48) ) τ
zz zz
temperatureprofilesnearthemidplane. where T is the gas temperature, T (R,z) is the initial
0
Standard constant-M˙ disk with zero central torque temperature profile described in §4.1, τ = 0.5Ω−K1 is the
cooling time. This cooling prescription ensures that vis-
represents just one particular example of a steady disk.
cous heating does not significantly influence our initial
Amoregeneraltime-independentstructureofthediskis
describedbyF (R)=M˙ l+F (0)(Rafikov2013,2016a), diskconfiguration,andthatoursimulationsreachsteady
J J state.
where F (0) represents the torque applied at the disk
J
center(whichisnaturallyabsentinastandardconstant-
4.1. Initial entropy and density profiles
M˙ case).
For simplicity, we assume the power-law behavior of
A particularly interesting case to explore is that of
the density and temperature at the mid-plane ρ (R) ∝
a disk with no mass accretion at the center, in which 0
the gas inflow is fully suppressed by a strong cen- Rδρ,T0(R)∝RδT. Oursimulationsadoptasetuptypical
tral torque (Lynden-Bell & Pringle 1974; Pringle 1991; for standard constant M˙ disks (i.e. we do not simulate
Rafikov 2016b). Such ”dead” disks were first studied by disks considered in §3.3), making different assumptions
7
Fig. 3.— Maps of radial velocity uR(R,z), averaged over azimuth and normalized by ν/R, for T(R) profile with δT = −1, constant
accretion rate M˙ (i.e. δF = 1/2) and two values of the second derivative of entropy with height, s(cid:48)z(cid:48)z: (a) s(cid:48)z(cid:48)z = 0.5 (entropy slowly
increasing with z), which exhibits outflow at the midplane and inflow at high altitudes; (b) s(cid:48)(cid:48) =1.5 (entropy rapidly growing with z),
zz
whichshowsgasinflowatallaltitudes.
Fig. 4.—SlicesofradialvelocityuR(R,z), averagedoverazimuthandnormalizedbyν/R forδT =−1andconstantaccretionrateM˙.
(a): Slices of normalized uR at different radii from the center at t = 100 for s(cid:48)z(cid:48)z = 0.75, demonstrating the self-similar character of the
meridionalcirculation(allthreecurvesfallontopofeachother);(b)NormalizeduR profilesatR=5fors(cid:48)z(cid:48)z =0.75atdifferentmoments
of time, showing the convergence of uR to a steady-state; (c): Same at R = 5 in steady state for different values of s(cid:48)z(cid:48)z, characterizing
the falloff of the vertical entropy profile. Theoretical predictions for uR at z =0 (equation 16) are shown as stars. One can see both the
excellentagreementwiththeoryandthefactthatsteeplyincreasingentropyprofilesresultingasinflowatallaltitudes(see§5).
about the radial profile of α. Close to the mid-plane, we urallydevelopstemperatureprofilerisingwithheight. As
consider the following entropy profile we show below, in this case one finds gas inflow at the
diskmidplane. Onthecontrary,lowvaluesofs(cid:48)(cid:48) <γ−1
s(cid:48)(cid:48) (cid:20) z (cid:21)2 result in T dropping with z (by design, at highzzaltitudes
s(R,z)=s (R)+ zz , (25)
0 2 H (R) our T profiles always converge to isothermal).
0
To check our analytical prediction (16), we run two
such that the definition (17) holds true; also s (R) = sets of simulations. In a first set, we adopt a radially
0
ln(P /ργ). In order to prevent large gradients of fluid constant α profile with α = 0.01, so that equations (21)
0 0
variables from developing high above the mid-plane, we and (22) hold. In a second set, we assume that α ∝ R,
modifythisentropybehaviorinsuchawaythatthedisk with α = 0.01 at R = 1. In this case, the constant
becomes locally isothermal at |z|>2H0. M˙ assumption leads to a different relation between the
To compute the vertical disk structure corresponding power-law indices of density and temperature profiles:
to this entropy behavior, we numerically integrate the
equation of hydrostatic equilibrium (2). The azimuthal
3
component of the gas velocity is then computed using δ =−4− δ . (26)
equation (1). We show disk profiles for δ = −1, con- ρ 2 T
T
stantM˙ andtwovaluesofs(cid:48)(cid:48) inFigure2. Itshowsthat
zz
forhighvaluesoftheentropyderivatives(cid:48)(cid:48) thedisknat- Plugging this δ in equation (20), we obtain a new crite-
zz ρ
8
forT(R)profilewithδ =−1andtwovaluesofs(cid:48)(cid:48) =0.5
T zz
and 1.5. The two disks have identical radial profiles of
all fluid variables in the midplane, but different vertical
profiles of ρ and T. Clearly, this leads to an important
difference: while s(cid:48)(cid:48) = 0.5 case shows outflow of gas
zz
in the midplane and inflow at high altitudes, the s(cid:48)(cid:48) =
zz
1.5 run shows inflow at all altitudes, in agreement with
analyticalexpectations(asshowninFigure2,inthiscase
T increases with height).
In Figure 4a we show meridional profiles of the ra-
dial velocity, normalized to ν/R, at different radii for
s(cid:48)(cid:48) =0.75. Given our thermodynamic prescription (25),
zz
we expect and find to hold in simulations the follow-
ing parabolic profile of the radial velocity: u (R,z) ∝
R
(cid:2)1+C×(z/H )2(cid:3), C < 0 (c.f. Jacquet 2013). In the
0
case of δ = −1, when H = R this behavior results in
T 0
themeridionalprofileofu normalizedbyν/R,whichis
R
independentofR. Thisself-similarityisindeedobserved
in our simulations, as demonstrated by the overlapping
curves in Figure 4a.
InFigure4bwedemonstratetheconvergenceoftheu
R
profile to a steady state. One can see that convergence
happens within several tens Ω−1, which is much faster
K
than the local viscous time scale in the disk, in agree-
ment with earlier studies (Siemiginowska 1988; Kley &
Lin 1992; Rozyczka et al. 1994).
InFigure4cweshowhowthesteady-stateverticalpro-
filesofu (atfixedRandδ )varyaswechangethevalue
R T
ofs(cid:48)(cid:48) . Onecanseethatincreasingtherateatwhichgas
zz
entropy grows with height uniformly shifts the u(R,z)
profile down. Stars at z =0 show theoretical prediction
(16) for the value of s(cid:48)(cid:48) corresponding to each curve.
zz
Oursimulationresultsobviouslymatchtheoryverywell.
To test the inflow criterion (22) we run a set of such
simulations for different values of δ = −1.5,−1,−0.5.
T
Ineachrun, wemeasuretheradialcomponentofthegas
velocity in the midplane u (R,0) and record its sign.
R
The outcome of these runs is shown in Figure 5a, where
we mark solutions with u (R,0)<0 as blue stars, with
R
u (R,0) > 0 as red, and simulations with u (R,0) ≈ 0
R R
as yellow. Analytical inflow criterion (22) is shown as a
solid line. We find that for each inspected value of δ ,
T
there is a critical value of the vertical entropy gradient,
s(cid:48)(cid:48) , above which the disk exhibits the inflow at all z.
Fig. 5.—SummaryofsimulationresultsfordifferentvaluesofδT Tzhzis Figure demonstrates that this critical value is in
ands(cid:48)(cid:48) ,anddifferentradialbehaviorsofviscousα: (a): α=const,
zz very good agreement with the criterion (22).
(b)α∝R. Red,yellowandbluestarsrepresentsimulationsshow-
ingoutflow,uR(R,0)≈0,andinflowatthemidplane,correspond- We also repeated the same set of simulations for a dif-
ingly. Solidcurvesrepresentanalyticalpredictionsgivenbyequa- ferent viscosity behavior α ∝ R. Their results, together
tions (a) (22) and (b) (27). One can see excellent agreement be- withtheanalyticalcriterion(27),aresummarizedinFig-
tween theory and simulations (e.g. yellow stars fall on top of the
ure 5b. As in the case of α =const, we find very good
curves).
rion for the inflow at all z, which replaces equation (22): agreement between theory and simulations.
2(7γ−4)+(6γ−1)δ 6. DISCUSSION
s(cid:48)(cid:48) > T. (27)
zz 8+δ Ourstudyilluminatesthedirecteffectofthediskther-
T
modynamics on the pattern of meridional circulation.
This constraint is again a function of the radial temper- Ourfindingthatinviscousdiskstheentropyrapidlyris-
ature profile only. ingwithheightmayreversetheradialoutflowatthemid-
Forbothprescriptionsforα-viscosity,weperformsim- plane, resulting in inflow at all heights, is a new result.
ulations for several sets of different values of s(cid:48)z(cid:48)z and δT. Itwasnotfoundinpreviousstudiesfocusedprimarilyon
the locally isothermal disks for the following reason.
5. SIMULATIONRESULTS According to equation (4), meridional circulation is
We start by presenting the results of simulations for driven by both radial and vertical components of the
α(R)=const case. In Figure 3 we show the steady state stress tensor. Starting with the former, large negative
distribution of the radial velocity u in the R−z plane radial gradient of density at the midplane makes the
R
9
net contribution of T to u (R,0) positive, promoting thus, sets the accretion rate M˙ . As found by many nu-
rφ R
amidplaneoutflow(Takeuchi&Lin2002). Ontheother merical studies for locally isothermal disks, T is nearly
rφ
hand, athighaltitudesthedensitygradientwithrespect constant at low altitudes and falls off sharply in the
to R changes sign, so that T drives inflow at high z. corona, where plasma β parameter is of the order of
rφ
This leads to the characteristic parabolic shape of the unity, and MRI is quenched (Fromang et al. 2011; Flock
vertical u profiles, changing sign at some intermediate et al. 2011; Bai & Stone 2013). Considerably less is
R
altitude. known about the behavior of T , although some nu-
zφ
Vertical stress is a bit more subtle. As the second line merical results were presented in Fromang et al. (2011).
of equation (7) shows, it contains a contribution (first GiventhattheMRIstressbehaviordiffersfromthatof
term) proportional to the radial density gradient that the shear stress, it is not surprising that in their isother-
drives an outflow. But the second term proportional to mal simulations with no vertical field Fromang et al.
∂2ρ/∂z2 < 0 promotes inflow in non-pathological cases (2011)foundaverydifferentcirculationpattern,namely
when P decreases with R. It arises because the gas at an outflow at all altitudes, so that u does not change
0 R
higher z rotates slower than the midplane (see equation sign. Jacquet (2013) explained this result by adopting
[3]) and viscously removes angular momentum from the a phenomenological prescription for the MRI stress be-
underlying layers, forcing them to flow inwards. The ef- haviorsimilartotheonefoundbyFromangetal.(2011).
ficiency of this coupling near the midplane is regulated Themajordifferenceofthisprescription,whencompared
primarilybytheeffectoftheverticalvariationofρonΩ: to the anzatz (5), is that neither T nor T was as-
rφ zφ
densityslowlydecayingwithzresultsinlower|∂2Ω/∂z2|, sumed8 by Jacquet (2013) to scale with density near the
see equation (A3). This results in less efficient angular midplane. As a result, the logic used above to explain
momentumremovalfromthemidplanelayers,whichcan- high-altutudeinflowinisothermaldiskswithshearstress
not oppose the outflow tendency at the midplane. Such fails, and it becomes possible to have pure outflow solu-
situation holds in locally isothermal disks, whose gaus- tions even with locally isothermal EOS.
siandensityprofiledoesnotdecayrapidlyenoughtopre- It should, however, be remembered that the picture
vent an outflow near the midplane. of the meridional circulation in MRI-active disks is still
However,whentheentropyincreaseswithheightmore far from complete. For example, in their simulations
rapidly,thedensityfallsoffwithzmoresteeply,resulting withverysimilarsetup(startingwithweaktoroidalfield)
in larger vertical shear. This leads to a stronger angular Flocketal.(2011)foundabehavioroftheu (R,z)differ-
R
momentum loss from the midplane to the upper layers, ent from Fromang et al. (2011): inflow at the midplane,
making it possible to suppress the positive contribution and outflow at high z, i.e. u (R,z) profile does change
R
of the Trφ to uR(R,0). As a result, gas can flow in at all sign.
heights. In a somewhat different setup — starting with a net
Good agreement found between our analytical theory verticalfield—Suzuki&Inutsuka(2014)observedweak
and the results of viscous hydrodynamical simulations, gas outflow at the midplane, with clear inflow higher up,
evidentinFigures4c&5,suggeststhatourresultscanbe i.e. the circulation pattern typical for disks with shear
used for code testing purposes. They can also be called viscosity. This finding has been recently confirmed by
upon for designing hydrodynamical viscous (not MHD) the non-ideal MHD simulations of B´ethune et al. (2016)
simulations, which do not exhibit an outflow at the mid- and ideal MHD simulations of Zhu & Stone (2017): in
plane. This setup can be useful e.g. in the long-term both studies the midplane outflow changing to an inflow
studies of the boundary layers of accretion disks. Such athighaltitudeswasobserved. Eventhoughthispattern
simulations are numerically expensive to be performed of circulation agrees with the prediction of the conven-
with full MHD in 3D, calling for a viscous hydro ap- tional viscous disk theory (Urpin 1984), the true rea-
proachinstead. Butonemightalsowanttosuppressthe son for such behavior in these simulations is likely quite
outward fluid motions anywhere in the disk to exclude complicated. Inparticular,Zhu&Stone(2017)findthat
the effect of the boundary layer on the boundary condi- magneticstressesplaymoreimportantrolethatthether-
tions at the outer edge of the simulation domain. Our malpressuregradientsindeterminingtheverticalprofile
results demonstrate that this can be naturally achieved of Ω near the midplane. Clearly, this feature cannot be
just by making the vertical profile of entropy in the disk capturedintheframeworkofourpurelyhydrodynamical
to rapidly increase with height. model.
All our results hold for a particular stress model — ThedifferencesbetweentheaforementionedMRIstud-
shearstress—representedbyequations(5). Inpractice, ies strongly suggest that further global stratified MRI
angularmomentumtransportinhot,wellionizeddisksis simulationswithdifferentinitialfieldgeometriesandvar-
expectedtobeeffectedbytheMRI(Balbus2003). Inthe ied thermodynamics are needed to clarify both the be-
weaklyionizedregionsoftheprotoplanetarydisks,where haviorofu (R,z)andtheverticalvariationofthestress
R
theMRImaybeinactive,theVSI(Urpin&Brandenburg tensor (including T !), as well as to establish the con-
zφ
1998; Urpin 2003) has been suggested to drive angular nection between them.
momentum transport. In both cases the stress behavior
may be different from that given by equations (5); we 6.2. Circulation in the VSI-dominated disks
discuss this possibility next.
Stoll & Kley (2014, 2016) found that the weak trans-
port associated with the VSI gives rise to a near-
6.1. Circulation in the MRI-dominated disks
8 We note, however, that the vertical profile of T found by
zφ
Significant effort has been invested in quantifying the Fromang et al. (2011) is not incompatible with the prescription
T behavior for the MRI, which determines (cid:104)u (cid:105) and, (5).
rφ R ρ
10
midplane inflow of gas, switching to an outflow at high dreich 1997), but this is insufficient to revert the radial
altitudes. Even though VSI is a purely hydrodynamical outflowinthenear-midplanepartofthedisk. Thus,pas-
instability,thisbehaviorisclearlydifferentfromthatex- siveprotoplanetarydisksshouldfeatureanear-midplane
pectedindisksgovernedbytheshearstresswithisotropic outflow, if their transport is governed by equations (5).
viscosity (i.e. when theradial and vertical stresscompo- Disksheatedpredominantlybyviscousdissipationand
nents Trφ and Tzφ are characterized by the same value accretingathighM˙ arelikelytobeopticallythick; their
of ν), regardless of their entropy profile, see Figure 4c. T(R,z) profile inevitably exhibits temperature dropping
As shown in the recent work of Stoll et al. (2017), this withheight. Onceagain,ourstudyshowsthatsuchdisks
behavior is caused by a strong anisotropy of the stress should exhibit radial outflow near the midplane. Disks
tensor in the VSI-dominated disks: Tzφ was found to with isentropic stratification (s(cid:48)z(cid:48)z = 0) may feature the
exceed TRφ by a factor of several hundred. Moreover, fastest radial outflow at z = 0, since entropy dropping
the results of Stoll et al. (2017) suggest that near the with z will likely drive efficient convection enforcing ver-
midplane the vertical variation of Tzφ can be reasonably tically homogeneous entropy.
well approximated by the shear stress prescription (5). However, viscously heated disks accreting at low M˙
The radial stress T certainly does not follow the scal-
Rφ should be optically thin to their own emission (Narayan
ing T ∝ρ postulated in equations (5), as Stoll & Kley
Rφ & Popham 1993). If the specific viscous dissipation rate
(2014)andStolletal.(2017)findT toincreasewithz
Rφ increases with height (which is expected, for example, if
near the midplane, switching to a decay only at high al-
heating is produced by the MRI (Miller & Stone 2000;
titudes. However,giventhenegligibleroleplayedbyT
Rφ Hirose et al. 2006)), then gas temperature will increase
in the VSI-dominated disks, we can nevertheless use our
with z at all altitudes. This can potentially provide the
results to understand the meridional circulation pattern
conditions favorable for the radial inflow at all heights,
observed in simulations of such disks.
as we showed in this work. Thus, optically thin, vis-
Indeed, let us look at the equation (16), in which we
couslyheatedaccretiondiskspresentthebestsettingfor
will drop the terms in the first line, as they result from
suppressing the radial outflow at all altitudes.
thenegligibleT . Wealsosets(cid:48)(cid:48) =γ−1asappropriate
Rφ zz
fortheverticallyisothermaldiskstructureoftenadopted 7. SUMMARY
intheVSIsimulations(Stolletal.2017). Thentheterms
We explored meridional circulation in accretion disks
inthesecondlineofequation(16),comingfromtheT ,
zφ with shear viscosity and varied thermodynamics. While
naturally result in u(R,0) = ν δ /R, where ν is the
z T z previous studies of this problem focused on the verti-
valueofthekinematicviscositycoefficientcharacterizing
cally isothermal disks, finding radial gas outflow at the
the vertical behavior of the T as found in Kley et al.
zφ midplane and inflow at high latitudes, we demonstrate
Given that the disk temperature decreases with R (i.e.
that different assumptions about thermal stratification
δ < 0), one finds that u(R,0) < 0 (meaning midplane
T can change this pattern. We show that the direction of
inflow) in the VSI-mediated disks dominated by the ver-
the flow at the midplane is intimately connected to the
tical stress. Thus, our analytical calculations provide a
behavior of the vertical stress T . Vertical density pro-
natural explanation for the meridional circulation pat- zφ
files steeply falling off with height induce significant ver-
tern found in simulations of such disks.
ticalshear, whichcanmakethegasatthediskmidplane
Despite all the complications related to stress
to flow inward. We derive analytical criterion relating
anisotropy, based on our results, we would still expect
the direction of the midplane flow to the thermal strati-
thermal stratification to have an important effect on the
fication inthe disk and show that theradial inflow atall
meridionalcirculationin3Ddiskseveniftheangularmo-
heights (without change of sign of u ) naturally sets in
mentumtransportismediatedbymechanismotherthan R
diskswithsteeplygrowingverticalprofilesofentropyand
the shear viscosity.
temperature. Such conditions can be naturally realized
in optically thin disks heated primarily by viscous dissi-
pation. Although our findings rely on the assumption of
6.3. Applicability to real astrophysical systems
shearviscosity,wealsocommentonothermechanismsof
Our study shows that whenever the angular momen- angular momentum transport. Our results can also be
tum transport in the disk is effected by shear viscos- used for code testing and designing simulations with the
ity, the radial inflow at all altitudes above the disk mid- prescribed pattern of meridional circulation.
planenecessarilyrequiresbothtemperatureandentropy
to rapidly increase with |z|. In light of this result, it
is natural to ask, in which astrophysical systems such We are grateful to Willy Kley for useful discussions.
conditions could hold. A.A.P.issupportedbyPorterOgdenJacobusFellowship,
Protoplanetary disks, heated predominantly by stellar awarded by the graduate school of Princeton Univer-
irradiation,haveroughlyisothermalverticaltemperature sity. Financial support for this study has been provided
profile within several scaleheights above the midplane by the NSF via grant AST-1515763, NASA via grants
(Chiang & Goldreich 1997). This expectation motivated 14-ATP14-0059, and The Ambrose Monell Foundation.
the adoption of the isothermal thermodynamic setup in Simulations presented in this article used computational
many previous studies of meridional circulation (Urpin resourcessupportedbythePICSciE-OITTIGRESSHigh
1984; Ciesla 2007; Jacquet 2013). Presence of the super- Performance Computing Center and by NSF through an
heateddustlayerhighabovethemidplaneshouldeventu- XSEDE computational time allocation TG-AST160008
ally result in temperature rise at some z (Chiang & Gol- on TACC Stampede and Ranch.