# On the angular momentum transport due to vertical convection in accretion discs

Geoffroy Lesur, G. Ogilvie

► To cite this version:

Geoffroy Lesur, G. Ogilvie. On the angular momentum transport due to vertical convection in accretion discs. Monthly Notices of the Royal Astronomical Society, 2010, 404, pp.L64. hal-00631988

HAL Id: hal-00631988

https://hal.science/hal-00631988v1

Submitted on 22 Nov 2021

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# On the angular momentum transport due to vertical convection in accretion discs

Geoffroy Lesur\* and Gordon I. Ogilvie

Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA

Accepted 2010 February 22. Received 2010 February 22; in original form 2010 January 8

## ABSTRACT

The mechanism of angular momentum transport in accretion discs has long been debated. Although the magnetorotational instability appears to be a promising process, poorly ionized regions of accretion discs may not undergo this instability. In this Letter, we revisit the possibility of transporting angular momentum by turbulent thermal convection. Using high-resolution spectral methods, we show that strongly turbulent convection can drive outward angular momentum transport at a rate that is, under certain conditions, compatible with observations of discs. We find, however, that the angular momentum transport is always much weaker than the vertical heat transport. These results indicate that convection might be another way to explain global disc evolution, provided that a sufficiently unstable vertical temperature profile can be maintained.

**Key words:** accretion, accretion discs – hydrodynamics – turbulence.

## 1 INTRODUCTION

As shown by the early models of Shakura & Sunyaev (1973) and Lynden-Bell & Pringle (1974), the structure and evolution of accretion discs are controlled by the process of angular momentum transport. It is often assumed that turbulence is responsible for this transport, although its nature and physical origin are barely understood. The most straightforward way to obtain a turbulent disc is to find a linear instability of the basic flow, which could drive turbulence in the non-linear regime. However, the presence of a turbulent flow does not necessarily imply an outward transport of angular momentum through the disc, if the turbulence is sustained by sources of energy other than the Keplerian shear.

The magnetorotational (MRI) instability, pointed out in the context of accretion discs by Balbus & Hawley (1991), is believed to play an important role in accretion discs, as it produces turbulence that efficiently transports angular momentum outwards (e.g. Balbus 2003). However, the behaviour of the MRI in poorly ionized gas is uncertain. In some regions of discs, the instability might cease to exist because the resistivity is too large, leading to the concept of ‘dead zones’ (Gammie 1996). On the other hand, purely hydrodynamical processes might also be at work in discs, such as the subcritical baroclinic instability (Klahr & Bodenheimer 2003; Petersen, Stewart & Julien 2007; Lesur & Papaloizou 2010), driven by radial entropy gradients.

The presence of a vertical convective instability in discs was initially suggested by Lin & Papaloizou (1980) as a way of transporting angular momentum. They pointed out that the temperature

dependence of the opacity in protostellar (or protoplanetary) discs naturally leads to convective instability in certain regimes. However, subsequent studies, using either linear approaches (Ryu & Goodman 1992) or non-linear simulations (Cabot 1996; Stone & Balbus 1996), predicted a weak inward transport associated with convection, leading to the conclusion that this process was irrelevant for the problem of accretion disc dynamics. Note that these simulations were performed in a ‘weakly non-linear’ regime with relatively low Reynolds and Rayleigh numbers (Cabot 1996).

In this Letter, we challenge these conclusions using high-resolution spectral methods and a simplified setup for the convection problem. In Section 2, we introduce our model, the relevant dimensionless parameters and the numerical setup. We present simulations of convection with free-slip vertical boundary conditions in Section 3 and of homogeneous convection in Section 4. Finally, we briefly discuss these results and their implications in Section 5.

## 2 MODEL AND EQUATIONS

### 2.1 Model

In the following, we consider a model flow in which rotation, shear and thermally induced convection are included. Instead of considering the full disc problem, with its particular density and temperature structure, we use the shearing box approximation (Hawley, Gammie & Balbus 1995; Balbus 2003; Regev & Umurhan 2008). As a simplification, we assume the flow is incompressible and we introduce vertical stratification within the Boussinesq approximation (Spiegel & Veronis 1960). Although this approximation is not fully justified in a real disc, it allows us to explore the regime of strongly turbulent

