Radial Transport and Meridional Circulation in Accretion Disks

# Radial Transport and Meridional Circulation in Accretion Disks

Alexander A. Philippov11affiliation: Department of Astrophysical Sciences, Princeton University, Ivy Lane, Princeton, NJ 08540 44affiliation: sashaph@princeton.edu and Roman R. Rafikov22affiliation: Institute for Advanced Study, Einstein Drive, Princeton, NJ 08540 33affiliation: Centre for Mathematical Sciences, Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK
###### Abstract

Radial transport of particles, elements and fluid driven by internal stresses in three-dimensional (3D) astrophysical accretion disks is an important phenomenon, potentially relevant for the outward dust transport in protoplanetary disks, origin of the refractory particles in comets, isotopic equilibration in the Earth-Moon system, etc. To gain better insight into these processes, we explore the dependence of meridional circulation in 3D disks with shear viscosity on their thermal stratification, and demonstrate strong effect of the latter on the radial flow. Previous locally isothermal studies have normally found a pattern of the radial outflow near the midplane, switching to inflow higher up. Here we show, both analytically and numerically, that a flow, which is inward at all altitudes, is possible in disks with entropy and temperature steeply increasing with height. Such thermodynamic conditions may be typical in the optically thin, viscously heated accretion disks. Disks in which these conditions do not hold should feature radial outflow near the midplane, as long as their internal stress is provided by the shear viscosity. Our results can also be used for designing hydrodynamical disk simulations with a prescribed pattern of the meridional circulation.

accretion, accretion disks — protoplanetary disks — hydrodynamics

## 1. Introduction

Evolution of astrophysical accretion disks is believed to be driven primarily by their internal stresses (Papaloizou & Lin, 1995). Gross features of this process can be understood by treating the disk as a geometrically thin structure and characterizing its properties using the vertically integrated (or properly averaged across the disk thickness) variables such as the surface density (Shakura & Sunyaev, 1973; Lynden-Bell & Pringle, 1974). In this approach the steady accretion disks with the radially constant mass accretion rate necessarily exhibit purely inward radial motion of gas, driven by the angular momentum redistribution in the disk due to its internal stresses. In other words, the density-weighted vertical average of the radial velocity is always negative, .

Situation becomes more complicated once one considers the full three-dimensional disk structure, in particular, the profile of the radial velocity as a function of the vertical coordinate . It was first shown by Urpin (1984) that when the angular momentum transport in the disk is effected by the effective shear viscosity (Shakura & Sunyaev, 1973), the radial velocity of the gas is actually positive at , implying a radial ouflow near the disk midplane. The magnitude of positive steadily decreases with the vertical distance from the midplane, and ultimately changes sign at some altitude of order the disk scale height . As a result, a meridional circulation pattern sets in the poloidal plane of the disk. Radial gas inflow at high carries more mass inward than is transported out by the midplane outflow, resulting in a net (vertically integrated) inflow of mass towards the accreting object (Urpin, 1984; Jacquet, 2013). Thus, even though is still negative, the radial outflow near the midplane actually transports mass out in 3D. This remarkable analytical prediction was later confirmed by numerical hydrodynamical calculations of Siemiginowska (1988), Kley & Lin (1992), and Rozyczka et al. (1994), which explicitly employed -prescription (Shakura & Sunyaev, 1973) to describe the internal stress.

The meridional circulation, if indeed present in accretion disks, should have very important implications. It may provide a natural way of coupling the inner disk to its outer parts: the material and information from the vicinity of an accreting object could potentially be directly transported in the advective fashion farther out in the disk by the near-midplane outflow.

To provide some illustrative examples, Takeuchi & Lin (2002) have suggested that meridional circulation can drive the outward transport of dust grains in protoplanetary disks. More recently, the samples collected by the Stardust mission from comet 81P/Wild 2 revealed the presence of large (m) grains composed of high-temperature minerals that appear to have formed in the inner regions of the Solar nebula (Brownlee et al., 2006; Zolensky et al., 2006). Their presence in comets that form at tens of AU indicates that global outward radial transport of dust particles could have operated in the Solar nebula (Ciesla, 2007, 2009; Hughes & Armitage, 2010; Jacquet, 2013).

Moon is thought to have formed in a collision of a Mars-size planetary embryo with the proto-Earth (Hartmann & Davis, 1975), primarily out of the material derived from the impactor (Canup, 2004). This should have resulted in the different compositions of the Moon and the Earth. However, the isotopic ratios of oxygen derived from lunar Apollo samples were found to be essentially identical111Terrestrial O isotopic ratios are quite distinct from e.g. their martian values. to their terrestrial analogs (Wiechert et al., 2001). To explain this puzzle Pahlevan & Stevenson (2007) suggested that equilibration of oxygen isotopes in the Earth-Moon system resulted from the rapid radial mixing between the terrestrial magma ocean and the proto-Lunar disk of vapor-melt debris produced by the impact. Pahlevan & Stevenson (2007) appealed to diffusive turbulent mixing as the mechanism of compositional equilibration inside the disk. However, meridional circulation can also be an important contributor to this process by directly advecting into the disk and mixing the material from the magma ocean via the midplane outflow.

Weakly magnetized accreting objects (cataclysmic variables, classical nova progenitors, neutron stars in low-mass X-ray binaries, etc.) are expected to develop a boundary layer between the inner edge of the accretion disk and the stellar surface (Popham et al., 1993; Popham & Narayan, 1995; Belyaev et al., 2012, 2013a, 2013b). If the material from the near-surface layers of the accretor can be dredged up (evidence for this exists in e.g. classical novae Truran & Livio 1986; Gehrz et al. 1998) into the boundary layer by some internal processes, then the near-midplane outflow could subsequently transport it to the outer parts of the disk. Such elemental pollution might spread the metals with unusual abundances across the disk, affecting its observational appearance (e.g. via spectral signatures).

