Anisotropic hydrodynamic turbulence in accretion disks

TL;DR

This paper investigates the anisotropic nature of turbulence caused by the vertical shear instability in accretion disks, revealing that vertical stresses dominate and influence meridional circulation, with implications for particle transport.

Contribution

It demonstrates that VSI turbulence is anisotropic, with vertical stresses exceeding radial stresses, and shows this explains the observed meridional circulation pattern in disks.

Findings

Vertical stress is significantly stronger than radial stress.

Good agreement between viscous simulations and VSI turbulence results.

Anisotropy impacts the transport of small particles in disks.

Abstract

Recently, the vertical shear instability (VSI) has become an attractive purely hydrodynamic candidate for the anomalous angular momentum transport required for weakly ionized accretion disks. In direct three-dimensional numerical simulations of VSI turbulence in disks, a meridional circulation pattern was observed that is opposite to the usual viscous flow behavior. Here, we investigate whether this feature can possibly be explained by an anisotropy of the VSI turbulence. Using three-dimensional hydrodynamical simulations, we calculate the turbulent Reynolds stresses relevant for angular momentum transport for a representative section of a disk. We find that the vertical stress is significantly stronger than the radial stress. Using our results in viscous disk simulations with different viscosity coefficients for the radial and vertical direction, we find good agreement with the VSI turbulence for the stresses and meridional flow; this provides additional evidence for the anisotropy. The results are important with respect to the transport of small embedded particles in disks.

Tables

Table 1. Table 1: Model parameter: Domain size and grid resolution

Parameter	$\alpha$-model	VSI
Radial range [$5.2\text{au}$]	$0.4$ - $2.5$	$0.4$ - $2.5$
Vertical range [H]	$\pm 4$	$\pm 5$
Phi range [rad]	2D	$0$ - $2\pi$
Radial grid size	$600$	$600$
Theta grid size	$121$	$128$
Phi grid size	2D	$1024$

Equations

$$\rho (r,z)={\rho }_0{\left(\frac{r}{r_0}\right)}^p\exp \left[\frac{GM}{c_{\mathrm{s}}^2}\left(\frac{1}{R}-\frac{1}{r}\right)\right],$$

$$T(r)=T_0{\left(\frac{r}{r_0}\right)}^q,$$

$$\Omega (r,z)={\Omega }_K{\left[(p+q){\left(\frac{H}{r}\right)}^2+(1+q)-\frac{qr}{\sqrt{r^2+z^2}}\right]}^{\frac{1}{2}},$$

$$R_{i,\varphi }(z)=<\rho u_i\delta u_{\varphi }{>}_{t,\varphi ,r},$$

$$\frac{\partial \rho r^2\Omega }{\partial t}+\nabla \cdot (\text{@vec}u\rho r^2\Omega )=-\frac{\partial P}{\partial \varphi }+\nabla \cdot (r\text{@vec}{S}_{\varphi }).$$

$$\text{@vec}{S}_{\varphi }=(S_{r\varphi },S_{z\varphi })=\left(\rho {\nu }_{r}r\frac{\partial \Omega }{\partial r},\rho {\nu }_{z}r\frac{\partial \Omega }{\partial z}\right).$$

$$u_r\propto \frac{1}{\rho r^2\Omega }\left(\frac{1}{r}\frac{\partial (r^2{S}_{r\varphi })}{\partial r}+\frac{\partial r{S}_{z\varphi }}{\partial z}\right).$$

$$\Omega (r,z)={\Omega }_K(r)\left(1+\frac{q}{4}\frac{z^2}{r^2}+\frac{q+p}{2}{h}^2\right),$$

$$\frac{\partial \Omega }{\partial z}={\Omega }_K\frac{q}{2}\frac{z}{r^2}$$

$$S_{z\varphi }=q\rho {\nu }_z\frac{{\Omega }_K}{2}\frac{z}{r}\text{mbox}and{S}_{r\varphi }=-\rho {\nu }_r\frac{3}{2}{\Omega }_K.$$

$$S_{z\varphi }=-\frac{q}{3}{S}_{r\varphi }\frac{{\nu }_z}{{\nu }_r}\frac{z}{r}.$$

$$u_r(z=0)$$

$$u_r(z=0)\propto \left[\frac{3}{2}{\alpha }_r-{\alpha }_z\right].$$

$$r=R\sin \theta \text{and}z=R\cos \theta ,$$

$$S_{R\varphi }$$

$$S_{\theta \varphi }$$

$$\frac{\partial \Omega }{\partial r}$$

$$\frac{\partial \Omega }{\partial z}$$

$$S_{R\varphi }$$