<footnote>\*E-mail: g.lesur@damtp.cam.ac.uk</footnote>
sheared convection in a more controlled way. In this Letter, we explore the way turbulent convection transports angular momentum, but we do not include the generation of an unstable temperature profile due to viscous heating. Instead, we impose a 'background' stratification profile as prescribed by the Boussinesq approximation. Therefore, one has to keep in mind that convection is not generated self-consistently in our model.

The shearing-box equations are found by considering a Cartesian box centred at $r = R_0$, rotating with the disc at angular velocity $\Omega = \Omega(R_0)$. Defining $r - R_0 \to x$ and $R_0\phi \to y$ as in Hawley et al. (1995), one obtains the following set of equations:

$$ \partial_t u + u \cdot \nabla u = -\nabla \Pi - 2\Omega \times u + 2\Omega Sx \mathbf{e}_x - \Lambda N^2 \theta \mathbf{e}_z + \nu \Delta u, \quad (1) $$

$$ \partial_t \theta + u \cdot \nabla \theta = u_z / \Lambda + \chi \Delta \theta, \quad (2) $$

$$ \nabla \cdot u = 0, \quad (3) $$

where $u$ is the total velocity, $\theta \equiv \delta\rho/\rho$ is the perturbation of the density logarithm (or entropy, as pressure perturbations are much smaller in the Boussinesq approximation), $\nu$ is the kinematic viscosity and $\chi$ is the thermal diffusivity. In these equations, we have defined the mean shear $S = -r\partial_t\Omega$, which is set to $S = (3/2)\Omega$ for a Keplerian disc. The generalized pressure $\Pi$ is calculated by solving a Poisson equation derived from the incompressibility condition. For dimensional consistency with the traditional Boussinesq approach, we have introduced a stratification length $\Lambda \equiv -g/N^2$, where $g$ is the vertical component of the gravity. Note, however, that $\Lambda$ disappears from the dynamical properties of this system of equations as one can renormalize the variables by defining $\theta' \equiv \Lambda\theta$. The stratification itself is controlled by the Brunt-Väisälä frequency $N$, defined for a perfect gas by

$$ N^2 = -\frac{1}{\gamma \rho} \frac{\partial P}{\partial z} \frac{\partial}{\partial z} \ln \left( \frac{P}{\rho^\gamma} \right), \quad (4) $$

where $P$ and $\rho$ are the pressure and density of the background equilibrium profile and $\gamma$ is the adiabatic index. Physically, $N$ corresponds to the oscillation frequency of a fluid particle displaced vertically from its equilibrium position. When the flow is convectively unstable, one has $N^2 < 0$. In our model, we assume $N$ is constant, which formally corresponds to a local model in $z$.

One can easily check that the velocity field $U = -Sxe_y$ is a steady solution of equations (1)–(2). In the following, we consider the evolution of the perturbations $v$ (not necessarily small) of this profile defined by $v = u - U$.

In the following, we use shearing-sheet boundary conditions (Hawley et al. 1995) in the $x$ direction and periodic boundary conditions in $y$. In the $z$ direction, we use either free-slip boundary conditions (Section 3), imposing $v_z = \theta = 0$ and $\partial_z v_x = \partial_z v_y = 0$, or periodic boundary conditions (Section 4). The free-slip vertical boundary conditions correspond to a rigid wall with no tangential stress and a fixed temperature. Similar boundary conditions for the velocity field were used by Cabot (1996).

## 2.2 Dimensionless numbers

To simplify the analysis of convection in sheared flow, we will use three dimensionless numbers defined as follows.

(i) The Richardson number $Ri = -N^2/S^2$ compares buoyancy forces to the mean shear. In a real disc, this number is determined by the balance between viscous heating and vertical heat transport.

Assuming a typical velocity $v$ due to convective motions, we define the turbulent viscosity using mixing length theory as $v_t \sim \delta vH$, where $\delta$ is an efficiency factor and $H$ is the disc thickness. The turbulent heat transport is on the other hand $F \sim \rho v^3$. Assuming thermal balance, we get $\rho v^3/H \sim \rho v_t \Omega^2$ leading to $v^2 \sim \delta c_s^2$ where $c_s$ is the local sound speed. Since turbulence is driven by the unstable entropy gradient, we have in first approximation $v \sim |N|H$, so that $Ri \sim \delta^{1/2}$. As we will see in the following, we generally get $\delta \ll 1$, meaning that the regime $Ri < 1$ is relevant to accretion discs.