Given the variety of situations, in which meridional circulation may be important, it is natural to ask how robust this phenomenon is. The majority of its analytical and numerical investigations focused on disks, in which angular momentum transport is accomplished by the conventional -viscosity (Shakura & Sunyaev, 1973). All of them invariably find implying the near-midplane outflow. However, in many types of disks the transport is more likely to be effected by the magnetorotational instability (MRI, Balbus 2003), for which the stress tensor is anisotropic. The anisotropy of has been shown (Jacquet, 2013) to play an important role for the meridional circulation, although the numerical evidence on the issue is rather mixed at this point (Fromang et al., 2011; Flock et al., 2011; Suzuki & Inutsuka, 2014), see the discussion in §6. But even without magnetic fields, global simulations by Stoll & Kley (2014, 2016) of the purely hydrodynamic vertical shear instability (VSI; Urpin & Brandenburg 1998; Urpin 2003) exhibit radial inflow at the midplane, changing to an outflow at high , contrary to the expectations for disks mediated by the shear stress (Urpin, 1984; Takeuchi & Lin, 2002).

In this study we address a different aspect of the problem and explore the effect of thermal properties of the disk on the meridional circulation, while still confining ourselves to the shear viscosity as the source of the angular momentum transport (which allows us to perform detailed analytical investigation). Existing analytical and numerical work has typically assumed the disk to have locally isothermal vertical structure (Urpin, 1984; Fromang et al., 2011; Jacquet, 2013; Suzuki & Inutsuka, 2014). This may be a reasonable assumption for the externally irradiated protoplanetary disks, but should not generally hold in accretion disks dominated by viscous dissipation. Some earlier numerical studies (Siemiginowska, 1988; Kley & Lin, 1992; Rozyczka et al., 1994) have attempted to use more refined treatments of the disk thermodynamics. However, they were limited in resolution and did not explore the full range of possible thermodynamic regimes. Our present goal is to provide a thorough analysis of the effects of thermal stratification on the meridional circulation in the disk.

This work is organized as follows. In §2 we present general theoretical description of meridional circulation in accretion disks. We explore conditions, under which the inflow occurs at the disk midplane in §3. In §4-5 we verify our theoretical results with 3D viscous numerical simulations. We discuss our findings in §6, focusing on applications to real systems in §6.3. Our results are summarized in §7.

## 2. Theoretical considerations

We consider a geometrically thin, axisymmetric disk orbiting in a central potential of a point mass ; we work in cylindrical coordinates . We are interested in the spatial structure of the radial velocity , which is non zero because of the internal stress operating in the disk. We are interested only in advective, laminar motions of the fluid and do not consider turbulent diffusion, which is often invoked in studies of particle transport in disks (Takeuchi & Lin, 2002; Hughes & Armitage, 2010). Thermal structure of the disk is specified in §2.1.

Equations of the radial and vertical balance, describing the equilibrium disk structure, are

 Ω2R−1ρ∂P∂R = ∂Φ∂R, (1) 1ρ∂P∂z = −∂Φ∂z, (2)

where and are the gas density and pressure, is the angular frequency, and is the central potential (we neglect the disk self-gravity). In these equations we also neglected viscous terms resulting from the fluid motions in the meridional plane; we comment on the significance of this simplification later.

Because of the variation of the central gravity with height and radial pressure support in the disk the angular frequency deviates from its Keplerian value :

 Ω(R,z)≈ΩK(1−34z2R2+12Ω2KR2Pρ∂lnP∂lnR). (3)

This approximate relation follows from our assumption of a geometrically thin disk, for which , , where is the local disk scaleheight, and is the isothermal sound speed. Note that the last term in parentheses is a function of both222In particular, one can show with the aid of equation (2) that is independent of for a barotropic equation of state. and . The deviations from purely Keplerian rotation in equation (3) are at the (second term in parentheses) and (third term) level.

The -component of the equation of motion describing the angular momentum balance can be reduced (under reasonable assumptions regarding the amplitudes of the radial and vertical fluid velocities, Takeuchi & Lin 2002; Fromang et al. 2011) to the following expression for the radial velocity of the disk fluid (Jacquet, 2013):

 uR(R,z) = −1ρ(∂l∂R)−1 (4) × [1R∂∂R(R2TRϕ)+∂∂z(RTzϕ)].

Here is the specific angular momentum of the disk fluid, while and represent the and components of the internal stress tensor. This expression is accurate at the level; it neglects certain terms proportional to the vertical velocity and possible time variability of the azimuthal velocity (Fromang et al., 2011). It demonstrates that the radial fluid motion of the 3D disk is determined not only by the horizontal but also by the vertical component of the stress tensor. This turns out to be of crucial importance for our work. On the other hand, multiplying equation (4) by and integrating it over , one finds that is determined by the behavior only (Jacquet, 2013).

In real astrophysical disks internal stress can be provided by a variety of mechanisms, with the MRI (Balbus, 2003) being one of the important possibilities. However, our present focus is still on standard shear viscosity prescription, for which one has

 TRϕ=−ρνR∂Ω∂R,    Tzϕ=−ρνR∂Ω∂z, (5)

where is a kinematic viscosity.