$$S_{\theta \varphi }$$

Full text

1 Introduction

The exact origin of the driving mechanism of accretion disks is still not fully understood. To accrete matter onto the central object, the matter needs to lose its angular momentum, and because molecular viscosity is by many orders of magnitudes too small to facilitate the required angular transport, it has been suggested that disks are driven by turbulence. The discovery of a linear magneto-rotational instability (MRI) for rotating flows with a negative angular velocity gradient has led to the suggestion that accretion disks are driven by magnetohydrodynamical (MHD) turbulence (Balbus & Hawley, 1998). However, the MRI only works efficiently for well-ionized media, for example, in disks around compact objects, but for lower ionization levels the non-ideal MHD effects become stronger and the operability of the MRI questionable. In particular for the cool and low-ionized regions, so-called dead zones with very low turbulent activity have been predicted (Gammie, 1996).

Hence, other alternatives are sought for. In the past years the purely hydrodynamical vertical shear instability (VSI) has attracted some attention because the only requirement is a vertical shear in the angular velocity profile, which is in fact a natural consequence of a radial temperature gradient in the disk, for example, induced by irradiation from the central object. Through numerical simulations and linear analysis it has been shown that the VSI operates efficiently for vertically isothermal disks (Nelson et al., 2013) as well as for fully radiative disks that include stellar irradiation (Stoll & Kley, 2014). This makes the VSI a promising candidate to bring at least some life back into the dead zones. For many astrophysical applications it is useful to parameterize the turbulence and describe the angular momentum transport by an effective viscous prescription, for example, the well-known ansatz by Shakura & Sunyaev (1973) where the uncertainties of the turbulence are summarized in one constant parameter $\alpha$ (Balbus & Papaloizou, 1999). Such an approach using one parameter is applicable for an isotropic turbulence, and useful when the interest is in the overall radial evolution of the disk. To analyze the internal flow field of the disk, which is important for the motion of embedded small particles, it is necessary to take possible non-isotropy effects into account.

Here, we demonstrate that the flow reversal found in recent VSI-turbulent disks (Stoll & Kley, 2016) can in fact be traced back to the intrinsic anisotropy of the VSI-turbulence. Using multidimensional hydrodynamical simulations, we calculate the effective radial and vertical transport coefficients and use them in viscous disk simulations.

2 Model setup

We use the pluto code (Mignone et al., 2007) for our simulations, where we model a section of a locally isothermal accretion disk. For the direct turbulent simulations we use a full 3D setup and spherical coordinates ($R$,$\theta$,$\phi$), while for comparison laminar viscous simulations ($\alpha$-disk), we use a 2D axisymmetric setup and cylindrical coordinates ($r,z$). The model parameters are given in Table 1.

Even though some of the simulations are performed in spherical coordinates, the results are analyzed in cylindrical coordinates. The initial disk is axisymmetric and extends from $0.4\ r_{\mathrm{0}}$ to $2.5\ r_{\mathrm{0}}$, where $r_{0}=5.2\,\mathrm{au}$. The initial density profile is given by vertical hydrostatic equilibrium

[TABLE]

where $\rho_{0}$ is the gas mid-plane density at $r=r_{0}$, and $p=-1.5$. In our locally isothermal approximation the temperature of the disk is a function of the cylindrical radius only,

[TABLE]

where we choose $q=-1,$ which causes the disk aspect ratio to be constant and $T_{0}$ such that $H/r=h=0.05$. The pressure is given by $P=c_{s}^{2}\rho$, where $c_{s}$ is the isothermal sound speed with $c_{s}\propto r^{-1/2}$. The gas moves initially with the angular velocity given by the Keplerian value, corrected by the pressure support (Nelson et al., 2013)

[TABLE]

with $\Omega^{2}_{K}(r)=GM_{*}/r^{3}$, and the meridional flow ($u_{r},u_{z}$) is set to zero. At the inner and outer boundary we use reflecting boundaries. To increase numerical stability for the turbulent runs, we damp the variables $\rho,u_{r},u_{z}$ to the initial values on a timescale of half a local orbit following the recipe by de Val-Borro et al. (2006), with a damping applied at the inner boundary from 0.4-0.5 $r_{0}$ and at the outer boundary from 2.3-2.5 $r_{0}$. At the vertical boundaries we use reflective boundaries when the flow is directed inward and zero gradient otherwise. The VSI model is inviscid, and for the $\alpha$-disk model the kinematic viscosity is given by $\nu=2/3\alpha c_{s}H$, with a constant $\alpha=5\cdot 10^{-4}$, which matches the outcome of VSI model.

3 Disk structure