(ii) The Prandtl number $Pr = v/\chi$ compares viscous dissipation to thermal diffusion. Using standard values for the viscosity and thermal conductivity of $H_2$, relevant to a protoplanetary disc, one has

$$ Pr = 1.5 \times 10^{-9} \left( \frac{T}{300 K} \right)^{-5/2} \left( \frac{\kappa}{1 \text{ cm}^2 \text{ g}^{-1}} \right) \left( \frac{\rho}{10^{-9} \text{ g cm}^{-3}} \right), $$

where $\kappa$ is the opacity.

(iii) The Rayleigh number $Ra = -N^2 L^4 / \chi \nu$ compares buoyancy forces to dissipation processes. In a disc, one gets

$$ \frac{Ra}{Ri} = 9 \times 10^{20} \left(\frac{L}{r}\right)^4 \left(\frac{r}{\text{AU}}\right) \left(\frac{\kappa}{1 \text{ cm}^2 \text{ g}^{-1}}\right) \left(\frac{M}{M_{\odot}}\right) \times \left(\frac{T}{300 \text{ K}}\right)^{-7/2} \left(\frac{\rho}{10^{-9} \text{ g cm}^{-3}}\right)^3. $$

In the following, we will assume $Pr = 1$ for simplicity. A more complete study including the effect of varying $Pr$ will be the subject of a subsequent paper. The length $L$ introduced in these definitions is assumed to be a typical dimension of the box and should be of the order of the disc thickness. In the following, we define $L$ as the smallest box length with $L = L_z$ in Section 3 and $L = L_x$ in Section 4. When not explicitly mentioned, the unit of time is the shear time-scale $S^{-1}$ which is equal to $(3\pi)^{-1}$ orbital periods. In the following, all of our simulations are computed in the shear-dominated regime ($Ri < 1$) which is expected in discs.

## 2.3 Turbulent transport and energy budget

By analogy with Stone & Balbus (1996), we define the dimensionless transport of angular momentum through the $\alpha$-like coefficient $\alpha \equiv \langle v_x v_y \rangle / S^2 L^2$, $\langle \cdot \rangle$ being a time and space average. Since $S = 1$ and $L = 1$ in all the simulations presented below, we will simply write $\alpha = \langle v_x v_y \rangle$ in our units.

Multiplying the $y$ component of (1) by $v_y$, integrating over the volume of the box and using the boundary conditions presented above, we obtain the energy budget of the fluctuations in the $y$ direction,

$$ \partial_t \overline{v_y^2} = \overline{\Pi \partial_y v_y} + (S - 2\Omega) \overline{v_x v_y} - \nu |\overline{\nabla v_y}|^2, \quad (5) $$