Previous studies (Fromang et al., 2011; Jacquet, 2013) clearly demonstrated that the possible meridional outflow is most pronounced at the midplane of the disk, at . For that reason, our primary goal will be to determine the factors that affect the midplane value of the radial velocity, . To that effect, in our analysis we will often take the limit in evaluating different expressions and account for the symmetry of the disk with respect to its midplane. The latter implies, in particular that the first derivatives of various fluid variables with respect to , e.g. , , etc., go to zero as .

In Appendix A we show that after some straightforward manipulations the midplane value of the radial velocity can be written as

 uR(R,0) = νR[152−3∂ln(ρ0ν0R3)∂R (6) − δρ−P0Ω2Kρ20∂lnP0∂lnR×∂2ρ∂z2∣∣∣z→0],

where

 δρ≡∂lnρ0/∂lnR, (7)

and subscript implies the value of a particular variable at the midplane, e.g. , , etc. Definition (7) implies that is closely related to the density power law exponent used in previous studies (Takeuchi & Lin, 2002; Fromang et al., 2011). The two coincide if , but the power law density profile is not necessary for our results to be valid. Note that the expression (6) does not make any assumptions about the vertical or radial dependence of the viscosity, although it does rely on the stress model (5).

The two terms in the first line of the equation (6) originate from , while the last two terms are due to the vertical stress . Since the midplane pressure is a decreasing function of , it is clear that a vertical density profile steeply declining with (i.e. high value of ) should be favorable for suppressing the midplane outflow in the disk and driving gas inflow () at all heights.

Note that equations (1)-(2) neglect contributions from viscous stresses arising from shear in the meridional plane. Such terms were considered in studies of Kley & Lin (1992), Kluzniak & Kita (2000), and Regev & Gitelman (2002), resulting in the additional relative contribution to the expression (3) for the angular frequency. Their inclusion would, in turn, lead to the emergence of an additional contribution in the expression (6) for . Given the expectation of in real disks, the omission of the viscous terms should not present a problem333Kluzniak & Kita (2000) and Regev & Gitelman (2002) find that inclusion of the viscous terms changes midplane outflow to an inflow only for , although Kley & Lin (1992) suggest a lower value of critical for this transition.. This expectation is subsequently confirmed in §5 by the good agreement between our analytical theory (which explicitly neglects viscous contributions to ) and numerical results.

### 2.1. Disk thermodynamics

So far our treatment was fully general. To make further progress we need to look in more detail into the disk thermodynamics. We will assume a rather general equation of state (EOS) for the disk fluid in the polytropic form

 P=esργ,   s(R,Z)=c−1VS(R,z), (8)

where is the polytropic index and is the scaled (dimensionless) gas entropy ; is the specific heat capacity. Note that in this work is specified explicitly, i.e. we do not attempt to calculate it, e.g. from the energy equation by solving for the radiation transfer in the disk.

An often used assumption of a locally isothermal disk structure (Takeuchi & Lin, 2002; Fromang et al., 2011; Jacquet, 2013) implies that , so that

 siso(R,z)=s0(R)+γ−12[zH0(R)]2, (9)

where . A different limit of the locally isentropic disks considered later in §3.2.1 has .

Using the ansatz (8) we show in Appendix A that the last term inside the brackets in the expression (6) can be written as

 δρ+δTγ(1+PΩ2Kρ∂2s∂z2∣∣∣z→0), (10)

where

 δT≡∂lnT0/∂lnR, (11)

and is the midplane value of the gas temperature. Definition of makes it closely related to the temperature power law index used in other studies of the meridional circulation in accretion disks (Takeuchi & Lin, 2002; Fromang et al., 2011), which assume .

Next we connect the behavior of the second term in the equation (6) to the global structure of the disk. We introduce the viscous angular momentum flux across a given radius , which is equal to the viscous torque exerted by the inner disk on the outer disk. The concept of is known to be very useful for describing viscous accretion disks in a variety of situations, especially a steady state (Lynden-Bell & Pringle, 1974; Rafikov, 2013, 2016a). According to the definition, can be calculated by multiplying in equation (5) by and integrating over . Using equation (5) this implies, to accuracy, that

 FJ=−∞∫−∞2πR3ρν∂Ω∂Rdz=3πl∞∫−∞ρνdz, (12)

where and we have assumed a Keplerian rotation profile. Also, multiplying equation (4) by and integrating over one finds that the mass accretion rate (defined to be positive for inflow) is related to as (Rafikov, 2013)

 ˙M(R)=−2πR∞∫−∞ρuRdz=∂FJ∂l, (13)

where the derivative is with respect to the specific angular momentum, so that .

We now make an additional assumption that the disk structure obeys a certain similarity property, namely that

 S(R,z)=S(R,zH0(R)),    H0(R)≡H(R,0), (14)

where is the local value of the scaleheight at the disk midplane. This is not a very constraining assumption and its adoption should not limit the applicability of our subsequent results.

With such entropy behavior it follows from equation (8) that and also depend on only in combination ; this makes it natural to expect that as well. Then equation (12) implies that

 FJ∝ρ0ν0H0l∝ρ0ν0T1/20R2, (15)

allowing us to tackle the second term in equation (6).

Plugging equations (10) and (15) into the expression (6) we finally find that

 uR(R,0) = νR[92−3(δF−δT2) (16) − δρ+δρ+δTγ(1+s′′zz)],

where we introduce the shorthand notation

 s′′zz≡H20∂2s∂z2∣∣∣z→0 (17)

for the (dimensionless) second derivative of the entropy at the disk midplane, and

 δF≡∂lnFJ/∂lnR. (18)

