Radial Transport and Meridional Circulation in Accretion Disks

TL;DR

This paper investigates how thermal stratification influences radial transport and meridional circulation patterns in 3D accretion disks, revealing conditions that lead to inward or outward flows relevant for astrophysical processes.

Contribution

It demonstrates analytically and numerically that steeply increasing entropy and temperature with height can cause inward radial flow at all altitudes in accretion disks.

Findings

Steep thermal gradients induce inward flow at all altitudes.

Disks with typical conditions exhibit midplane outflow and upper inflow.

Results inform simulation design for prescribed circulation patterns.

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.

Equations

$${\Omega }^{2}R-\frac{1}{\rho }\frac{\partial P}{\partial R}$$

$$\frac{1}{\rho }\frac{\partial P}{\partial z}$$

$$\Omega (R,z)\approx {\Omega }_{\mathrm{K}}\left(1-\frac{3}{4}\frac{z^2}{R^2}+\frac{1}{2{\Omega }_{\mathrm{K}}^2{R}^2}\frac{P}{\rho }\frac{\partial \ln P}{\partial \ln R}\right).$$

$$u_R(R,z)$$

$$T_{R\varphi }=-\rho \nu R\frac{\partial \Omega }{\partial R},T_{z\varphi }=-\rho \nu R\frac{\partial \Omega }{\partial z},$$

$$u_R(R,0)$$

$${\delta }_{\rho }\equiv \partial \ln {\rho }_0/\partial \ln R,$$

$$P=e^s{\rho }^{\gamma },s(R,Z)=c_V^{-1}S(R,z),$$

$$s^{\mathrm{iso}}(R,z)=s_0(R)+\frac{\gamma -1}{2}{\left[\frac{z}{H_0(R)}\right]}^2,$$

$$\displaystyle\frac{\delta_{\rho}+\delta_{T}}{\gamma}\left(1+\frac{P}{\Omega_{\rm K}^{2}\rho}\frac{\partial^{2}s}{\partial z^{2}}\Bigg{|}_{z\to 0}\right),$$

$${\delta }_T\equiv \partial \ln T_0/\partial \ln R,$$

$$F_J=-\underset{-\infty }{\overset{\infty }{\int }}2\pi R^3\rho \nu \frac{\partial \Omega }{\partial R}dz=3\pi l\underset{-\infty }{\overset{\infty }{\int }}\rho \nu dz,$$

$$\dot{M}(R)=-2\pi R\underset{-\infty }{\overset{\infty }{\int }}\rho u_{R}dz=\frac{\partial F_J}{\partial l},$$

$$S(R,z)=S\left(R,\frac{z}{H_0(R)}\right),H_0(R)\equiv H(R,0),$$

$$F_J\propto {\rho }_0{\nu }_0{H}_{0}l\propto {\rho }_0{\nu }_0{T}_0^{1/2}{R}^2,$$

$$u_R(R,0)$$

$$\displaystyle s^{\prime\prime}_{zz}\equiv H_{0}^{2}\frac{\partial^{2}s}{\partial z^{2}}\Bigg{|}_{z\to 0}$$

$${\delta }_F\equiv \partial \ln F_J/\partial \ln R.$$

$$s_{zz}^{\prime \prime }$$

$$s_{zz}^{\prime \prime }>-\frac{\gamma }{{\delta }_{\rho }+{\delta }_T}\left(3+\frac{3\gamma +2}{2\gamma }{\delta }_T+\frac{1-\gamma }{\gamma }{\delta }_{\rho }\right).$$

$${\delta }_{\rho }=-3-\frac{3}{2}{\delta }_T.$$

$$s_{zz}^{\prime \prime }>\frac{6(2\gamma -1)+(6\gamma -1){\delta }_T}{6+{\delta }_T}.$$

$${\delta }_{\rho }>\frac{6\gamma +(3\gamma +2){\delta }_T}{2(\gamma -1)}.$$

$$\Lambda =-\rho \frac{T-T_0}{\tau },$$

$$s(R,z)=s_0(R)+\frac{s_{zz}^{\prime \prime }}{2}{\left[\frac{z}{H_0(R)}\right]}^2,$$

$${\delta }_{\rho }=-4-\frac{3}{2}{\delta }_T.$$