where $\overline{X}$ denotes the volume average of $X$. Considering a quasi-steady turbulence and assuming the pressure-strain correlation $\overline{\Pi \partial_y v_y}$ to be negligible, Stone & Balbus (1996) concluded from this equation that $(S - 2\Omega)\overline{v_x v_y} > 0$, leading to $\alpha < 0$ in a disc. This hypothesis was motivated by numerical simulations (Stone & Balbus 1996; Cabot 1996) in which the turbulence was weakly non-axisymmetric, leading to the conclusion that $y$ derivatives should be negligible. It is, however, well known that the pressure-strain correlation tensor is responsible for the redistribution of energy in fully developed anisotropic turbulence (Speziale 1991; Pope 2000). Therefore, when the flow becomes strongly non-linear (i.e. at large $Ra$), there is no reason to assume that this term is negligible.
<table><caption>**Table 1.** List of simulations with free-slip boundary conditions. The turbulent kinetic energy (sixth column), angular momentum transport (seventh column) and heat transport (eighth column) are computed from simulations.</caption><thead><tr><th>Run</th><th>Resolution<br>n<sub>x</sub> × n<sub>y</sub> × n<sub>z</sub></th><th>Box size<br>L<sub>x</sub> × L<sub>y</sub> × L<sub>z</sub></th><th>Ra</th><th>Ri</th><th>⟨½v<sup>2</sup>⟩/S<sup>2</sup>L<sup>2</sup></th><th>⟨v<sub>x</sub>v<sub>y</sub>⟩/S<sup>2</sup>L<sup>2</sup> ≡ α</th><th>Λ⟨θv<sub>z</sub>⟩/SL<sup>2</sup></th></tr></thead><tbody><tr><td>W1</td><td>256 × 256 × 128</td><td>2 × 4 × 1</td><td>1 × 10<sup>5</sup></td><td>0.4</td><td>5.7 × 10<sup>-3</sup></td><td>-3.9 × 10<sup>-4</sup></td><td>7.7 × 10<sup>-3</sup></td></tr><tr><td>W2</td><td>256 × 256 × 128</td><td>2 × 4 × 1</td><td>3 × 10<sup>5</sup></td><td>0.4</td><td>7.0 × 10<sup>-3</sup></td><td>-2.3 × 10<sup>-4</sup></td><td>7.4 × 10<sup>-3</sup></td></tr><tr><td>W3</td><td>256 × 256 × 128</td><td>2 × 4 × 1</td><td>6 × 10<sup>5</sup></td><td>0.4</td><td>6.9 × 10<sup>-3</sup></td><td>-1.4 × 10<sup>-4</sup></td><td>7.0 × 10<sup>-3</sup></td></tr><tr><td>W4</td><td>256 × 256 × 128</td><td>2 × 4 × 1</td><td>1.2 × 10<sup>6</sup></td><td>0.4</td><td>7.2 × 10<sup>-3</sup></td><td>-4.6 × 10<sup>-5</sup></td><td>6.2 × 10<sup>-3</sup></td></tr><tr><td>W5</td><td>256 × 256 × 128</td><td>2 × 4 × 1</td><td>2.5 × 10<sup>6</sup></td><td>0.4</td><td>7.2 × 10<sup>-3</sup></td><td>5.8 × 10<sup>-5</sup></td><td>5.7 × 10<sup>-3</sup></td></tr><tr><td>W6</td><td>256 × 256 × 128</td><td>2 × 4 × 1</td><td>5 × 10<sup>6</sup></td><td>0.4</td><td>7.1 × 10<sup>-3</sup></td><td>1.4 × 10<sup>-4</sup></td><td>5.1 × 10<sup>-3</sup></td></tr><tr><td>W7</td><td>256 × 256 × 128</td><td>2 × 4 × 1</td><td>1 × 10<sup>7</sup></td><td>0.4</td><td>6.7 × 10<sup>-3</sup></td><td>2.1 × 10<sup>-4</sup></td><td>4.6 × 10<sup>-3</sup></td></tr><tr><td>W8</td><td>256 × 256 × 128</td><td>2 × 4 × 1</td><td>1.4 × 10<sup>7</sup></td><td>0.4</td><td>6.4 × 10<sup>-3</sup></td><td>2.4 × 10<sup>-4</sup></td><td>4.3 × 10<sup>-3</sup></td></tr><tr><td>W9</td><td>256 × 256 × 128</td><td>2 × 4 × 1</td><td>3 × 10<sup>5</sup></td><td>0.2</td><td>2.1 × 10<sup>-3</sup></td><td>-1.6 × 10<sup>-4</sup></td><td>3.9 × 10<sup>-3</sup></td></tr><tr><td>W10</td><td>256 × 256 × 128</td><td>2 × 4 × 1</td><td>1.2 × 10<sup>6</sup></td><td>0.2</td><td>2.2 × 10<sup>-3</sup></td><td>-1.0 × 10<sup>-4</sup></td><td>3.7 × 10<sup>-3</sup></td></tr><tr><td>W11</td><td>256 × 256 × 128</td><td>2 × 4 × 1</td><td>5 × 10<sup>6</sup></td><td>0.2</td><td>2.2 × 10<sup>-3</sup></td><td>-5.3 × 10<sup>-5</sup></td><td>3.1 × 10<sup>-3</sup></td></tr><tr><td>W12</td><td>256 × 256 × 128</td><td>2 × 4 × 1</td><td>2 × 10<sup>7</sup></td><td>0.2</td><td>2.0 × 10<sup>-3</sup></td><td>-8.7 × 10<sup>-6</sup></td><td>2.5 × 10<sup>-3</sup></td></tr><tr><td>W13</td><td>512 × 512 × 256</td><td>2 × 4 × 1</td><td>8 × 10<sup>7</sup></td><td>0.2</td><td>1.8 × 10<sup>-3</sup></td><td>3.0 × 10<sup>-5</sup></td><td>1.9 × 10<sup>-3</sup></td></tr><tr><td>W14</td><td>512 × 512 × 256</td><td>2 × 4 × 1</td><td>3.2 × 10<sup>8</sup></td><td>0.2</td><td>1.7 × 10<sup>-3</sup></td><td>5.8 × 10<sup>-5</sup></td><td>1.6 × 10<sup>-3</sup></td></tr></tbody></table>