Once again, in equation (16) the first two terms arise from , while the terms in the second line originate from .

Equation (16) represents the main analytical result of this work. It provides a connection between the amplitude (and direction) of the radial velocity at the disk midplane and the vertical thermal stratification characterized by . This link has not been established in previous studies of meridional circulation.

## 3. Conditions for inflow at all z

Equation (16) allows us to determine the conditions under which the disk will exhibit inflow at all , including the midplane. Setting one finds a necessary criterion for this to be the case:

 s′′zz >−γδρ+δT (19) ×(92−3δF+3γ+22γδT+1−γγδρ),

where we assumed that decreases with so that . When this inequality is fulfilled meridional circulation is unable to convey information and material from the inner disk to the outer disk — only inward propagation is allowed.

Equation (19) is the most general form of the inflow criterion that does not make any assumptions about the values of , , and — they can take arbitrary values allowing one to explore meridional circulation even in evolving disks. Moreover, even though it does assume a particular stress model given by the equation (5), it makes only a weak assumption about the actual behavior of the viscosity (that depends on only in combination ).

### 3.1. Inflow criterion for a standard constant ˙M= disk

In practice, one is often interested in accretion disks that have reached a steady state. One of the most popular assumptions used in many studies is that of the radially constant mass accretion rate through the disk (const), with no torque applied at its center. In this case, integrating equation (13) one finds (Rafikov, 2013), so that . This transforms inequality (19) into the following inflow criterion:

 s′′zz>−γδρ+δT(3+3γ+22γδT+1−γγδρ). (20)

Note that the constant assumption does not constrain or because of the freedom in choosing the radial variation of the viscosity . This constraint is illustrated in Figure 1 for two values of ( and ) and three values of . One can see that steeper decay of with (more negative ) makes it easier to achieve inflow at the disk midplane, i.e. requires less extreme values of .

We can now go one step further and adopt a particular viscosity ansatz, namely the -model of Shakura & Sunyaev (1973), in which . If we additionally assume that the effective viscosity parameter is independent of , then locally . Plugging this and into equation (15) one obtains that

 δρ=−3−32δT. (21)

This results in yet another version of the inflow criterion for disks, which depends only on the radial temperature profile (i.e. ):

 s′′zz>6(2γ−1)+(6γ−1)δT6+δT. (22)

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 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 isentropic disk, in which entropy (and ) does not depend on . Such disks naturally have temperature decreasing with height. The isentropic vertical stratification may arise e.g. if the disk is convectively unstable.

In this case the left hand side of equation (20) becomes zero and the condition of inflow turns into a constraint on the radial behavior of and in the form

 δρ>6γ+(3γ+2)δT2(γ−1). (23)

For protoplanetary disks with this becomes . In hotter accretion disks with the constraint is .

It is clear that even if decays with radius as rapidly as (i.e. ) the inflow at the disk midplane would require midplane density increasing outwards. This is a pretty unusual arrangement, which makes midplane inflow essentially impossible in locally isentropic disks, including the disks which are convectively unstable.

#### 3.2.2 Locally isothermal disk

A popular assumption of the locally isothermal disk structure () turns the left hand side of the inflow criterion into , see equation (9). Then it is easy to see from the inequality (20) that the inflow at all heights is possible only if . This result is completely independent of either or , thus it should hold for an arbitrary radial profile of the midplane density .

Finding conditions in which a disk might have decaying with radius faster than is not easy. Based on this we can conclude that the locally isothermal constant disks are predetermined to naturally exhibit an outflow at the midplane, in agreement with a number of previous studies (Urpin, 1984; Jacquet, 2013).

Summarizing the results of this and previous subsection (§3.2.1), we conclude that the inflow at requires the disk to have temperature increasing with height, starting from the midplane. Such disks naturally have entropy rising more steeply with than in the isothermal case and may satisfy the conditions (20) and (22) without making unrealistic assumptions about the behavior of and .

### 3.3. Inflow criterion for Fj= const disk

Standard constant- disk with zero central torque represents just one particular example of a steady disk. A more general time-independent structure of the disk is described by (Rafikov, 2013, 2016a), where represents the torque applied at the disk center (which is naturally absent in a standard constant- case).

A particularly interesting case to explore is that of a disk with no mass accretion at the center, in which the gas inflow is fully suppressed by a strong central torque (Lynden-Bell & Pringle, 1974; Pringle, 1991; Rafikov, 2016b). Such ”dead” disks were first studied by 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 const, implying that .

Equation (16) demonstrates that reaching inflow at the midplane of such a disk is more difficult than in the standard constant- disk. Indeed, the inflow criterion (20) remains essentially the same, however the first term inside the parentheses changes from 3 to 9/2. As a result, a higher value of is needed to guarantee inflow at all heights (i.e. ) in the constant- disk, as compared to the standard constant disk (§3.1). Inequality (22) gets modified in a similar fashion.

Moreover, repeating the calculations of §3.2.2 one finds that the inflow at in a locally isothermal, constant- disk would require . Such steeply declining radial profiles of temperature are unlikely in real disks.

## 4. Numerical setup

In order to confirm the analytical result (16), we performed viscous hydrodynamical simulations using new Godunov code Athena++ code (Stone et. al. 2016, in preparation). Compared with its predecessor Athena (Gardiner & Stone, 2008), Athena++ is highly optimized and uses flexible grid structures, significantly facilitating global numerical simulations. In this work we perform our simulations in spherical coordinates, using uniform grid in and constraining ourselves to small range in . We use the domain in these coordinates with numerical resolution of cells. We verified that simulations in the full range lead to same results. We also verified that our results are converged with regard to numerical resolution.