$$s_{zz}^{\prime \prime }>\frac{2(7\gamma -4)+(6\gamma -1){\delta }_T}{8+{\delta }_T}.$$

$$u_R(R,0)$$

$$\displaystyle\frac{\partial^{2}\Omega}{\partial z^{2}}\Bigg{|}_{z\to 0}$$

$$\approx$$

$$\displaystyle\frac{P}{\Omega_{\rm K}^{2}\rho^{2}}\frac{\partial^{2}\rho}{\partial z^{2}}\Bigg{|}_{z\to 0}=-\frac{1}{\gamma}\left(1+\frac{P}{\Omega_{\rm K}^{2}\rho}\frac{\partial^{2}s}{\partial z^{2}}\Bigg{|}_{z\to 0}\right).$$

Full text

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 $\Sigma$ (Shakura & Sunyaev, 1973; Lynden-Bell & Pringle, 1974). In this approach the steady accretion disks with the radially constant mass accretion rate $\dot{M}$ 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 $u_{R}$ is always negative, $\langle u_{R}\rangle_{\rho}=\Sigma^{-1}\int\rho u_{R}dz<0$.

Situation becomes more complicated once one considers the full three-dimensional disk structure, in particular, the profile of the radial velocity $u_{R}$ as a function of the vertical coordinate $z$. 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 $z=0$, implying a radial ouflow near the disk midplane. The magnitude of positive $u_{R}$ steadily decreases with the vertical distance from the midplane, and $u_{R}$ ultimately changes sign at some altitude of order the disk scale height $H$. As a result, a meridional circulation pattern sets in the poloidal plane of the disk. Radial gas inflow at high $z$ 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 $\langle u_{R}\rangle_{\rho}=\Sigma^{-1}\int\rho u_{R}dz$ 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 $\alpha$-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 ($>1\mu$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 $\alpha$-viscosity (Shakura & Sunyaev, 1973). All of them invariably find $u_{R}(R,0)>0$ 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 $T_{ij}$ is anisotropic. The anisotropy of $T_{ij}$ 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 $z$, 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 $M_{\star}$; we work in cylindrical coordinates $(R,z)$. We are interested in the spatial structure of the radial velocity $u_{R}(R,z)$, which is non zero because of the internal stress $T_{ij}$ 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

[TABLE]

where $\rho$ and $P$ are the gas density and pressure, $\Omega$ is the angular frequency, and $\Phi(R,z)=-GM_{\star}(R^{2}+z^{2})^{-1/2}$ 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 $\Omega_{\rm K}\equiv(GM_{\star}/R^{3})^{1/2}$:

[TABLE]

This approximate relation follows from our assumption of a geometrically thin disk, for which $H/R\ll 1$, $z/R\ll 1$, where $H(R,z)\equiv c_{s}/\Omega_{\rm K}$ is the local disk scaleheight, and $c_{s}^{2}(R,z)\equiv P/\rho=k_{B}T(R,z)/\mu$ 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 $\Omega$ is independent of $z$ for a barotropic equation of state. $R$ and $z$. The deviations from purely Keplerian rotation in equation (3) are at the $(z/R)^{2}\ll 1$ (second term in parentheses) and $(H/R)^{2}\ll 1$ (third term) level.

The $\phi$-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 $u_{R}$ of the disk fluid (Jacquet, 2013):

[TABLE]

Here $l\equiv\Omega R^{2}$ is the specific angular momentum of the disk fluid, while $T_{R\phi}$ and $T_{z\phi}$ represent the $R-\phi$ and $z-\phi$ components of the internal stress tensor. This expression is accurate at the $(z/R)^{2},(H/R)^{2}\ll 1$ 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 $\rho$ and integrating it over $z$, one finds that $\langle u_{R}\rangle_{\rho}$ is determined by the $T_{R\phi}$ 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

[TABLE]

where $\nu$ 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 $z=0$. For that reason, our primary goal will be to determine the factors that affect the midplane value of the radial velocity, $u_{R}(R,0)$. To that effect, in our analysis we will often take the limit $z\to 0$ 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 $z$, e.g. $\partial\rho/\partial z$, $\partial\Omega/\partial z$, etc., go to zero as $z\to 0$.

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

[TABLE]

where

[TABLE]

and subscript [math] implies the value of a particular variable at the midplane, e.g. $\rho_{0}(R)\equiv\rho(R,0)$, $P_{0}(R)\equiv P(R,0)$, etc. Definition (7) implies that $\delta_{\rho}$ is closely related to the density power law exponent $p$ used in previous studies (Takeuchi & Lin, 2002; Fromang et al., 2011). The two coincide if $\rho_{0}(R)\propto R^{p}$, 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 $T_{R\phi}$, while the last two terms are due to the vertical stress $T_{z\phi}$. Since the midplane pressure $P_{0}$ is a decreasing function of $R$, it is clear that a vertical density profile steeply declining with $z$ (i.e. high value of $-\partial^{2}\rho/\partial z^{2}$) should be favorable for suppressing the midplane outflow in the disk and driving gas inflow ($u_{R}<0$) 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 $O(\alpha^{2})$ relative contribution to the expression (3) for the angular frequency. Their inclusion would, in turn, lead to the emergence of an additional $O(\alpha^{3})$ contribution in the expression (6) for $u_{R}$. Given the expectation of $\alpha\ll 1$ 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 $\alpha\gtrsim 0.7$, although Kley & Lin (1992) suggest a lower value of critical $\alpha\approx 0.06$ for this transition.. This expectation is subsequently confirmed in §5 by the good agreement between our analytical theory (which explicitly neglects viscous contributions to $\Omega(R,z)$) 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

[TABLE]

where $\gamma$ is the polytropic index and $s(R,z)$ is the scaled (dimensionless) gas entropy $S(R,z)$; $c_{V}\equiv(\gamma-1)^{-1}k_{B}/\mu$ is the specific heat capacity. Note that in this work $S(R,z)$ 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 $P/P_{0}=\rho/\rho_{0}=\exp(-z^{2}/2H_{0}^{2})$, so that

[TABLE]

where $s_{0}=\ln\left(P_{0}/\rho_{0}^{\gamma}\right)$. A different limit of the locally isentropic disks considered later in §3.2.1 has $s(R,z)=s(R)$.

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

[TABLE]

where

[TABLE]

and $T_{0}(R)\equiv T(R,0)$ is the midplane value of the gas temperature. Definition of $\delta_{T}$ makes it closely related to the temperature power law index $q$ used in other studies of the meridional circulation in accretion disks (Takeuchi & Lin, 2002; Fromang et al., 2011), which assume $T_{0}(R)\propto R^{q}$.

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 $F_{J}$, which is equal to the viscous torque exerted by the inner disk on the outer disk. The concept of $F_{J}$ 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, $F_{J}$ can be calculated by multiplying $T_{R\phi}$ in equation (5) by $2\pi R^{2}$ and integrating over $z$. Using equation (5) this implies, to $(H/R)^{2}$ accuracy, that

[TABLE]

where $l\equiv\Omega_{K}R^{2}$ and we have assumed a Keplerian rotation profile. Also, multiplying equation (4) by $2\pi R\rho$ and integrating over $z$ one finds that the mass accretion rate $\dot{M}$ (defined to be positive for inflow) is related to $F_{J}$ as (Rafikov, 2013)

[TABLE]

where the derivative is with respect to the specific angular momentum, so that $\partial/\partial l=2(\Omega_{\rm K}R)^{-1}\partial/\partial R$.

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

[TABLE]

where $H_{0}=c_{s}(R,0)/\Omega_{\rm K}$ 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 $\rho$ and $P$ also depend on $z$ only in combination $z/H_{0}$; this makes it natural to expect that $\nu=\nu(R,z/H_{0})$ as well. Then equation (12) implies that

[TABLE]

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

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

[TABLE]

where we introduce the shorthand notation

[TABLE]

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

[TABLE]

Once again, in equation (16) the first two terms arise from $T_{R\phi}$, while the terms in the second line originate from $T_{z\phi}$.

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 $s^{\prime\prime}_{zz}$. 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 $z$, including the midplane. Setting $u_{R}(R,0)<0$ one finds a necessary criterion for this to be the case:

[TABLE]

where we assumed that $P_{0}(R)$ decreases with $R$ so that $\delta_{\rho}+\delta_{T}<0$. 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 $\delta_{F}$, $\delta_{\rho}$, and $\delta_{T}$ — 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 $\nu$ depends on $z$ only in combination $z/H_{0}$).