We compare the structure of the VSI unstable disk to viscous disks described by an $\alpha$-parameter, and study the main differences. The VSI disk is evolved to a quasi-equilibrium state before the analysis is performed. We focus in particular on the flow field and stress tensor, and refer to Nelson et al. (2013) and Stoll & Kley (2014, 2016) for an analysis of the turbulent flow structure.

3.1 Radial velocity

While the density distributions of the viscous and turbulent disks in equilibrium are very similar due to the necessary pressure equilibrium, there is an important difference in the meridional flow within the disk in particular for the radial velocity. We show the azimuthally averaged radial velocity $u_{r}$ at $r=r_{0}$ as a function of the vertical distance for the turbulent and the standard viscous disk in Fig. 1. For the turbulent disk we averaged the velocity in time over 50 orbits and in space around $r_{0}$ in the region $(0.8-1.25)r_{0}$. For better visibility we rescaled the viscous case by a factor of $200$. Obviously, the two radial velocity profiles have an opposite behavior. The standard viscous disk using a single (isotropic) value of $\alpha$ shows the typical outflow in the midplane that has been predicted analytically by Urpin (1984) and was later shown in fully time-dependent numerical simulations (Kley & Lin, 1992). This behavior of $u_{r}(z)$ can be derived from the equilibrium angular momentum equation that contains the vertical disk structure (Urpin, 1984). In spite of the outwardly directed flow in the midplane, the total vertically integrated mass flow is nevertheless directed inward in case of accreting disks (Kley & Lin, 1992).

In contrast, the mean flow field for the turbulent flow (labeled 3D VSI and shown also with a dotted line for a model with double resolution in all directions) is fully reversed; it is not only negative in the midplane and positive in the corona, but also much stronger, as indicated by the different scaling of the curves in Fig. 1. This special feature of the meridional flow field in the VSI case has been found for isothermal as well as radiative disks in Stoll & Kley (2016), but was not analyzed with respect to anisotropic turbulence. From the direct comparison to the viscous case, it is clear that with a standard shear viscosity prescription using a constant $\alpha$-value or a constant kinematic viscosity, no agreement can be obtained because this will always lead to an inverted parabolic type of profile (Jacquet, 2013). This raises the general question whether the VSI turbulence in disks can be described by a standard Navier-Stokes approach to model the angular momentum diffusion. In Sect. 4 we show that we obtain a good match of a fully turbulent and viscous flow for a non-anisotropic turbulent viscosity where the radial and vertical parts enter with a different strength (see curve 2D anisotropic stress), which might be expected for the clearly non-isotropic character of the flow structure in VSI turbulence (Stoll & Kley, 2014).

3.2 Turbulent stresses

To analyze the effect of the turbulence with respect to angular momentum transport in the disk, it is necessary to calculate the turbulent stresses of the VSI disk. For the overall mass flow in accretion disks, which is driven by outward angular momentum transport, the $r\phi$-component of the Reynolds stress tensor, $R$, is the most important component because it is generated by the strong shear in the azimuthal velocity. In the case of VSI turbulent disks, it is clear that the vertical dependence of the turbulent stresses ($z\phi$-component) may be of importance as well. Hence, from our simulations we calculate the following turbulent Reynolds stresses

[TABLE]

where $u_{i}$ with $i$ in $(r,z)$ denote the radial and vertical flow velocity ($u_{r},u_{z}$) and $\delta u_{\phi}$ is the deviation of the azimuthal velocity, $u_{\phi}=r\Omega$, from unperturbed equilibrium, as given by Eq. (3). For the later analysis it is beneficial that the stresses are calculated in cylindrical coordinates ($r,z$), see Sect. 4. Since we are interested in particular in the vertical dependence of the stresses, $R_{i,\phi}$ is calculated by averaging over the full azimuth ($2\pi$), over a small radial range around a reference radius, and over time. Here, we time-average from orbit 80 to orbit 200 using a series of over 60 snap shots and space-average in radius from 0.75 - 1.35 $r_{0}$.

The results of this averaging procedure are shown in Fig. 2. The solid curves refer to the specific stresses $R_{i,\phi}(z)/\rho(r_{0},z)$ ($i=r,z$) in units of $c_{s0}^{2}$, where $\rho(r_{0},z)$ is the equilibrium density distribution, given by Eq. (1), and $c_{s0}$ is the sound speed, both at the reference radius $r_{0}$.