For 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 in all our simulations. The simulation time is always measured in units of , where is the Keplerian angular velocity at the inner boundary.

In order for the disk structure to remain unchanged and maintain its initially prescribed entropy profile during the simulation, we use the optically thin cooling function:

 Λ=−ρT−T0τ, (24)

where is the gas temperature, is the initial temperature profile described in §4.1, is the cooling time. This cooling prescription ensures that viscous heating does not significantly influence our initial disk configuration, and that our simulations reach steady state.

### 4.1. Initial entropy and density profiles

For simplicity, we assume the power-law behavior of the density and temperature at the mid-plane , . Our simulations adopt a setup typical for standard constant disks (i.e. we do not simulate disks considered in §3.3), making different assumptions about the radial profile of . Close to the mid-plane, we consider the following entropy profile

 s(R,z)=s0(R)+s′′zz2[zH0(R)]2, (25)

such that the definition (17) holds true; also . In order to prevent large gradients of fluid variables from developing high above the mid-plane, we modify this entropy behavior in such a way that the disk becomes locally isothermal at .

To compute the vertical disk structure corresponding to this entropy behavior, we numerically integrate the equation of hydrostatic equilibrium (2). The azimuthal component of the gas velocity is then computed using equation (1). We show disk profiles for , constant and two values of in Figure 2. It shows that for high values of the entropy derivative the disk naturally develops temperature profile rising with height. As we show below, in this case one finds gas inflow at the disk midplane. On the contrary, low values of result in dropping with (by design, at high altitudes our profiles always converge to isothermal).

To check our analytical prediction (16), we run two sets of simulations. In a first set, we adopt a radially constant profile with , so that equations (21) and (22) hold. In a second set, we assume that , with at . In this case, the constant assumption leads to a different relation between the power-law indices of density and temperature profiles:

 δρ=−4−32δT. (26)

Plugging this in equation (20), we obtain a new criterion for the inflow at all , which replaces equation (22):

 s′′zz>2(7γ−4)+(6γ−1)δT8+δT. (27)

This constraint is again a function of the radial temperature profile only.

For both prescriptions for -viscosity, we perform simulations for several sets of different values of and .

## 5. Simulation results

We start by presenting the results of simulations for const case. In Figure 3 we show the steady state distribution of the radial velocity in the plane for profile with and two values of and 1.5. The two disks have identical radial profiles of all fluid variables in the midplane, but different vertical profiles of and . Clearly, this leads to an important difference: while case shows outflow of gas in the midplane and inflow at high altitudes, the run shows inflow at all altitudes, in agreement with analytical expectations (as shown in Figure 2, in this case increases with height).

In Figure 4a we show meridional profiles of the radial velocity, normalized to , at different radii for . Given our thermodynamic prescription (25), we expect and find to hold in simulations the following parabolic profile of the radial velocity: , (c.f. Jacquet 2013). In the case of , when this behavior results in the meridional profile of normalized by , which is independent of . This self-similarity is indeed observed in our simulations, as demonstrated by the overlapping curves in Figure 4a.

In Figure 4b we demonstrate the convergence of the profile to a steady state. One can see that convergence happens within several tens , which is much faster than the local viscous time scale in the disk, in agreement with earlier studies (Siemiginowska, 1988; Kley & Lin, 1992; Rozyczka et al., 1994).

In Figure 4c we show how the steady-state vertical profiles of (at fixed and ) vary as we change the value of . One can see that increasing the rate at which gas entropy grows with height uniformly shifts the profile down. Stars at show theoretical prediction (16) for the value of corresponding to each curve. Our simulation results obviously match theory very well.

To test the inflow criterion (22) we run a set of such simulations for different values of . In each run, we measure the radial component of the gas velocity in the midplane and record its sign. The outcome of these runs is shown in Figure 5a, where we mark solutions with as blue stars, with as red, and simulations with as yellow. Analytical inflow criterion (22) is shown as a solid line. We find that for each inspected value of , there is a critical value of the vertical entropy gradient, , above which the disk exhibits the inflow at all . This Figure demonstrates that this critical value is in very good agreement with the criterion (22).

We also repeated the same set of simulations for a different viscosity behavior . Their results, together with the analytical criterion (27), are summarized in Figure 5b. As in the case of const, we find very good agreement between theory and simulations.

## 6. Discussion

Our study illuminates the direct effect of the disk thermodynamics on the pattern of meridional circulation. Our finding that in viscous disks the entropy rapidly rising with height may reverse the radial outflow at the midplane, resulting in inflow at all heights, is a new result. It was not found in previous studies focused primarily on the locally isothermal disks for the following reason.

According to equation (4), meridional circulation is driven by both radial and vertical components of the stress tensor. Starting with the former, large negative radial gradient of density at the midplane makes the net contribution of to positive, promoting a midplane outflow (Takeuchi & Lin, 2002). On the other hand, at high altitudes the density gradient with respect to changes sign, so that drives inflow at high . This leads to the characteristic parabolic shape of the vertical profiles, changing sign at some intermediate altitude.