## 2.4 Numerical model

We use the Snoopy code to solve equations (1)–(3). This code uses a 3D Fourier representation of the flow with a third-order time-stepping scheme to have an excellent accuracy at all scales. It has been used in several studies of hydrodynamic and MHD turbulence (see e.g. Lesur & Longaretti 2005, 2007). When working with free-slip boundary conditions, a sine–cosine decomposition is used instead of the Fourier basis in the z direction (see e.g. Cattaneo, Emonet & Weiss 2003).

## 3 CONVECTION WITH WALLS

In the first set of simulations, we introduce free-slip boundary conditions in the z direction. These boundary conditions prevent the development of large-scale vertical flows which are unrealistic for a disc (see Section 4). In the particular setup considered here, convection starts at $Ra \approx 4100$ for $Ri = 0.4$ and $Ra \approx 6900$ for $Ri = 0.2$.

In order to study the dependence of the turbulent transport of angular momentum on our dimensionless parameters, we have carried out several simulations which are presented in Table 1. Average values are computed from simulations run for $500 S^{-1}$, excluding the first $100 S^{-1}$ from the averaging procedure to eliminate the influence of the initial conditions. As shown in Fig. 1, the simulations are sufficiently homogeneous on these time-scales for this average to be meaningful.

### 3.1 Rayleigh number and transport

The most striking observation in these simulations is the correlation between $Ra$ and the turbulent transport $\alpha$. In particular, for $Ra < 10^6$ we recover the result of Stone & Balbus (1996) of a small and negative transport (runs W1–W4 and W9–W12). However, for sufficiently large $Ra$, the sign of the transport is reversed. To demonstrate this effect, we show in Fig. 2 the dependence of the turbulent transport on $Ra$, for $Ri = 0.4$ and $Ri = 0.2$.

It is possible to deduce a general scaling for the transport from these results:

$$ \alpha \approx -1.8 \times 10^{-3} + 2.9 \times 10^{-4} \log_{10}(Ra) \quad (6) $$

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Turbulent transport as a function of time for run W1 (lower curve) and W8 (upper curve). The modification of the average turbulent transport is due to the Rayleigh number (see the text).</figcaption></figure>

for $Ri = 0.4$ and

$$ \alpha \approx -5.4 \times 10^{-4} + 7.2 \times 10^{-5} \log_{10}(Ra) \quad (7) $$

for $Ri = 0.2$. Surprisingly, these scalings predict a transport of the order of $10^{-3}$ for values of $Ra$ which could be expected in a disc. However, this conclusion is challenged by the fact that at large $Ra$, the dependence of $\alpha$ on $Ra$ should become weaker since the length-scales associated with transport are then much larger than the diffusive scales and the diffusion coefficients may be expected to become less important. A weak effect of this kind is already observed in Fig. 2 for $Ra > 10^7$.

Interestingly, changing the Rayleigh number does not significantly modify the turbulent kinetic energy of the system. This clearly demonstrates that the correlation observed above is not due to a modification of the turbulent activity in the flow but results instead from a change in the anisotropy of the turbulence. Note also that the turbulent heat transport $\langle\theta v_z\rangle$ decreases as $Ra$ increases. However, we always find $\langle\theta v_z\rangle \gg \langle v_x v_y \rangle$, showing that convection is inefficient at transporting angular momentum. This point will be examined further in Section 5.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Correlation between the turbulent transport and Ra for Ri = 0.4 (plain line) and Ri = 0.2 (dashed line).</figcaption></figure>

## 3.2 Energy budget

To check if the pressure–strain correlation is important in our simulations, we plot in Fig. 3 the time history of each term of equation (5) for run W7. As shown in this plot, the energy balance at high Ra is dominated by an equilibrium between the pressure–strain correlation and the dissipation term. Moreover, the turbulent transport term is negligible in this situation. This clearly indicates that at high enough Ra the hypothesis of Stone & Balbus (1996) is no longer valid and the turbulent transport of angular momentum can be directed outwards.

The presence of a large pressure–strain correlation term in equation (5) shows that non-axisymmetric structures are playing an important role at large Ra. To isolate the role played by these non-axisymmetric structures, we plot in Fig. 4 the 1D spectrum of the pressure–strain correlation. As expected, we find no contribution from axisymmetric modes ($k_y = 0$). However, as Ra becomes larger, the contribution from smaller-scale non-axisymmetric modes becomes more important. This effect is evident in run W7 where contributions from scales as small as $L_y/20$ are still important.