3.1. Inflow criterion for a standard constant $\dot{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 ($\dot{M}=$const), with no torque applied at its center. In this case, integrating equation (13) one finds $F_{J}=\dot{M}l$ (Rafikov, 2013), so that $\delta_{F}=1/2$. This transforms inequality (19) into the following inflow criterion:

[TABLE]

Note that the constant $\dot{M}$ assumption does not constrain $\delta_{\rho}$ or $\delta_{T}$ because of the freedom in choosing the radial variation of the viscosity $\nu_{0}$. This constraint is illustrated in Figure 1 for two values of $\gamma$ ($5/3$ and $7/5$) and three values of $\delta_{\rho}$. One can see that steeper decay of $\rho_{0}$ with $R$ (more negative $\delta_{\rho}$) makes it easier to achieve inflow at the disk midplane, i.e. requires less extreme values of $s^{\prime\prime}_{zz}$.

We can now go one step further and adopt a particular viscosity ansatz, namely the $\alpha$-model of Shakura & Sunyaev (1973), in which $\nu=\alpha c_{s}^{2}/\Omega_{\rm K}$. If we additionally assume that the effective viscosity parameter $\alpha$ is independent of $R$, then locally $\nu_{0}\propto R^{\delta_{T}+3/2}$. Plugging this and $F_{J}=\dot{M}l$ into equation (15) one obtains that

[TABLE]

This results in yet another version of the inflow criterion for $\alpha$ disks, which depends only on the radial temperature profile (i.e. $\delta_{T}$):

[TABLE]

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 $\dot{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 isentropic disk, in which entropy (and $s$) does not depend on $z$. 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 $\rho_{0}(R)$ and $T_{0}(R)$ in the form

[TABLE]

For protoplanetary disks with $\gamma=7/2$ this becomes $\delta_{\rho}>(42+31\delta_{T})/4$. In hotter accretion disks with $\gamma=5/3$ the constraint is $\delta_{\rho}>3(10+7\delta_{T})/4$.

It is clear that even if $T_{0}(R)$ decays with radius as rapidly as $R^{-1}$ (i.e. $\delta_{T}=-1$) 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 ($\partial T(R,z)/\partial z=0$) turns the left hand side of the inflow criterion into $\gamma-1$, see equation (9). Then it is easy to see from the inequality (20) that the inflow at all heights is possible only if $\delta_{T}<-6/5$. This result is completely independent of either $\delta_{\rho}$ or $\gamma$, thus it should hold for an arbitrary radial profile of the midplane density $\rho_{0}(R)$.

Finding conditions in which a disk might have $T_{0}(R)$ decaying with radius faster than $R^{-6/5}$ is not easy. Based on this we can conclude that the locally isothermal constant $\dot{M}$ 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 $z=0$ requires the disk to have temperature increasing with height, starting from the midplane. Such disks naturally have entropy rising more steeply with $z$ than in the isothermal case and may satisfy the conditions (20) and (22) without making unrealistic assumptions about the behavior of $\rho_{0}(R)$ and $T_{0}(R)$.

3.3. Inflow criterion for $F_{J}=$ const disk

Standard constant-$\dot{M}$ 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 $F_{J}(R)=\dot{M}l+F_{J}(0)$ (Rafikov, 2013, 2016a), where $F_{J}(0)$ represents the torque applied at the disk center (which is naturally absent in a standard constant-$\dot{M}$ 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 $F_{J}(R)=$const, implying that $\delta_{F}=0$.

Equation (16) demonstrates that reaching inflow at the midplane of such a disk is more difficult than in the standard constant-$\dot{M}$ 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 $s^{\prime\prime}_{zz}$ is needed to guarantee inflow at all heights (i.e. $u_{R}(R,0)<0$) in the constant-$F_{J}$ disk, as compared to the standard constant $\dot{M}$ 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 $z=0$ in a locally isothermal, constant-$F_{J}$ disk would require $\delta_{T}<-9/5$. 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 $r,\theta,\phi$ coordinates, using uniform grid in $\log r$ and constraining ourselves to small range in $\phi$. We use the domain $[1,10]\times[\pi/2-0.5,\pi/2+0.5]\times[0,0.1]$ in these coordinates with numerical resolution of $1024\times 512\times 16$ cells. We verified that simulations in the full $\phi$ range $[0,2\pi]$ lead to same results. We also verified that our results are converged with regard to numerical resolution.

For $r$ and $\theta$ boundaries, we use the so-called ”do-nothing” boundary condition, i.e. all fluid variables are fixed at their initial values. In $\phi$-direction, we impose a periodic boundary condition. The Mach number of the orbital flow at the inner radius is $10$ in all our simulations. The simulation time is always measured in units of $\Omega_{K}^{-1}$, where $\Omega_{K}$ 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:

[TABLE]

where $T$ is the gas temperature, $T_{0}(R,z)$ is the initial temperature profile described in §4.1, $\tau=0.5\Omega_{K}^{-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 $\rho_{0}(R)\propto R^{\delta_{\rho}}$, $T_{0}(R)\propto R^{\delta_{T}}$. Our simulations adopt a setup typical for standard constant $\dot{M}$ disks (i.e. we do not simulate disks considered in §3.3), making different assumptions about the radial profile of $\alpha$. Close to the mid-plane, we consider the following entropy profile

[TABLE]

such that the definition (17) holds true; also $s_{0}(R)=\ln\left(P_{0}/\rho_{0}^{\gamma}\right)$. 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 $|z|>2H_{0}$.

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 $\delta_{T}=-1$, constant $\dot{M}$ and two values of $s^{\prime\prime}_{zz}$ in Figure 2. It shows that for high values of the entropy derivative $s^{\prime\prime}_{zz}$ 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 $s^{\prime\prime}_{zz}<\gamma-1$ result in $T$ dropping with $z$ (by design, at high altitudes our $T$ 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 $\alpha$ profile with $\alpha=0.01$, so that equations (21) and (22) hold. In a second set, we assume that $\alpha\propto R$, with $\alpha=0.01$ at $R=1$. In this case, the constant $\dot{M}$ assumption leads to a different relation between the power-law indices of density and temperature profiles:

[TABLE]

Plugging this $\delta_{\rho}$ in equation (20), we obtain a new criterion for the inflow at all $z$, which replaces equation (22):

[TABLE]

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

For both prescriptions for $\alpha$-viscosity, we perform simulations for several sets of different values of $s^{\prime\prime}_{zz}$ and $\delta_{T}$.

5. Simulation results

We start by presenting the results of simulations for $\alpha(R)=$const case. In Figure 3 we show the steady state distribution of the radial velocity $u_{R}$ in the $R-z$ plane for $T(R)$ profile with $\delta_{T}=-1$ and two values of $s^{\prime\prime}_{zz}=0.5$ and 1.5. The two disks have identical radial profiles of all fluid variables in the midplane, but different vertical profiles of $\rho$ and $T$. Clearly, this leads to an important difference: while $s^{\prime\prime}_{zz}=0.5$ case shows outflow of gas in the midplane and inflow at high altitudes, the $s^{\prime\prime}_{zz}=1.5$ run shows inflow at all altitudes, in agreement with analytical expectations (as shown in Figure 2, in this case $T$ increases with height).

In Figure 4a we show meridional profiles of the radial velocity, normalized to $\nu/R$, at different radii for $s^{\prime\prime}_{zz}=0.75$. Given our thermodynamic prescription (25), we expect and find to hold in simulations the following parabolic profile of the radial velocity: $u_{R}(R,z)\propto\left[1+C\times(z/H_{0})^{2}\right]$, $C<0$ (c.f. Jacquet 2013). In the case of $\delta_{T}=-1$, when $H_{0}=R$ this behavior results in the meridional profile of $u_{R}$ normalized by $\nu/R$, which is independent of $R$. 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 $u_{R}$ profile to a steady state. One can see that convergence happens within several tens $\Omega_{K}^{-1}$, 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 $u_{R}$ (at fixed $R$ and $\delta_{T}$) vary as we change the value of $s^{\prime\prime}_{zz}$. One can see that increasing the rate at which gas 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^{\prime\prime}_{zz}$ 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 $\delta_{T}=-1.5,-1,-0.5$. In each run, we measure the radial component of the gas velocity in the midplane $u_{R}(R,0)$ and record its sign. The outcome of these runs is shown in Figure 5a, where we mark solutions with $u_{R}(R,0)<0$ as blue stars, with $u_{R}(R,0)>0$ as red, and simulations with $u_{R}(R,0)\approx 0$ as yellow. Analytical inflow criterion (22) is shown as a solid line. We find that for each inspected value of $\delta_{T}$, there is a critical value of the vertical entropy gradient, $s^{\prime\prime}_{zz}$, above which the disk exhibits the inflow at all $z$. 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 $\alpha\propto R$. Their results, together with the analytical criterion (27), are summarized in Figure 5b. As in the case of $\alpha=$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 $T_{r\phi}$ to $u_{R}(R,0)$ positive, promoting a midplane outflow (Takeuchi & Lin, 2002). On the other hand, at high altitudes the density gradient with respect to $R$ changes sign, so that $T_{r\phi}$ drives inflow at high $z$. This leads to the characteristic parabolic shape of the vertical $u_{R}$ 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 $\partial^{2}\rho/\partial z^{2}<0$ promotes inflow in non-pathological cases when $P_{0}$ decreases with $R$. It arises because the gas at higher $z$ 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 $\rho$ on $\Omega$: density slowly decaying with $z$ results in lower $|\partial^{2}\Omega/\partial z^{2}|$, 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 $z$ 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 $T_{r\phi}$ to $u_{R}(R,0)$. 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 $T_{r\phi}$ behavior for the MRI, which determines $\langle u_{R}\rangle_{\rho}$ and, thus, sets the accretion rate $\dot{M}$. As found by many numerical studies for locally isothermal disks, $T_{r\phi}$ is nearly constant at low altitudes and falls off sharply in the corona, where plasma $\beta$ 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 $T_{z\phi}$, 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 $u_{R}$ 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 $T_{r\phi}$ nor $T_{z\phi}$ was assumed444We note, however, that the vertical profile of $T_{z\phi}$ 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 $u_{R}(R,z)$ different from Fromang et al. (2011): inflow at the midplane, and outflow at high $z$, i.e. $u_{R}(R,z)$ 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 $\Omega$ 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 $u_{R}(R,z)$ and the vertical variation of the stress tensor (including $T_{z\phi}$!), 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 $T_{r\phi}$ and $T_{z\phi}$ are characterized by the same value of $\nu$), 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: $T_{z\phi}$ was found to exceed $T_{R\phi}$ by a factor of several hundred. Moreover, the results of Stoll et al. (2017) suggest that near the midplane the vertical variation of $T_{z\phi}$ can be reasonably well approximated by the shear stress prescription (5). The radial stress $T_{R\phi}$ certainly does not follow the scaling $T_{R\phi}\propto\rho$ postulated in equations (5), as Stoll & Kley (2014) and Stoll et al. (2017) find $T_{R\phi}$ to increase with $z$ near the midplane, switching to a decay only at high altitudes. However, given the negligible role played by $T_{R\phi}$ 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 $T_{R\phi}$. We also set $s^{\prime\prime}_{zz}=\gamma-1$ 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 $T_{z\phi}$, naturally result in $u(R,0)=\nu_{z}\delta_{T}/R$, where $\nu_{z}$ is the value of the kinematic viscosity coefficient characterizing the vertical behavior of the $T_{z\phi}$ as found in Kley et al. Given that the disk temperature decreases with $R$ (i.e. $\delta_{T}<0$), one finds that $u(R,0)<0$ (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 $|z|$. 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 $z$ (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 $\dot{M}$ are likely to be optically thick; their $T(R,z)$ 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 ($s^{\prime\prime}_{zz}=0$) may feature the fastest radial outflow at $z=0$, since entropy dropping with $z$ will likely drive efficient convection enforcing vertically homogeneous entropy.

However, viscously heated disks accreting at low $\dot{M}$ 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 $z$ 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 $T_{z\phi}$. 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 $u_{R}$) 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 $u_{R}$.

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

[TABLE]

Here we approximated $l=\Omega_{\rm K}R^{2}$, which is accurate at the $(H/R)^{2}$ level. To the same degree of accuracy we can replace $\Omega$ with $\Omega_{\rm K}$ in all terms with the radial derivatives of $\Omega$. However, the last term in the brackets must be treated more carefully. Using the relation (3) one finds

[TABLE]

where $\delta_{\rho}$ 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 $z$, and using the symmetry property at $z=0$ one finds

[TABLE]