Vertical stress is a bit more subtle. As the second line of equation (7) shows, it contains a contribution (first term) proportional to the radial density gradient that drives an outflow. But the second term proportional to promotes inflow in non-pathological cases when decreases with . It arises because the gas at higher rotates slower than the midplane (see equation [3]) and viscously removes angular momentum from the underlying layers, forcing them to flow inwards. The efficiency of this coupling near the midplane is regulated primarily by the effect of the vertical variation of on : density slowly decaying with results in lower , see equation (A3). This results in less efficient angular momentum removal from the midplane layers, which cannot oppose the outflow tendency at the midplane. Such situation holds in locally isothermal disks, whose gaussian density profile does not decay rapidly enough to prevent an outflow near the midplane.

However, when the entropy increases with height more rapidly, the density falls off with more steeply, resulting in larger vertical shear. This leads to a stronger angular momentum loss from the midplane to the upper layers, making it possible to suppress the positive contribution of the to . As a result, gas can flow in at all heights.

Good agreement found between our analytical theory and the results of viscous hydrodynamical simulations, evident in Figures 4c & 5, suggests that our results can be used for code testing purposes. They can also be called upon for designing hydrodynamical viscous (not MHD) simulations, which do not exhibit an outflow at the midplane. This setup can be useful e.g. in the long-term studies of the boundary layers of accretion disks. Such simulations are numerically expensive to be performed with full MHD in 3D, calling for a viscous hydro approach instead. But one might also want to suppress the outward fluid motions anywhere in the disk to exclude the effect of the boundary layer on the boundary conditions at the outer edge of the simulation domain. Our results demonstrate that this can be naturally achieved just by making the vertical profile of entropy in the disk to rapidly increase with height.

All our results hold for a particular stress model — shear stress — represented by equations (5). In practice, angular momentum transport in hot, well ionized disks is expected to be effected by the MRI (Balbus, 2003). In the weakly ionized regions of the protoplanetary disks, where the MRI may be inactive, the VSI (Urpin & Brandenburg, 1998; Urpin, 2003) has been suggested to drive angular momentum transport. In both cases the stress behavior may be different from that given by equations (5); we discuss this possibility next.

### 6.1. Circulation in the MRI-dominated disks

Significant effort has been invested in quantifying the behavior for the MRI, which determines and, thus, sets the accretion rate . As found by many numerical studies for locally isothermal disks, is nearly constant at low altitudes and falls off sharply in the corona, where plasma parameter is of the order of unity, and MRI is quenched (Fromang et al., 2011; Flock et al., 2011; Bai & Stone, 2013). Considerably less is known about the behavior of , although some numerical results were presented in Fromang et al. (2011).

Given that the MRI stress behavior differs from that of the shear stress, it is not surprising that in their isothermal simulations with no vertical field Fromang et al. (2011) found a very different circulation pattern, namely an outflow at all altitudes, so that does not change sign. Jacquet (2013) explained this result by adopting a phenomenological prescription for the MRI stress behavior similar to the one found by Fromang et al. (2011). The major difference of this prescription, when compared to the anzatz (5), is that neither nor was assumed444We note, however, that the vertical profile of found by Fromang et al. (2011) is not incompatible with the prescription (5). by Jacquet (2013) to scale with density near the midplane. As a result, the logic used above to explain high-altutude inflow in isothermal disks with shear stress fails, and it becomes possible to have pure outflow solutions even with locally isothermal EOS.

It should, however, be remembered that the picture of the meridional circulation in MRI-active disks is still far from complete. For example, in their simulations with very similar setup (starting with weak toroidal field) Flock et al. (2011) found a behavior of the different from Fromang et al. (2011): inflow at the midplane, and outflow at high , i.e. profile does change sign.

In a somewhat different setup — starting with a net vertical field — Suzuki & Inutsuka (2014) observed weak gas outflow at the midplane, with clear inflow higher up, i.e. the circulation pattern typical for disks with shear viscosity. This finding has been recently confirmed by the non-ideal MHD simulations of Béthune et al. (2016) and ideal MHD simulations of Zhu & Stone (2017): in both studies the midplane outflow changing to an inflow at high altitudes was observed. Even though this pattern of circulation agrees with the prediction of the conventional viscous disk theory (Urpin, 1984), the true reason for such behavior in these simulations is likely quite complicated. In particular, Zhu & Stone (2017) find that magnetic stresses play more important role that the thermal pressure gradients in determining the vertical profile of near the midplane. Clearly, this feature cannot be captured in the framework of our purely hydrodynamical model.

The differences between the aforementioned MRI studies strongly suggest that further global stratified MRI simulations with different initial field geometries and varied thermodynamics are needed to clarify both the behavior of and the vertical variation of the stress tensor (including !), as well as to establish the connection between them.

### 6.2. Circulation in the VSI-dominated disks

Stoll & Kley (2014, 2016) found that the weak transport associated with the VSI gives rise to a near-midplane inflow of gas, switching to an outflow at high altitudes. Even though VSI is a purely hydrodynamical instability, this behavior is clearly different from that expected in disks governed by the shear stress with isotropic viscosity (i.e. when the radial and vertical stress components and are characterized by the same value of ), regardless of their entropy profile, see Figure 4c.

As shown in the recent work of Stoll et al. (2017), this behavior is caused by a strong anisotropy of the stress tensor in the VSI-dominated disks: was found to exceed by a factor of several hundred. Moreover, the results of Stoll et al. (2017) suggest that near the midplane the vertical variation of can be reasonably well approximated by the shear stress prescription (5). The radial stress certainly does not follow the scaling postulated in equations (5), as Stoll & Kley (2014) and Stoll et al. (2017) find to increase with near the midplane, switching to a decay only at high altitudes. However, given the negligible role played by in the VSI-dominated disks, we can nevertheless use our results to understand the meridional circulation pattern observed in simulations of such disks.