This result clearly indicates that as the flow becomes turbulent, it develops small-scale motions which tend to redistribute the energy

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Energy budget in the y direction (equation 5) from run W7.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** 1D spectra of the pressure–strain correlation from run W1 (plain line), W5 (dashed line) and W7 (dash-dotted line).</figcaption></figure>

fluctuations in the three directions. This is consistent with the picture of anisotropic turbulence becoming isotropic at small scales thanks to the pressure–strain correlation. It also indicates that the results of Cabot (1996) and Stone & Balbus (1996) were due to a too small Ra, an issue already mentioned by Cabot (1996).

## 4 HOMOGENEOUS CONVECTION

To test the possible limitations of the free-slip boundary conditions used above, we have also considered the problem of homogeneous convection, using periodic boundary conditions in the vertical direction applied to $v$ and $\theta$ (Borue & Orszag 1997). In this context, the flow becomes convectively unstable for $Ra > 1558$, developing first ‘elevator’ modes with a purely vertical motion $v_z \propto \cos(k_x x + \phi)$. Since these modes are also exact non-linear solutions of equations (1)–(3), they can grow exponentially for long times, leading to ‘spikes’ in the turbulent energy and transport. These exponential growth phases are unrealistic, as they should be broken by secondary instabilities or the effects of density stratification, which are neglected in the Boussinesq approximation. A way to limit the formation and amplification of these modes is to consider ‘tall’ boxes with $L_z > L_x$, so that secondary instabilities are more efficient at breaking up the elevator modes. Note that the phenomenology of these elevator modes is very similar to that of channel mode solutions found in the context of the MRI instability (see Goodman & Xu 1994).

We have used this tall box setup with several values of Ra and $L_y$ in the regime of weak convection, $Ri = 0.1$ (Table 2). Looking at the time series (Fig. 5), we note the presence of large peaks of transport at small Ra, with both positive and negative values. These peaks are associated with the growth of elevator modes to high amplitude and their break-up through secondary instabilities. We note, however, that this behaviour disappears at larger Ra, although large fluctuations remain, leading to a positive-definite turbulent transport with an average value $\alpha \sim 10^{-3}$.

Modifying the size of the box in the y direction does not change dramatically the results, and on average, we find a transport that is stronger in the homogeneous case than with free-slip boundary conditions. Note also that the average turbulent kinetic energy is
<table><caption>**Table 2.** List of simulations with periodic boundary conditions in the z direction.</caption><thead><tr><th>Run</th><th>Resolution<br>n<sub>x</sub> × n<sub>y</sub> × n<sub>z</sub></th><th>Box size<br>L<sub>x</sub> × L<sub>y</sub> × L<sub>z</sub></th><th>Ra</th><th>Ri</th><th>⟨1/2 v<sup>2</sup>⟩/S<sup>2</sup>L<sup>2</sup></th><th>⟨v<sub>x</sub>v<sub>y</sub>⟩/S<sup>2</sup>L<sup>2</sup> ≡ α</th><th>Λ⟨θv<sub>z</sub>⟩/SL<sup>2</sup></th></tr></thead><tbody><tr><td>H1</td><td>64 × 64 × 128</td><td>1 × 1 × 2</td><td>10<sup>4</sup></td><td>0.1</td><td>2.8 × 10<sup>-2</sup></td><td>-3.4 × 10<sup>-4</sup></td><td>1.1 × 10<sup>-1</sup></td></tr><tr><td>H2</td><td>64 × 64 × 128</td><td>1 × 1 × 2</td><td>4 × 10<sup>4</sup></td><td>0.1</td><td>1.8 × 10<sup>-2</sup></td><td>9.6 × 10<sup>-4</sup></td><td>6.2 × 10<sup>-2</sup></td></tr><tr><td>H3</td><td>128 × 128 × 256</td><td>1 × 1 × 2</td><td>4 × 10<sup>5</sup></td><td>0.1</td><td>1.6 × 10<sup>-2</sup></td><td>9.0 × 10<sup>-4</sup></td><td>4.2 × 10<sup>-2</sup></td></tr><tr><td>H4</td><td>128 × 128 × 256</td><td>1 × 1 × 2</td><td>4 × 10<sup>6</sup></td><td>0.1</td><td>1.9 × 10<sup>-2</sup></td><td>1.6 × 10<sup>-3</sup></td><td>3.7 × 10<sup>-2</sup></td></tr><tr><td>H5</td><td>128 × 256 × 256</td><td>1 × 2 × 2</td><td>4 × 10<sup>5</sup></td><td>0.1</td><td>1.5 × 10<sup>-2</sup></td><td>5.3 × 10<sup>-4</sup></td><td>4.0 × 10<sup>-2</sup></td></tr><tr><td>H6</td><td>128 × 256 × 256</td><td>1 × 2 × 2</td><td>4 × 10<sup>6</sup></td><td>0.1</td><td>1.3 × 10<sup>-2</sup></td><td>9.0 × 10<sup>-4</sup></td><td>2.8 × 10<sup>-2</sup></td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Turbulent transport as a function of time for run H2 (dashed line) and H4 (plain line).</figcaption></figure>