We compare it to the viscous shear stress prescription (see Eq. 6 with $\nu=\nu_{r}=\nu_{z}$) using $\alpha=5\cdot 10^{-4}$ (dotted lines), which is close to the average of the $r\phi$-component of the specific VSI stress, $4\cdot 10^{-4}$. We rescaled the $z\phi$-component of the viscous shear stress by a factor of 650 to match the value from the VSI stress. From this large rescaling factor it is immediately clear that $R_{z\phi}$ is far larger than expected from an isotropic viscous shear prescription. We note that the deviation of the vertical $R_{r\phi}/\rho$-profile from the constant $S_{r\phi}/\rho$-profile is not important for our argument, which is why for simplicity we chose a constant $\alpha$. For the angular momentum transport only the vertical average of this stress component plays a role. In our case, the very large $z\phi$-component of the stress dominates the meridional flow and the influence of $R_{r\phi}$ is small.

4 Anisotropic viscosity of the VSI turbulence

From our numerical studies of the VSI turbulence, in particular the vertical component of the Reynolds stress tensor, we can infer that the turbulence is non-isotropic. Even though we have integrated the hydrodynamical equations using spherical coordinates, we calculated and display $R_{i\phi}$ in cylindrical coordinates because this is simpler to analyze, as we explain in the following. In our approach we investigate whether the turbulent Reynolds stresses, $R$, can be modeled by a viscous ansatz where $R$ is replaced by the standard viscous stress tensor, $S$, with an effective turbulent eddy viscosity as introduced by Boussinesq. In cylindrical coordinates the change in angular momentum is given by the following evolution equation:

[TABLE]

For an axisymmetric flow, that is, when $\partial/\partial\phi$ vanishes, the vector $\@vec{S}_{\phi}$ (the $\phi$-row of the viscous stress tensor) is given by (Tassoul, 1978)

[TABLE]

Axisymmetry is expected for accretion disk flows, and we found this in our simulations. In Eq.(6) we allowed for the option of an anisotropic viscosity by splitting the kinematic viscosity into two components $\nu_{r}$ and $\nu_{z},$ where $\nu_{r}$ refers to the radial part, which typically is the main contribution in accretion disks in driving the angular momentum transport and mass accretion. The $\nu_{z}$ part connects to the vertical variation of the angular velocity $\Omega$ and has commonly been assumed to be on the same order as $\nu_{r}$. Indeed, using $\nu_{z}=\nu_{r}$ and performing viscous axisymmetric, two-dimensional simulations, we find the typical meridional flow field in disks with outflow in the midplane and inflow in upper layers of the disk, as shown in Fig. 1 by the curve labeled ’isotropic stress’, as found in classic studies (Urpin, 1984; Kley & Lin, 1992). The outflow in the midplane can easily be derived by an analysis of the radial flow that can be obtained from the angular momentum equation (see also Fromang et al., 2011; Jacquet, 2013). From Eq. (5) we find for the equilibrium state that

[TABLE]

Expanding Eq. (3) around the midplane (small $z$), we find

[TABLE]

where $q$ and $p$ are the exponents in the radial power laws of the density and temperature, respectively, and $h=H/r$ denotes the relative scale height of the disk. Using this relation for $\Omega(r,z),$ we find

[TABLE]

and then, neglecting terms of order $z^{2}/r^{2}$

[TABLE]

Combining this with the $S_{z\phi}$ relation, we find

[TABLE]

This is plotted in the lower panel in Fig. 2 using $\nu_{z}=650\nu_{r}$, which demonstrates that the linear dependence of the specific $S_{z\phi}$-stress is a direct consequence of Eq. (11). From Eq. (7) we find for the relation determining the sign of $u_{r}$ in the midplane the following relation:

[TABLE]

where for the last step we assumed that the viscosity has an $\alpha$-type behavior with $\nu\sim\alpha c_{s}^{2}/\Omega_{K}$ and $\rho\propto r^{p}$ and $c_{s}^{2}\propto r^{q}$. We note that the direction of flow in the disk midplane can be influenced by the slopes ($p,q$) in the disk stratification (Fromang et al., 2011; Philippov & Rafikov, 2017). From the above relation we find directly that for our disk with slopes $q=-1$ and $p=-3/2$ that

[TABLE]

For isotropic turbulence with $\alpha_{r}=\alpha_{z}$ we therefore have outflow in the disk midplane. For the turning point we find that $\alpha_{z}$ must be larger than $1.5\alpha_{r}$. Upon increasing $\alpha_{z}$ over $\alpha_{r},$ the midplane radial velocity becomes more and more negative, and the entire vertical flow profile eventually reverses. In our numerical simulations of viscous disks using different values for $\nu_{z}$ we could indeed find the observed flow reversal for moderate values of $\nu_{z}/\nu_{r}$. By increasing $\nu_{z}$ further, we found that the VSI turbulence can be modeled by an anisotropic eddy viscosity with $\nu_{z}$ over two magnitudes larger than $\nu_{r}$ (650 for our chosen parameter, as shown in Fig. 2).