Indeed, let us look at the equation (16), in which we will drop the terms in the first line, as they result from the negligible . We also set as appropriate for the vertically isothermal disk structure often adopted in the VSI simulations (Stoll et al., 2017). Then the terms in the second line of equation (16), coming from the , naturally result in , where is the value of the kinematic viscosity coefficient characterizing the vertical behavior of the as found in Kley et al. Given that the disk temperature decreases with (i.e. ), one finds that (meaning midplane inflow) in the VSI-mediated disks dominated by the vertical stress. Thus, our analytical calculations provide a natural explanation for the meridional circulation pattern found in simulations of such disks.

Despite all the complications related to stress anisotropy, based on our results, we would still expect thermal stratification to have an important effect on the meridional circulation in 3D disks even if the angular momentum transport is mediated by mechanism other than the shear viscosity.

### 6.3. Applicability to real astrophysical systems

Our study shows that whenever the angular momentum transport in the disk is effected by shear viscosity, the radial inflow at all altitudes above the disk midplane necessarily requires both temperature and entropy to rapidly increase with . In light of this result, it is natural to ask, in which astrophysical systems such conditions could hold.

Protoplanetary disks, heated predominantly by stellar irradiation, have roughly isothermal vertical temperature profile within several scaleheights above the midplane (Chiang & Goldreich, 1997). This expectation motivated the adoption of the isothermal thermodynamic setup in many previous studies of meridional circulation (Urpin, 1984; Ciesla, 2007; Jacquet, 2013). Presence of the superheated dust layer high above the midplane should eventually result in temperature rise at some (Chiang & Goldreich, 1997), but this is insufficient to revert the radial outflow in the near-midplane part of the disk. Thus, passive protoplanetary disks should feature a near-midplane outflow, if their transport is governed by equations (5).

Disks heated predominantly by viscous dissipation and accreting at high are likely to be optically thick; their profile inevitably exhibits temperature dropping with height. Once again, our study shows that such disks should exhibit radial outflow near the midplane. Disks with isentropic stratification () may feature the fastest radial outflow at , since entropy dropping with will likely drive efficient convection enforcing vertically homogeneous entropy.

However, viscously heated disks accreting at low should be optically thin to their own emission (Narayan & Popham, 1993). If the specific viscous dissipation rate increases with height (which is expected, for example, if heating is produced by the MRI (Miller & Stone, 2000; Hirose et al., 2006)), then gas temperature will increase with at all altitudes. This can potentially provide the conditions favorable for the radial inflow at all heights, as we showed in this work. Thus, optically thin, viscously heated accretion disks present the best setting for suppressing the radial outflow at all altitudes.

## 7. Summary

We explored meridional circulation in accretion disks with shear viscosity and varied thermodynamics. While previous studies of this problem focused on the vertically isothermal disks, finding radial gas outflow at the midplane and inflow at high latitudes, we demonstrate that different assumptions about thermal stratification can change this pattern. We show that the direction of the flow at the midplane is intimately connected to the behavior of the vertical stress . Vertical density profiles steeply falling off with height induce significant vertical shear, which can make the gas at the disk midplane to flow inward. We derive analytical criterion relating the direction of the midplane flow to the thermal stratification in the disk and show that the radial inflow at all heights (without change of sign of ) naturally sets in disks with steeply growing vertical profiles of entropy and temperature. Such conditions can be naturally realized in optically thin disks heated primarily by viscous dissipation. Although our findings rely on the assumption of shear viscosity, we also comment on other mechanisms of angular momentum transport. Our results can also be used for code testing and designing simulations with the prescribed pattern of meridional circulation.

We are grateful to Willy Kley for useful discussions. A.A.P. is supported by Porter Ogden Jacobus Fellowship, awarded by the graduate school of Princeton University. Financial support for this study has been provided by the NSF via grant AST-1515763, NASA via grants 14-ATP14-0059, and The Ambrose Monell Foundation. Simulations presented in this article used computational resources supported by the PICSciE-OIT TIGRESS High Performance Computing Center and by NSF through an XSEDE computational time allocation TG-AST160008 on TACC Stampede and Ranch.

## Appendix A Derivation of uR.

Plugging the expressions (5) into the equation (4) and using the disk symmetry with respect to , one finds that the radial velocity at the disk midplane is

 uR(R,0) = 2νRΩK[∂2Ω∂R2+∂Ω∂R×∂ln(ρνR3)∂R∣∣∣z→0+∂2Ω∂z2∣∣∣z→0]. (A1)

Here we approximated , which is accurate at the level. To the same degree of accuracy we can replace with in all terms with the radial derivatives of . However, the last term in the brackets must be treated more carefully. Using the relation (3) one finds

 ∂2Ω∂z2∣∣∣z→0 ≈ ΩK2R2[−3+1Ω2K∂2∂z2(1ρ∂P∂lnR)∣∣∣z→0] (A2) ≈ −ΩK2R2[δρ+PΩ2Kρ2∂lnP∂lnR∂2ρ∂z2∣∣∣z→0], (A3)

where is defined by equation (7). In going from (A2) to (A3) we used equation (2) and disk symmetry with respect to its midplane. Plugging result (A3) into the equation (A1) one arrives at the expression (6).

Further progress involves the knowledge of the EOS of the disk fluid. Plugging the anzatz (8) into the equation (2), taking a derivative of both sides with respect to , and using the symmetry property at one finds

 PΩ2Kρ2∂2ρ∂z2∣∣∣z→0=−1γ(1+PΩ2Kρ∂2s∂z2∣∣∣z→0). (A4)