larger in homogeneous runs (Table 2) which partly explains the higher turbulent transport observed in those cases.

## 5 DISCUSSION

In this Letter, we have shown that turbulent convection, if present in accretion discs, could lead to an *outward* transport of angular momentum, provided that the Rayleigh number is large enough (i.e. molecular viscosity and thermal diffusivity are small enough). We also conjectured that the previous negative results (Cabot 1996; Stone & Balbus 1996) were due to a too small Rayleigh number, leading to nearly axisymmetric turbulence which cannot transport angular momentum outwards.

Despite this positive result, this Letter does not demonstrate that convection is actually present in discs and is responsible for the observed accretion rates. To demonstrate this point, one probably has to include viscous heating and a realistic model of disc stratification, which are beyond the scope of this Letter.

We note that the use of the Boussinesq approximation with constant stratification is a strong simplification of the problem. However, we always find $v < 0.1SL$, showing that our incompressible approximation is valid. Moreover, we believe that the Rayleigh number effect found in this work will remain valid with a more

realistic vertical profile, as small scales are weakly influenced by the global (vertical) structure.

Interestingly, we find that the turbulent transport of angular momentum is always much weaker than the turbulent heat transport (typically $\langle\theta v_z\rangle \sim 20 \alpha$). This is to be expected as turbulent convection is not shear driven, but it also indicates that convection is inefficient at transporting angular momentum. This leads us to conjecture that it might not be possible to generate vertical convection self-consistently with viscous heating only.

## ACKNOWLEDGMENTS

GL thanks Steven Balbus and John Papaloizou for useful and stimulating discussions. Some of the simulations presented in this Letter were performed using the Darwin Supercomputer of the University of Cambridge. The authors acknowledge support by STFC.

## REFERENCES

Balbus S. A., 2003, ARA&A, 41, 555
Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214
Borie V., Orszag S. A., 1997, J. Sci. Comput., 12, 305
Cabot W., 1996, ApJ, 465, 874
Cattaneo F., Emonet T., Weiss N., 2003, ApJ, 588, 1183
Gammie C. F., 1996, ApJ, 457, 355
Goodman J., Xu G., 1994, ApJ, 432, 213
Hawley J. F., Gammie C. F., Balbus S. A., 1995, ApJ, 440, 742
Klahr H. H., Bodenheimer P., 2003, ApJ, 582, 869
Lesur G., Longaretti P.-Y., 2005, A&A, 444, 25
Lesur G., Longaretti P.-Y., 2007, MNRAS, 378, 1471
Lesur G., Papaloizou J. C. B., 2010, A&A, in press, doi:10.1051/0004-6361/200913594
Lin D. N. C., Papaloizou J., 1980, MNRAS, 191, 37
Lynden-Bell D., Pringle J. E., 1974, MNRAS, 168, 603
Petersen M. R., Stewart G. R., Julien K., 2007, ApJ, 658, 1252
Pope S. B., 2000, Turbulent flows. Cambridge Univ. Press, Cambridge
Regev O., Umurhan O. M., 2008, A&A, 481, 21
Ryu D., Goodman J., 1992, ApJ, 388, 438
Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
Speziale C. G., 1991, Ann. Rev. Fluid Mech., 23, 107
Spiegel E. A., Veronis G., 1960, ApJ, 131, 442
Stone J. M., Balbus S. A., 1996, ApJ, 464, 364

This paper has been typeset from a TeX/LATEX file prepared by the author.