In performing the comparison simulations, we initially tried to use the same spherical coordinate system that we used for the turbulent VSI simulations by just increasing $\nu_{\theta}$ over $\nu_{R}$ in the corresponding components of the stress tensor in spherical coordinates. However, this did not lead to the expected results, in particular, the inversion of the parabola of $u_{r}(z)$. The results displayed in Fig. 1 were therefore obtained using a cylindrical coordinate system. Because the usage of a spherical coordinate system is beneficial in disk simulations, the equations need to be transformed from $r,z$ to $R,\theta$.

This means we need to transform the viscous part, $\nabla\cdot(r\@vec{S}_{\phi}$), to spherical coordinates. The transformation equations are

[TABLE]

and the components of the stress tensor transform according to

[TABLE]

The derivatives of $\Omega$ transform then as

[TABLE]

Using Eq. (6) in Eqs. (15) and (16) and substituting this into Eqs. (17) and (18), we obtain for the stress tensor components in spherical coordinates

[TABLE]

As can be inferred from these equations, the relatively simple anisotropic relation in cylindrical coordinates leads to complex coupled equations in spherical coordinates with cross terms in the derivatives for $\Omega$. For the isotropic case we can set $\nu_{r}=\nu_{z}=\nu$ and obtain the standard relations for the stress tensor components in $R\theta$-coordinates (Tassoul, 1978). Now it becomes clear why merely increasing the $\theta$ part of the viscosity in $R\theta$-coordinates did not yield the correct answer when compared to the VSI case. Using the full viscosity terms as described in Eqs. (19) and (20) in the numerical simulations yields the correct behavior. However, this comes with a serious drawback because in this case the numerical integration required much smaller time steps for numerical stability than the cylindrical simulation with the same viscosity coefficients. Performing multidimensional simulations of viscous disks mimicking the non-isotropic behavior of the VSI turbulence in $R\theta$-coordinates therefore requires a considerable numerical effort. One solution to this problem may be the usage of an implicit solver for the viscous terms.

5 Conclusions

From direct three-dimensional simulations of locally isothermal accretion disks we observed that the eddies introduced by the VSI generate stresses that are strongly anisotropic. Specifically, we found that the vertical $z\phi$-component of the Reynolds stress is enhanced by a factor of 650 over the standard $r\phi$-part. By performing viscous disk simulations using a non-isotropic viscosity with $\alpha_{z}$ highly enlarged over $\alpha_{r},$ we could obtain the same flow reversal as seen in the VSI disk, which verifies the non-isotropy of the viscosity. Hence, the reversal of the radial flow profile compared to the usual $\alpha$-model is a clear consequence of the anisotropy. This will have an effect on the dust migration processes (Stoll & Kley, 2016) that need to be contrasted to the outward drift in the midplane viscous models (Takeuchi & Lin, 2002). From this we can conclude that we need to be careful with turbulence models imposed on accretions disks when we adopt viscous models to describe them.

The meridional circulation of MRI-turbulent disks has been analyzed by Fromang et al. (2011), who found a similar mean flow dynamics in the disk but did not attribute it to a non-isotropic turbulence (with large $S_{z\phi}$) but rather to a radial variation in the magnitude of the viscosity. However, their turbulent $R_{z\phi}$ profile (in Fig. 5) indicates that it may be enhanced over the standard viscous value as well, which could also be a reason for the flow reversal.

Even though the analysis is performed for simplicity for a locally isothermal disk, our results are quite general as simulations for fully radiative disks show the same behavior although they have different radial and vertical temperature profiles. For this purpose, we reanalyzed our radiative simulations in Stoll & Kley (2016) and found a similar anisotropy. Concerning numerical resolution, simulations with double resolution show the same results as shown in Fig. 1. Nevertheless, further exploration needs to be done in order to check how the anisotropy factor varies for different disk parameters and to explore the possibility of anisotropic stresses in MRI models.

Acknowledgements.

Moritz Stoll received financial support from the Landesgraduiertenförderung of the state of Baden-Württemberg and through the Ger man Research Foundation (DFG) grant KL 650/16. G. Picogna acknowledges the support through DFG-grant KL 650/21 within the collaborative research program “The first 10 Million Years of the Solar System”. Some simulations were performed on the bwGRiD cluster in Tübingen, funded by the state of Baden-Württemberg and the DFG. We thank Roman Rafikov for providing us with a preprint and very useful discussions.