## References

• Bai & Stone (2013) Bai, X.-N., & Stone, J. M. 2013, ApJ, 767, 30
• Balbus (2003) Balbus, S. A. 2003, ARA&A, 41, 555
• Belyaev et al. (2012) Belyaev, M. A., Rafikov, R. R., & Stone, J. M. 2012, ApJ, 760, 22
• Belyaev et al. (2013a) —. 2013a, ApJ, 770, 67
• Belyaev et al. (2013b) —. 2013b, ApJ, 770, 68
• Béthune et al. (2016) Béthune, W., Lesur, G., & Ferreira, J. 2016, ArXiv e-prints, arXiv:1612.00883
• Brownlee et al. (2006) Brownlee, D., Tsou, P., Aléon, J., et al. 2006, Science, 314, 1711
• Canup (2004) Canup, R. M. 2004, ARA&A, 42, 441
• Chiang & Goldreich (1997) Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368
• Ciesla (2007) Ciesla, F. J. 2007, Science, 318, 613
• Ciesla (2009) —. 2009, Icarus, 200, 655
• Flock et al. (2011) Flock, M., Dzyurkevich, N., Klahr, H., Turner, N. J., & Henning, T. 2011, ApJ, 735, 122
• Fromang et al. (2011) Fromang, S., Lyra, W., & Masset, F. 2011, A&A, 534, A107
• Gardiner & Stone (2008) Gardiner, T. A., & Stone, J. M. 2008, Journal of Computational Physics, 227, 4123
• Gehrz et al. (1998) Gehrz, R. D., Truran, J. W., Williams, R. E., & Starrfield, S. 1998, PASP, 110, 3
• Hartmann & Davis (1975) Hartmann, W. K., & Davis, D. R. 1975, Icarus, 24, 504
• Hirose et al. (2006) Hirose, S., Krolik, J. H., & Stone, J. M. 2006, ApJ, 640, 901
• Hughes & Armitage (2010) Hughes, A. L. H., & Armitage, P. J. 2010, ApJ, 719, 1633
• Illarionov & Sunyaev (1975) Illarionov, A. F., & Sunyaev, R. A. 1975, A&A, 39, 185
• Jacquet (2013) Jacquet, E. 2013, A&A, 551, A75
• Kley & Lin (1992) Kley, W., & Lin, D. N. C. 1992, ApJ, 397, 600
• Kluzniak & Kita (2000) Kluzniak, W., & Kita, D. 2000, ArXiv Astrophysics e-prints, astro-ph/0006266
• Lynden-Bell & Pringle (1974) Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603
• Miller & Stone (2000) Miller, K. A., & Stone, J. M. 2000, ApJ, 534, 398
• Narayan & Popham (1993) Narayan, R., & Popham, R. 1993, Nature, 362, 820
• Pahlevan & Stevenson (2007) Pahlevan, K., & Stevenson, D. J. 2007, Earth and Planetary Science Letters, 262, 438
• Papaloizou & Lin (1995) Papaloizou, J. C. B., & Lin, D. N. C. 1995, ARA&A, 33, 505
• Popham & Narayan (1995) Popham, R., & Narayan, R. 1995, ApJ, 442, 337
• Popham et al. (1993) Popham, R., Narayan, R., Hartmann, L., & Kenyon, S. 1993, ApJ, 415, L127
• Pringle (1991) Pringle, J. E. 1991, MNRAS, 248, 754
• Rafikov (2013) Rafikov, R. R. 2013, ApJ, 774, 144
• Rafikov (2016a) —. 2016a, ApJ, 827, 111
• Rafikov (2016b) —. 2016b, ApJ, 830, 7
• Regev & Gitelman (2002) Regev, O., & Gitelman, L. 2002, A&A, 396, 623
• Rozyczka et al. (1994) Rozyczka, M., Bodenheimer, P., & Bell, K. R. 1994, ApJ, 423, 736
• Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
• Siemiginowska (1988) Siemiginowska, A. 1988, Acta Astronomica, 38, 21
• Stoll & Kley (2014) Stoll, M. H. R., & Kley, W. 2014, A&A, 572, A77
• Stoll & Kley (2016) —. 2016, A&A, 594, A57
• Stoll et al. (2017) Stoll, M. H. R., Picogna, G., & Kley, W. 2017, ArXiv e-prints, arXiv:1702.00334
• Suzuki & Inutsuka (2014) Suzuki, T. K., & Inutsuka, S.-i. 2014, ApJ, 784, 121
• Syunyaev & Shakura (1977) Syunyaev, R. A., & Shakura, N. I. 1977, Soviet Astronomy Letters, 3, 138
• Takeuchi & Lin (2002) Takeuchi, T., & Lin, D. N. C. 2002, ApJ, 581, 1344
• Truran & Livio (1986) Truran, J. W., & Livio, M. 1986, ApJ, 308, 721
• Urpin (2003) Urpin, V. 2003, A&A, 404, 397
• Urpin & Brandenburg (1998) Urpin, V., & Brandenburg, A. 1998, MNRAS, 294, 399
• Urpin (1984) Urpin, V. A. 1984, AZh, 61, 84
• Wiechert et al. (2001) Wiechert, U., Halliday, A. N., Lee, D.-C., et al. 2001, Science, 294, 345
• Zhu & Stone (2017) Zhu, Z., & Stone, J. M. 2017, ArXiv e-prints, arXiv:1701.04627
• Zolensky et al. (2006) Zolensky, M. E., Zega, T. J., Yano, H., et al. 2006, Science, 314, 1735
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters