Hydrodynamic convection in accretion discs

# Hydrodynamic convection in accretion discs

Loren E. Held, Henrik N. Latter
Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Centre for Mathematical Sciences,
Wilberforce Road, Cambridge CB3 0WA, UK
E-mail: leh50@cam.ac.uk (LEH)
Accepted XXX. Received YYY; in original form ZZZ
###### Abstract

The prevalence and consequences of convection perpendicular to the plane of accretion discs have been discussed for several decades. Recent simulations combining convection and the magnetorotational instability have given fresh impetus to the debate, as the interplay of the two processes can enhance angular momentum transport, at least in the optically thick outburst stage of dwarf novae. In this paper we seek to isolate and understand the most generic features of disc convection, and so undertake its study in purely hydrodynamical models. First, we investigate the linear phase of the instability, obtaining estimates of the growth rates both semi-analytically, using one-dimensional spectral computations, as well as analytically, using WKBJ methods. Next we perform three-dimensional, vertically stratified, shearing box simulations with the conservative, finite-volume code PLUTO, both with and without explicit diffusion coefficients. We find that hydrodynamic convection can, in general, drive outward angular momentum transport, a result that we confirm with ATHENA, an alternative finite-volume code. Moreover, we establish that the sign of the angular momentum flux is sensitive to the diffusivity of the numerical scheme. Finally, in sustained convection, whereby the system is continuously forced to an unstable state, we observe the formation of various coherent structures, including large-scale and oscillatory convective cells, zonal flows, and small vortices.

###### keywords:
accretion, accretion discs – convection – hydrodynamics – instabilities, turbulence
pubyear: 2018pagerange: Hydrodynamic convection in accretion discsB

## 1 Introduction

Thermal convection – the bulk motion of fluid due to an entropy gradient – transports heat, mass, and angular momentum. There is compelling analytical, numerical, experimental, and observational evidence for convection both in the inner regions of stars and in the outer core and mantle of the Earth (Spruit et al., 1990; Aurnou et al., 2015), but its application to astrophysical discs is less clear. In particular, convection perpendicular to the plane of a disc bears marked differences to convection in most stars and planets. For example, the gravitational acceleration reverses sign at the disc midplane, and the gas supports a strong background shear flow. A third difference concerns the heat source necessary to drive the unstable entropy gradient. In stars this heat is provided by nuclear fusion, while in the Earth’s outer core convection is driven by both thermal and compositional gradients. In discs heating at the mid-plane is normally assumed to be supplied by accretion (and the associated viscous dissipation of heat). Turbulence arising from the magnetorotational instability and dissipation of large-scale density waves are two mechanisms that might convert orbital energy into the required heat. Note importantly that the fluid motions associated with these heating processes may in fact impede the onset of convection, and at the very least interact with it in non-trivial ways.

It has long been speculated that convection perpendicular to the plane of the disc influences the optically thick, high-accreting outburst phase of dwarf novae (for a succinct review, see Cannizzo (2011)). More recently, simulations of magnetorotational turbulence in stratified discs have shown that the MRI is capable of generating a convectively unstable entropy gradient. Moreover, three-dimensional, shearing box simulations with zero-net vertical magnetic flux reveal an interplay between convection and the MRI that might enhance angular momentum transport, typically quantified by the dimensionless parameter (Bodo et al., 2013; Hirose et al., 2014). A similar interplay between magnetorotational instability and convection might operate in the ionized inner regions of protoplanetary discs, or during FU Orionis outbursts (Bell & Lin, 1994; Hirose, 2015).

In addition to dwarf novae and the inner regions of protoplanetary discs, hydrodynamic convection might drive activity in the weakly ionized ‘dead zones’ of protoplanetary discs where the magnetorotational instability is unlikely to operate (Armitage, 2011). Although irradiation by the central star usually results in an isothermal vertical profile in protoplanetary discs (Chiang & Goldreich, 1997; D’Alessio et al., 1998), heating by spiral density waves driven by a planet might nevertheless generate regions exhibiting a convectively unstable background equilibrium (Boley & Durisen, 2006; Lyra et al., 2016). So too might the heating from the dissipation of the strong zonal magnetic fields associated with the Hall effect (Lesur, Kunz and Fromang 2014).

As mentioned, accretion may instigate convection. But convection might drive accretion itself. Motivated by considerations of the primordial solar nebula, Lin & Papaloizou (1980) constructed a simple hydrodynamic model of a cooling young protoplanetary disk contracting towards the mid-plane. They showed that vertical convection arises if the opacity increases sufficiently fast with temperature, with the heating source initially provided by the gravitational contraction. But if convection is modeled via a mixing length theory, with an effective viscosity, the process might self-sustain: convective eddies could extract energy from the background orbital shear, and viscous dissipation of that energy might replace the initial gravitational contraction as a source of heat for maintaining convection.

In order to go beyond mixing length theory, Ruden et al. (1988) analysed linear axisymmetric modes in a thin, polytropic disc in the shearing box approximation and estimated that mixing of gas within convective eddies could result in values of . They cautioned, however, that axisymmetric convective cells could not by themselves exchange angular momentum: such an exchange would have to be facilitated either by non-axisymmetric modes, or by viscous dissipation of axisymmetric modes.

Following these early investigations, a debate ensued that centered not so much on the size of but rather on its sign. Dissipation of non-linear axisymmetric convective cells was investigated by Kley et al. (1993) who performed quasi-global simulations (spanning about 100 stellar radii) of axisymmetric viscous discs with radiative heating and cooling. They measured an inward flux of angular momentum, but warned that this might be due to the imposed axial symmetry and to their relatively high viscosity. Convective shearing waves were first investigated by Ryu & Goodman (1992), who concluded that linear non-axisymmetric perturbations would result in a net inward angular momentum flux at sufficiently large time. These results were questioned however by Lin et al. (1993), who examined analytically and numerically a set of localised linear non-axisymmetric disturbances in global geometry. They demonstrated that these modes could transport angular momentum outwards in some cases, and opined that the inward transport of angular momentum reported by Ryu & Goodman (1992) was an artifact of the shearing box approximation.

Interest in hydro convection waned after local non-linear 3D compressible simulations showed that it resulted in inward rather than outward angular momentum transport. Using the finite-difference code ZEUS and rigid, isothermal vertical boundaries, Stone and Balbus (1996, heareafter SB96) initialized inviscid, fully compressible, and vertically stratified shearing box simulations with a convectively unstable vertical temperature profile, and measured a time-averaged value of . Inward angular momentum transport was also reported by Cabot (1996) who ran simulations similar to those of SB96 but included full radiative transfer and a relatively high explicit viscosity. Analytical arguments for the inward transport of angular momentum by convection were presented by SB96 which, crucially, assumed axisymmetry, especially in the pressure field. However, in a rarely cited paper Klahr et al. (1999) presented fully compressible, three-dimensional, global simulations including explicit viscosity and radiative transfer of hydrodynamic convection in discs that gave some indication that non-linear convection in discs actually assumes non-axisymmetric patterns.

Some fifteen years later, the claims of SB96 and of Cabot (1996) were called into question: fully local shearing box simulations of Boussinesq hydrodynamic convection in discs using the spectral code SNOOPY and employing explicit diffusion coefficients indicated that the sign of angular momentum transport due to vertical convection can be reversed provided the Rayleigh number (the ratio of buoyancy to kinematic viscosity and thermal diffusivity) is sufficiently large. It would appear then that vertical convection can possibly drive outward angular momentum transport after all (Lesur & Ogilvie, 2010).

Our aim in this paper is to explore the competing results and claims concerning hydrodynamic convection in accretion discs, in particular to determine the sign and magnitude of angular momentum transport. We also isolate and characterize other generic features of convection that should be shared by multiple disc classes (dwarf novae, protoplanetary disks, etc) and by different driving mechanisms (MRI turbulence, spiral shock heating, etc). We do so both analytically and through numerical simulations, working in the fully compressible, vertically stratified shearing box approximation. We restrict ourselves to an idealised hydrodynamic set-up, omitting magnetic fields and complicated opacity transitions, and in so doing ensure our results are as general as possible.

The structure of the paper is as follows: first, in Section 2 we introduce the basic equations and provide a brief overview of the numerical code, the numerical parameters and set-up, and our main diagnostics. In Section 3 we investigate the linear behavior of the convective instability, employing both analytical WKBJ and semi-analytical spectral methods to calculate the growth rates and the eigenfunctions. In Section 4 we explore the non-linear regime through unforced simulations, in which the convection is not sustained and is permitted to decay after non-linear saturation. Because this unforced convection is a transient phenomenon which might depend on the initial conditions, in Section 5 we explore the non-linear regime through simulations of forced convection, in which the internal energy is relaxed to its initial, convectively unstable, state to mimic, possibly, the action of background MRI turbulent dissipation and optically thick radiative cooling. Finally, Section 6 summarizes and discusses the results.

## 2 Methods

### 2.1 Governing equations

We work in the shearing box approximation (Goldreich & Lynden-Bell, 1965; Hawley et al., 1995; Latter & Papaloizou, 2017), which treats a local region of a disc as a Cartesian box located at some fiducial radius and orbiting with the angular frequency of the disc at that radius . A point in the box has Cartesian coordinates which are related to the cylindrical coordinates through , and . The equations of gas dynamics are now

 ∂tρ+∇⋅(ρu)=0, (1)
 ∂tu+u⋅∇u=−1ρ∇P−2Ω0ez×u+geff+1ρ∇⋅T, (2)
 ∂t(ρe)+u⋅∇(ρe)=−γρe∇⋅u+T:∇u+ξ∇2T+Λrelax, (3)

corresponding to conservation of mass, momentum and thermal energy, respectively. These equations are closed with a perfect gas equation of state , or, equivalently, in terms of the specific internal energy . Here is the thermal pressure of the fluid, is the temperature, is mass density, is the gas constant, and is the mean molecular weight. The adiabatic index (ratio of specific heats) is denoted by and is typically taken to be in all our simulations.111Technically is valid for ionized hydrogen. Since all our simulations are hydrodynamic it might be more appropriate to use , which is valid for molecular hydrogen, but the differences are small and do not affect our results. The effective gravitational potential is embodied in the tidal acceleration , where is the shear rate. For Keplerian disks a value we adopt throughout the paper.

Other terms in the equations include the viscous stress tensor . Here is the kinematic viscosity, and is the traceless shear tensor. The thermal conductivity is , but in our simulations we specify thermal diffusivity rather than . The former is related to the latter via where is the specific heat capacity at constant pressure.

In some of our simulations we include heating and cooling through a thermal relaxation term . This relaxes the thermal energy in each cell back to its initial state on a time-scale given by . We refer to simulations run with as simulations of forced compressible convection. Conversely, simulations with are referred to as unforced. The thermal relaxation time-scale is taken to be the inverse of the linear growth rate of the convective instability.

### 2.2 Important parameters and instability criteria

The onset of thermal convection is controlled by the local buoyancy frequency, defined as

 N2B=z[1γ∂lnP∂z−∂lnρ∂z]Ω20. (4)

But convection is opposed by the effects of viscosity and thermal diffusivity, with the ratio of destabilising and stabilising processes quantified by the Rayleigh number

 Ra≡N2BH4νχ, (5)

where is the disk’s scaleheight (formally defined later). On the other hand, the ratio of buoyancy to shear is expressed through the Richardson number

 Ri≡N2Bq2Ω20. (6)

But for convection perpendicular to the plane of an accretion disc, the buoyancy force is perpendicular to the direction of the background shear. Thus the Richardson number in our set-up is more an expression for the scaled intensity of . Finally, the Prandtl number is defined as the ratio of kinematic viscosity to thermal diffusivity

 Pr≡νχ. (7)

When the explicit diffusion coefficients are neglected (), convective instability occurs when the buoyancy frequency is negative (i.e. in regions where ). When explicit viscosity and thermal diffusivity are included, however, convective instability requires both that , and that the Rayleigh number exceed some critical value . Though in real astrophysical disks the viscosity is negligible, the inclusion of magnetic fields complicates the stability criterion, as does the presence of pre-existing turbulence (as might be supplied by the MRI) which itself may diffuse momentum and heat. In the latter case it may be possible to define a turbulent ‘Rayleigh number’, which must be sufficiently large so that convection resists the disordered background flow. In any case, the sign of is certainly insufficient to assign convection to MRI-turbulent flows, as is often done in recent work.

### 2.3 Numerical set-up

#### 2.3.1 Codes

Most of our simulations are run with the conservative, finite-volume code PLUTO (Mignone et al., 2007). Unless stated otherwise, we employ the Roe Riemann-solver, 3rd order in space WENO interpolation, and the 3rd order in space Runge-Kutta algorithm. Other configurations are explored in Section 4. To allow for longer time-steps, we employ the FARGO scheme (Mignone et al., 2012). When explicit viscosity and thermal diffusivity are included, we further reduce the computational time by using the Super-Time-Stepping (STS) scheme (Alexiades et al., 1996). Note that due to the code’s conservative form kinetic energy is not lost to the grid but converted directly into thermal energy. Ghost zones are used to implement the boundary conditions.

We use the built-in shearing box module in PLUTO (Mignone et al., 2012). Rather than solving Equations (1)-(3) (primitive form), PLUTO solves the governing equations in conservative form, evolving the total energy equation rather than the thermal energy equation, where the total energy density of the fluid is given by (kinetic + internal). In PLUTO the thermal relaxation term is not implemented directly in the total energy equation. Instead, it is included on the right-hand-side of the thermal energy equation, which is then integrated (in time) analytically.

Finally, we have verified the results of our fiducial simulation (presented in Section 4.1) with the conservative finite-volume code ATHENA (Stone et al., 2008).

#### 2.3.2 Initial conditions and units

We use two different convectively unstable profiles to initialize our simulations. Nearly all employ an equilibrium exhibiting a Gaussian temperature profile

 T=T0exp[−z2βH2], (8)

where is midplane temperature and is a dimensionless tuning parameter. See Figure 1 for all associated thermal profiles. For comparison we also use the profile introduced by SB96 in which the temperature follows a power law

 T=T0−Azp, (9)

where and are parameters, with usually (see Figure 2). Further details of both profiles are given in Appendix A. The equilibria are convectively unstable within a confined region (of size ) about the disc mid-plane, and convectively stable outside this region. We have checked that these vertical profiles are indeed good equilibrium solutions in our numerical set-up by running a series of simulations without perturbations; these show that after 55 orbits the initial profiles are unchanged with velocity fluctuation amplitudes typically less than at the vertical boundaries and less than at the midplane.

The background velocity is given by . At initialization we usually perturb all the velocity components with random noise exhibiting a flat power spectrum. The perturbations have maximum relative amplitude of about and can be either positive or negative. In order to investigate specifically the nature of linear axisymmetric convective modes, we initialize several PLUTO simulations with linear axisymmetric modes calculated semi-analytically rather than with random noise (see Section 3.5).

Time units are selected so that . From now the subscript on the angular frequency is dropped if it appears. The length unit is chosen so that the initial midplane sound speed , which in turn defines a constant reference scaleheight . Note, however, that the sound speed is generally a function of both space and time.

#### 2.3.3 Box size and resolution

We measure box size in units of scaleheight , defined above. We employ resolutions of , and in boxes of size , which correspond to resolutions of 16, 32, 64 and 128 grid cells per scaleheight in all directions. However, in simulations in which we force convection we employ a ‘wide-box’ () with a resolution of which corresponds to about grid cells per in the horizontal directions, and 64 grid cells per in the vertical direction.

#### 2.3.4 Boundary conditions

We use shear-periodic boundary conditions in the -direction (see Hawley et al. (1995)) and periodic boundary conditions in the -direction. In the vertical direction, we keep the ghost zones associated with the thermal variables in isothermal hydrostatic equilibrium (the temperature of the ghost zones being kept equal to the temperature of the vertical boundaries at initialization) in the manner described in Zingale et al. (2002). For the velocity components we use mostly standard outflow boundary conditions in the vertical direction, whereby the vertical gradients of all velocity components are zero, and variables in the ghost zones are set equal to those in the active cells bordering the ghost zones. In a handful of our simulations we also use free-slip and periodic boundary conditions to test the robustness of our results.

### 2.4 Diagnostics

#### 2.4.1 Averaged quantities

We follow the time-evolution of various volume-averaged quantities (e.g. kinetic energy density, thermal energy density, Reynolds stress, mass density, thermal pressure, and temperature). For a quantity the volume-average of that quantity is denoted and is defined as

 ⟨X⟩(t)≡1V∫VX(x,y,z,t)dV (10)

where is the volume of the box.

If we are interested only in the vertical structure of a quantity then we average over the - and -directions, only. The horizontal average of that quantity is denoted and is defined as

 ⟨X⟩xy(z,t)≡1A∫AX(x,y,z,t)dA, (11)

where is the horizontal area of the box. Horizontal averages over different coordinate directions (e.g. over the - and -directions) are defined in a similar manner.

We are also interested in averaging certain quantities (e.g. the Reynolds stress) over time. The temporal average of a quantity is denoted and is defined as

 ⟨X⟩t(x,y,z)≡1Δt∫tftiX(x,y,z,t)dt, (12)

where we integrate from some initial time to some final time and .

#### 2.4.2 Reynolds stress, alpha, and energy densities

In accretion discs, the radial transport of angular momentum is related to the -component of the Reynolds stress, defined as

 Rxy≡ρuxδuy (13)

where is the fluctuating part of the y-component of the total velocity . corresponds to positive (i.e. radially outward) angular momentum transport, while corresponds to negative (i.e. radially inward) angular momentum transport.

The Reynolds stress is related to the classic dimensionless parameter , but the reader should note that various definitions of exist in the literature and care must be taken when comparing measurements of quoted in different papers. We define as

 α≡Ω0⟨Rxy⟩qΩ⟨P⟩. (14)

The kinetic energy density of a fluid element is defined as

 Ekin≡12ρu2, (15)

where is the magnitude of the total velocity of a fluid element. Often we will plot the vertical kinetic energy density, in which case in the above is replaced by .

The total energy (kinetic plus thermal) is not conserved in the shearing box even when there are no (explicit or numerical) dissipative effects: work done by the effective gravitational potential energy on the fluid and flux of angular momentum through the -boundaries due to the Reynolds stress act as sources for the total energy. The upper and lower vertical boundaries also let energy (and other quantities) leave the box.

#### 2.4.3 Vertical profiles of horizontally averaged quantities

Finally, we track the vertical profiles of horizontally averaged pressure, density, temperature, and Reynolds stress. Of special interest is the buoyancy frequency, which is calculated from the pressure and density data by finite differencing the formula

 N2BΩ2=z[1γdln⟨P⟩xydz−dln⟨ρ⟩xydz]. (16)

We also calculate the (horizontally averaged) vertical profiles of mass and heat flux. We define the mass flux as

 Fmass=⟨ρuz⟩xy, (17)

and the heat flux as

 Fheat=⟨ρuzT⟩xy. (18)

## 3 Linear theory

In this section we investigate the linear phase of the axisymmetric convective instability both semi-analytically by solving a 1D boundary value / eigenvalue problem using spectral methods, and also analytically using WKBJ methods.

Our aim is to determine the eigenvalues (growth rates) and eigenfunctions (density, velocity and pressure perturbations) of the axisymmetric convective instability in the shearing box as a function of radial wavenumber . The convectively unstable background vertical profile used in all calculations is shown in Figure 1. The profile is discussed in detail in Section A.2. All calculations were carried out with profile parameters and an adiabatic index of . For our chosen background equilibrium, this corresponds to a Richardson number of .

We proceed as follows: first, we produce linearized equations for the perturbed variables; second, these are Fourier analyzed so that for each of the perturbed fluid variables and . Here can be complex (if is real and positive it is referred to as a growth rate). The linearized equations are given by

 σρ′=−ikxρ0u′x−∂z(ρ0u′z), (19)
 σu′x=2Ωu′y−ikxρ0P′, (20)
 σu′y=(q−2)Ωu′x, (21)
 σu′z=−1ρ0∂zP′+ρ′ρ20∂zP0, (22)
 σP′=−u′z∂zP0−γP0(ikxu′x+∂zu′z). (23)

Background fluid quantities are denoted . For equilibrium and see Appendix A.

### 3.1 Spectral methods

We discretize the fluid variables on a vertical grid containing points by expanding each eigenfunction as a linear superposition of Whittaker cardinal (sinc) functions. These are approprate as all perturbations should decay to zero far from the convectively unstable region. This choice also saves us from explicitly imposing boundary conditions. The discretised equations are next gathered up into a single algebraic eigenvalue problem (Boyd, 2001). Finally, we solve this matrix equation numerically using the QR algorithm to obtain the growth rates and eigenfunctions for a given radial wavenumber . The code which computes the eigensolution we often refer to as our ‘spectral eigensolver’.

### 3.2 WKBJ approach

Alternatively, Equations (19)-(23) can be combined into a second-order ODE of the form

 d2Udz2+k2z(z)U=0, (24)

where

 U=(γP0σ2+κ2+k2xc20)1/2u′z, (25)

in which is the equilibrium sound speed, and the ‘vertical wavenumber’ of the perturbations is given by

 k2z(z)≈−k2xN2B(z)+σ2κ2+σ2, (26)

with the epicyclic frequency (Ruden et al., 1988).

We obtain approximate solutions to Equation (24) analytically using WKBJ methods in the limit large. Equation (24) has two turning points at , which occur when . Physically, one turning point occurs at the disc mid-plane (), while the other turning point occurs at the boundary of the convectively unstable region . When then so the solutions of Equation (24) are spatially oscillatory, otherwise they are evanescent. The unstable modes, in which we are interested, are confined and oscillatory within the convectively unstable region, and spatially decay exponentially outside.

The standard WKBJ solution is given by

 U(z)∼c1exp(+i∫zkz(z′)dz′)+c2exp(−i∫zkz(z′)dz′),

where and are constants. By matching the interior (oscillatory) solution onto the exterior (exponentially decaying) solution, we obtain the eigenvalue equation (dispersion relation):

 ∫z1z0kz(z)dz=(n+12)πn=0,1,2,3,⋯, (27)

where lower and upper bounds of integration and , respectively, are related by the implicit equation , and is the mode number. Substituting for using Equation (26), we obtain an implicit equation relating the radial wavenumber and the growth rate which we solve numerically via a root-finding algorithm.

### 3.3 Eigenvalues

For a given radial wavenumber , the eigenvalues calculated with the spectral eigensolver come in pairs that lie close together in the complex plane. These correspond physically to the even and odd modes of the convective instability. As described in Section 3.4, even modes are symmetric about the mid-plane, whereas the odd modes are antisymmetric about the mid-plane. In Fig. 3 we plot as black dashed curves the first few eigenvalue branches as functions of , these corresponding to vertical quantum numbers and . Note that in the figure has been normalized by the width of the convectively unstable region (given by Equation (32) in Appendix A.2).

As the radial wavenumber is increased, the difference between the growth rates of the even and odd modes vanishes, and the growth rates tend asymptotically to the maximum absolute value of the buoyancy frequency (indicated by a horizontal dot-dashed line). In other words, the maximum growth rate is limited by the depth of the buoyancy-frequency profile. The fastest growing modes are at arbitrarily large (small radial wavenumber ), and thus manifest as thin elongated structures. Thicker structures are not favoured as they comprise greater radial displacements that are resisted by the radial angular momentum gradient.

Superimposed on the spectrally computed values are the WKBJ growth rates (solid lines). WKBJ methods cannot distinguish between even or odd modes, nonetheless we find very close agreement between the WKBJ and semi-analytical results for all radial wavenumbers, even at low where the WKBJ approximation is, strictly speaking, invalid. As the radial wavenumber increases, the semi-analytical and WKBJ results converge. Note that associated with each mode branch is a minimum radial wavenumber below which both the numerical and WKBJ growth rates are zero. This feature is just visible in the bottom left-hand side of Figure 3.

We must stress that, because the fastest growing modes are on the shortest possible scales, inviscid simulations of convection are problematic, at least in the early (linear) phase of the evolution. It is here that the simulated behaviour may depend on the varying ways that different numerical schemes deal with grid diffusion, something we explore in Section 4.3. Only with resolved physical viscosity can this problem be overcome. However, in the fully developed nonlinear phase of convection the short scale linear axisymmetric modes may be subordinate and this is less of an issue.

### 3.4 Vertical structure of the eigenfunctions

The reader is reminded that here we have an unstably stratified fluid in a medium in which the gravitational field smoothly changes sign. This contrasts to the more commonly studied situation of convection at the surface of the Earth (or within the solar interior). It hence is of interest to study the structure of the convective cells in some detail. We consider first the fundamental () odd mode with radial wavenumber .

In Figure 4 we plot the vertical profiles of the nonzero components of the eigenfunction. To provide greater visual clarity, we have rescaled the perturbations so that the maximum amplitude of each perturbation is unity: thus we can more easily observe the vertical structure of the perturbations, but not deduce their relative amplitudes. In Figure 5 we plot their real parts in the -plane with the correct amplitudes. We find that . Thus vertical velocity perturbations are greatest in magnitude, while pressure perturbations are smaller than vertical velocity perturbations by two-orders of magnitude, indicating that the convective cells are roughly in pressure balance with the surroundings, i.e that they are behaving adiabatically.

In Figure 4, (solid black curve) is antisymmetric about the mid-plane, whereas the remaining perturbations are symmetric about the mid-plane. Figs 4 and 5 clearly show that both odd and even modes consist of a chain of convective cells above and below the midplane and localized near the most convectively unstable point (denoted by the large red dot in the former Figure). The peak in the amplitude of that occurs near the this point tells us that the fluid elements reach their maximum acceleration where the buoyancy force is greatest. As the fluid perturbations behave adiabatically, they adjust their pressure to maintain a balance with the background pressure: thus where cool elements begin to rise (higher background pressure), the pressure perturbations (solid blue line) increase (i.e. ), and where hot elements begin to sink (lower background pressure), the pressure perturbations decrease (i.e. ). The vertical velocity perturbation (black solid line) is out of phase with the radial and azimuthal velocity perturbations (dashed lines), since vertical motion is converted into radial motion where the convective cells turn over (picture fluid motion at the top of a fountain).

The even modes possess a nonzero, but relatively small, at the midplane, and thus weakly couple the two sides of the disk. Such modes may permit some exchange of mass, momentum and thermal energy across the mid-plane when reaching large amplitudes.

As the radial wavenumber is increased from to (very small radial wavelengths), the perturbations become increasingly localized about the most convectively unstable point(s). The localisation in however is not as narrow as the radial wavelength. For larger , in fact, each pair of even and odd modes change their character and become entirely localised to either the upper or lower disk. Activity in the two halves of the disk are hence completely decoupled for such small-scale (but fast growing) disturbances.

### 3.5 Comparison of theory with simulations

In this section we compare the growth rates previously calculated with those measured from PLUTO simulations initialized with the exact linear modes taken from the spectral eigensolver.

Each PLUTO simulation is run in a box of size . We maintain a fixed number of grid cells per scaleheight in the -direction, grid cells per scaleheight in the -direction and grid cells per radial wavelength in the -direction. Thus the radial wavelength of the modes is always well resolved in our simulations. No explicit diffusion coefficients (i.e. viscosity or thermal diffusivity) are used.

The results are plotted in Figure 3 as squares (and one circle). Numerical values appear in Table 1. We find excellent agreement between simulation and theory: the relative error is for the fundamental mode in the interval and for the first harmonic at . As the radial wavenumber is increased the simulation growth rates begin to deviate from the theoretical curves, although the percentage error remains less than even at large wavenumbers. This behavior is expected due to the effects of numerical diffusion which acts to decrease the growth rates as the radial wavelength approaches the size of the vertical grid. Although the radial wavenumber is always well resolved in the simulations, the vertical wavelength becomes increasingly less well resolved.

The numerical diffusion could be reduced by increasing the vertical resolution of the simulations. To check this we have rerun the PLUTO simulation initialized with eigenfunctions corresponding to the first harmonic even mode with at twice the vertical resolution, i.e. instead of . At a resolution of we measured a growth rate of , corresponding to a percentage error of with the growth rate calculated semi-analytically (. At a vertical resolution of , however, we measured a growth rate of , corresponding to a percentage error of just . Thus, the growth rates measured in the simulations appear to converge to those calculated analytically and semi-analytically as the resolution of the simulations is increased.

Finally, in all simulations of the axisymmetric modes we observed inward angular momentum transport. This is in agreement with the analytical argument presented in (Stone & Balbus, 1996) regarding axisymmetric flow.

## 4 Simulations of Unforced Compressible Convection

In this section we describe fully compressible, three-dimensional simulations of hydrodynamic convection in the shearing box carried out with PLUTO. These simulations are ‘unforced’ in that they are initialized with a convectively unstable profile, but this unstable profile is not maintained by a separate process (e.g. turbulent heating via the magnetorotational instability). Thus convection rearranges the profile, shifting it to a marginally stable state, and thus ultimately hydrodynamical activity dies out. While the unforced convection studied in this section is essentially a transient phenomenon that depends on initial conditions, and is sensitive to the numerical scheme (as we discuss in Sections 4.2 and 4.3), it serves as a good starting point for investigating the problem and illuminates a number of interesting features.

The key result of this section is that, even without the inclusion of explicit viscosity, we observe that hydrodynamic convection in the shearing box generally produces a positive Reynolds stress, and thus can drive outward transport of angular momentum. This is in direct contradiction to simulations carried out in ZEUS by SB96. In Section 4.1 we describe our fiducial simulations, and in Section 4.2 we compare ATHENA and PLUTO simulations. In Section 4.3 we describe the sensitivity of the sign of angular momentum transport to the numerical scheme. Finally, in Section 4.4 we investigate the effects of including explicit diffusion coefficients in our simulations, thus connecting to the Boussinesq simulations of Lesur & Ogilvie (2010).

### 4.1 Fiducial simulations

#### 4.1.1 Set-up and initialization of fiducial simulations

The simulations described in this section were run at resolutions of , and in boxes of size , where is the scaleheight.222A simulation run at a resolution is included in the list of fiducial simulations in Table 2 (see Appendix B) but is not described in this section for brevity. All the simulations were initialized with the convectively unstable vertical profiles for density and pressure described in Appendix A.2 with profile parameters and an adiabatic index of (see Figure 1). The Stone and Balbus vertical profile SB96 was also trialled yielding very similar results, but we have omitted most of these in the interests of space. Vertical outflow conditions were employed for our fiducial simulations but periodic and free-slip boundary conditions in the vertical direction produced the same behaviour. Random perturbations to the velocity components were seeded at initialization, with a maximum amplitude . Finally, we adopt a very small but finite thermal diffusivity of to facilitate conduction of thermal energy through the vertical boundaries and to aid code stability. No explicit non-adiabatic heating, cooling, or thermal relaxation are included in the simulations described in this section: therefore convection gradual dies away after non-linear saturation.

#### 4.1.2 Time-evolution of averaged quantities

In the left panel of Figure 6 we track the time-evolution of the volume-averaged kinetic energy density associated with the vertical velocity for simulations run at resolutions of and , respectively.

Initially all simulations exhibit small-amplitude oscillations due to internal gravity waves excited in the convectively stable region at initialization. After some three orbits, the linear phase of the convective instability begins in earnest, characterised by exponential growth of the perturbation amplitudes. During this phase, internal energy is converted into the kinetic energy of the rising and sinking fluid motion that comprises the convective cells. As the resolution is increased, shorter scale modes are permitted to grow. Because they are the most vigorous the growth rate (proportional to the slope of the kinetic energy density) is slightly larger at better resolution.

The linear phase ends some 10 orbits into the simulation, and the flow becomes especially disordered. The peak in the vertical kinetic energy density occurs at this point, but this peak decreases with resolution, something we discuss below. After non-linear saturation, the kinetic energy decreases gradually: the convective cells redistribute thermal energy and mass, thus shifting the thermal profile from a convectively unstable to a marginal state, and ultimately the convective motions die out. After about 38 orbits, the kinetic energy density has dropped to about one hundredth of its value at non-linear saturation in the lower resolution runs. The level of numerical diffusivity has an appreciable effect in damping activity after non-linear saturation, with the decrease in vertical kinetic energy successively smaller over the same period of time in the , and simulations, respectively.

The behaviour of the kinetic energy aligns relatively well with our expectations. The Reynolds stress, on the other hand, is more interesting. The evolution of the -component of the Reynolds stress is plotted in the right panel of Fig. 6. The stress is small, but perhaps its most striking feature is its positivity over the entire duration of all three simulations.

Perhaps most surprising is the positivity of the stress during the linear phase of the instability. We might expect that the axisymmetric modes (which send angular momentum inwards) dominate this period of the evolution, as non-axisymmetric disturbances only have a finite window of growth (about an orbit) before they are sheared out and dissipated by the grid. But the positivity of the stress suggests that instead it is the shearing waves that are the dominant players in the linear phase. Visual inspection of the velocity fields confirms that the flow is significantly non-axisymmetric, and we find several examples of strong shearing waves ‘shearing through’ during this phase. It would appear these waves transport angular momentum primarily outward as they evolve from leading to trailing, behaviour that is in fact consistent with Ryu & Goodman (1992), who find that inward transport only occurs at sufficiently long times when the waves are strongly trailing and hence effectively axisymmetric. (Even then the stress is extremely oscillatory.)

Why do the shearing waves dominate over the axisymmetric modes? One hypothesis is that 3D white noise seeds fast growing shearing waves that can outcompete the axisymmetric convective instability over their permitted window of growth. If this were to be true then the simulated nonlinear regime is achieved before the dominant shearing waves become strongly trailing and begin to send angular momentum inward. Given that the typical timescale for shearing out is roughly a few orbits at best, this is marginal but not impossible. An alternative, more plausible, and more troubling hypothesis is that our inviscid numerical code misrepresents trailing shearing waves: more specifically, the code is artifically reseeding fresh leading waves from strongly trailing waves (‘aliasing’; Geoffroy Lesur, private communication). As a consequence, shearing waves outcompete the axisymmetric modes because they can shear through several times. If true, this is certainly concerning. But we hasten to add that this phenomena should only be problematic in the low amplitude linear phase; once the perturbations achieve large amplitudes the aliasing will be subsumed under physical mode-mode interactions.

The linear phase ends in a spike in the stress. Time-averages of the volume-averaged value of from orbit 5 to orbit 15, a period spanning non-linear saturation, show that the Reynolds stress decreases as the resolution increases (see Table 2 in Appendix B). The time- and volume-averaged value of over a period spanning non-linear saturation (orbits 5 to 15) is , and in the , and simulations, respectively. In terms of the turbulent alpha-viscosity this corresponds (from Equation (14)) to , and , respectively. The dependence on the numerical viscosity can be explained by appealing to secondary shear instabilities (see next subsection).

The volume-averaged density remains roughly constant during the linear phase of the instability, followed by a drop after non-linear saturation as mass is lost through the lower and upper boundaries. Although the convectively unstable region at initialization is confined to the region , and is therefore well within the box, convective overshoot deposits mass (and heat) outside the convectively unstable region. The overall decrease in density over the duration of the simulation is small, but appears to increase with greater resolution. The overall percentage change in the density is , and in the , and simulations, respectively.

#### 4.1.3 Structure of the flow

The development of convective instability and the associated convective cells are best observed through snapshots of the vertical component of the velocity in the -plane, shown at resolutions of and just after non-linear saturation in Figure 7. A full set of convective cells is clearly visible at all resolutions. They are thin, filamentary structures several grid cells wide comprising alternating negative and positive velocities (updrafts and downdrafts).

The maximum vertical Mach number of the flow around non-linear saturation in the run is about , with the largest vertical Mach numbers generally being measured near the vertical boundaries (where the temperature – and therefore the sound speed – is lowest).

The higher resolution simulations indicate that the plumes develop a wavy or buckling structure as they rise or sink, indicating the onset of a secondary shear instability. It is likely that the buckling of the convective plumes by these ‘parasitic modes’ limits the amplitudes of the linear modes, and ultimately leads to their breakdown. At lower resolution numerical diffusion smooths out the secondary shear modes, and so the convective plumes reach large amplitudes before breaking down (blue curve in right panel of Fig. 6). At high resolution the shear modes are not so impeded and make short work of these structures (black curve in the same figure).

#### 4.1.4 Vertical heat and mass flux

Figure 8 (a) shows the vertical profiles of horizontally-averaged heat and mass fluxes taken from snapshots just after non-linear saturation from the simulation with resolution . For clarity, we have time-averaged the horizontally-averaged vertical mass and heat flux profiles between orbits 5 and 10, a period spanning non-linear saturation (see Figure 6).

Negative (positive) values for the fluxes for and positive (negative) values for the fluxes for correspond to transport of heat and mass away from (towards) the mid-plane. Overall, the heat flux is away from the mid-plane, peaking in the vicinity of the most convectively unstable points at initialization (indicated by the vertical dashed red lines in Figure 8). In addition, there is some mass flux towards the mid-plane within . Thus, convection is transporting mass and heat such as to erase the convectively unstable stratification, as expected. Note that the positive peaks in the heat flux near the most convectively unstable points are similar to those observed by Hirose et al. (2014) in the thermally dominated region of the disc (i.e. ; c.f. the middle-panel of Figure 5 in their paper).

Beyond , the mass flux is directed away from the mid-plane, peaking just beyond the point which marks the boundary of the convectively unstable region. This outward mass flux might be due to convective overshoot, although we expect this effect to vanish if averaged over a suitable time interval. Alternatively, it is possible that heat transported towards the corona by the convective cells causes matter in the corona to heat up and become buoyant, or else that this heat generates a weak thermal wind.

Finally, in Figure 8 (b) we show the vertical profile of the horizontally-averaged taken from the simulation and time-averaged between orbits 5 and 10. The Reynolds stress is clearly positive over all of the vertical domain, peaking just beyond the most convectively unstable points.333Note however that at any given instant in time, we generally do observe regions where the Reynolds stress is negative. Thus the bulk of outward angular momentum transport occurs where the convective cells begin to turn-over, resulting in radial mixing of the gas. Note, however, that we also observe rather large positive stresses near the vertical boundaries, with the stress at the vertical boundaries about 1.5 times the peak stress in the remainder of the domain.

### 4.2 Simulation of compressible hydro convection in Athena.

In the previous section we found that hydrodynamic convection in a vertically stratified shearing box in PLUTO (without explicit viscosity) can drive outward angular momentum transport. Here we verify this result using the finite-volume code ATHENA (Stone et al., 2008; Stone & Gardiner, 2010). To facilitate as close a comparison as possible between the two codes and also to the ZEUS runs of SB96, we initialize both codes with Stone and Balbus’s convectively unstable profile (see Appendix A.1). Explicit diffusion coefficients were omitted and vertical periodic boundary conditions implemented. Both simulations were run in a shearing box of size and at a resolution of . PLUTO and ATHENA offer somewhat different suites of numerical schemes: here we settle on a combination that is slightly more diffusive than that employed in our fiducial simulations because this combination allows for as close a match as possible. Specifically, the numerical scheme employed in ATHENA is (second-order) piecewise linear interpolation on primitive variables, the HLLC Riemann solver and MUSCL-Hancock integration. In PLUTO we use (second-order) piecewise linear interpolation on primitive variables together with a Van Leer limiter function, the HLLC Riemann solver and MUSCL-Hancock integration. The angular frequency, and the sound speed at initialization, in both simulations were set to and , respectively.

We find that angular momentum transport is directed outwards in both codes, demonstrating that the outward transport of angular momentum by hydrodynamic convection in the non-linear phase is robust to a change of code. Figure 9 compares the time-evolution of the volume-averaged Reynolds stress taken from the ATHENA simulation with that taken from the PLUTO simulation. Both simulations exhibit exponential growth followed by non-linear saturation, together with the development of convective cells (not shown). These results also contrast with the SB96 runs with ZEUS and thus demonstrate that the positive transport reported in previous sections is not special to the Gaussian temperature profile.

Although for this particular combination of schemes the overall result is the same, we noticed that different schemes resulted in qualitative differences between the two codes. For example, when we use the HLLC solver, piecewise parabolic reconstruction (PPM) and Corner-Transport-Upwind (CTU) integration, ATHENA exhibits delayed onset of instability, a more gradual drop in kinetic energy density following non-linear saturation, and a slower drop in angular momentum transport compared to PLUTO. When we use the Roe solver, PPM reconstruction and CTU integration, the Reynolds stress is highly oscillatory in time in both codes, indicative, perhaps, of numerical instability. The different behavior (in both codes) based on which combination of numerical schemes is chosen is worrying, although we emphasize that the differences are probably accentuated by the transient nature of the evolution and its sensitivity to the initial conditions, and the fact that the linear phase is partially controlled by grid diffusion. The robustness of the results to changes is numerical scheme is investigated in more detail in the following section.

### 4.3 Sensitivity of sign of angular momentum transport to numerical scheme

An important result of Sections 4.1-4.2 is that purely hydrodynamic convection in PLUTO and in ATHENA resulted in , i.e. outward angular momentum transport. This is in disagreement with the ZEUS results of SB96, who reported a Reynolds stress of . Given that ZEUS is a non-conservative, finite-difference code, our hypothesis for explaining the discrepancy is that ZEUS run at comparatively low resolution (as in SB96) is sufficiently diffusive that it imposes an artificial near-axisymmetry on the flow (which will send angular momentum inward). In this section we report our attempts to assess the numerical diffusivity of various algorithms and their impact on convection.

First we reran the fiducial simulations described in Section 4.1.1 but with different combinations of numerical schemes. We found that the sign of angular momentum transport is indeed sensitive to our choice of scheme. This is best illustrated in Figure 10: the left-hand panel shows the time-evolution of the volume-averaged taken from a simulation run with the HLLC Riemann solver, while the right-hand panel shows the same quantity but taken from a simulation run with the (more diffusive) HLL solver. Both simulations exhibit similar exponential growth in the vertical kinetic energy during the linear phase of the instability, and the velocity field shows the development of convective cells in both cases, but the sign of the Reynolds stress is radically different. The HLL run is also far more laminar and axisymmetric. We next repeated the HLL runs with higher resolutions, up to , but with no change in the sign of . It must be stressed that going to higher resolutions does not necessarily help in the problem of convection; this is because the fastest growing linear modes are always near the grid scale and hence it is impossible (in the linear phase at least) to escape grid effects.

The results of different combinations of schemes is summarized in Table 3. They indicate that the sign of appears to be robust to changes in the interpolation and time-stepping schemes, but sensitive to the Riemann solver. In particular, the less diffusive Riemann solvers (Roe and HLLC) gave , while the more diffusive Riemann solvers (HLL and a simple Lax-Friedrichs solver) gave . Altogether we explored twelve different configurations of Riemann solver, interpolation scheme, and time-stepping method. These ranged from the most accurate (least diffusive) set-up which was identical to that employed in the simulations described in the previous section (a Roe Riemann solver, third-order-in-space WENO interpolation and third-order-in-time Runge-Kutta time-stepping), to the least accurate (most diffusive) set-up (a simple Lax-Friedrich’s Riemann solver, second-order-in-space linear interpolation and second-order-in-time Runge-Kutta time-stepping).

Thus we have compelling evidence that, as suspected, angular momentum transport due to convection can be sensitive to the diffusivity of the underlying numerical scheme. It appears that over-smoothing of the flow by diffusive Riemann solvers, such as HLL, or by codes such as ZEUS that employ artificial viscosity and finite-differencing of the pressure terms, impose a spurious axisymmetry on the flow. This axisymmetry picks out the axisymmetric convective modes, which in turn transport angular momentum inwards.

### 4.4 Viscous simulations

Given the concerns raised in the last section regarding numerical schemes, as well as the fact that the fastest growing inviscid modes are on the shortest scales, we expand our study to include explicit viscosity (and thermal diffusivity). A properly resolved viscous simulation should exhibit none of the numerical problems encountered above, and the fastest growing mode occurs on a well defined scale above the (resolved) viscous scale. Our main aim in this section is to test whether the results of our fiducial simulations are solid: mainly, if angular momentum transport can be positive in the presence of viscosity. Additionally, the inclusion of explicit diffusion coefficients enables us to investigate the Rayleigh number dependence of fully compressible hydrodynamic convection in the shearing box, and thus connect to previous work by Lesur and Ogilvie (2010).

We carry out a suite of simulations at a resolution mainly of investigating the effects on hydrodynamic convection when the Rayleigh number Ra is increased, but the Prandtl number Pr is fixed at 2.5.444We have also repeated some of the simulations at a Prandtl number of unity. The differences with the simulations are nominal. Thus viscosity and thermal diffusivity each have the same magnitude in any given simulation, and we decrease both in order to increase the Rayleigh number. The Rayleigh numbers of the simulations are , and . (Note that the simulations undertaken at the two highest Ra are probably underresolved, as explored later.) The Richardson number at initialization is fixed by the initial vertical profile, which is described in detail in Section 4.1.1 (and shown in Figure 1). For the profile parameters chosen, the Richardson number at initialization is . Further details of the simulations are given in Table 4 in Appendix B.

#### 4.4.1 Rayleigh number dependence

As the Rayleigh number is increased from low to high values the system proceeds through the same sequence of states found by Lesur and Ogilvie (2010). We observe no instability for . At instability occurs but the flow appears relatively laminar and axisymmetric; in particular the Reynolds stress is negative througout the linear and nonlinear phases. We conclude that the critical Rayleigh number for the onset of convection lies in the range . At the instability is more vigorous and the flow field significantly more chaotic and non-axisymmetric in the nonlinear phase, at which point the Reynolds stress has become positive. We conclude that the critical Ra at which the sign of flips lies between and . At higher Ra the flow appears even more turbulent and non-axisymmetric. It is a relief that in the nonlinear phase of the instability we find agreement with the inviscid simulations of previous sections at sufficiently high Ra.

Our critical Rayleigh number for the onset of convection is considerably higher than in Lesur & Ogilvie (2010), who report a value of , but this is perhaps explained by their much larger Richardson number which remains constant throughout their box. On the other hand, the critical Ra at which changes sign is much closer to ours, which suggests that the breakdown into non-axisymmetry may be controlled by secondary shear instabilities that are more sensitive to the viscosity than the convective driving.

While the transport of angular momentum is unambigously outward after nonlinear saturation when , it is almost always inwards during the linear phase, as is illustrated by the red curve in Fig. 11. Moreover, during this early stage, the simulation is dominated by the axisymmetric modes of Section 3. This disagrees with our inviscid fiducial simulations (see discussion in Section 4.1.2), which are non-axisymmetric from the outset. We confess that our viscous simulations are more in tune with our physical intuition: (a) in the linear phase growing axisymmetric modes outcompete shearing waves, (b) at sufficiently large amplitudes the modes are subject to secondary shear instabilities in the -plane that buckle the rising and falling plumes, (c) perhaps concurrently or on a short time later (and viscosity permitting), secondary non-axisymmetric instabilities also attack the plumes because they exhibit significant shear in the -plane as well, (d) at this point, the flow degenerates into something more disordered, and importantly, non-axisymmetric and the Reynolds stress flips sign. We conclude that viscosity preferentially damps shearing waves vis-a-vis axisymmetric modes, or effectively kills off the artificial aliasing of shearing waves in inviscid simulations (if this is present).

#### 4.4.2 Convergence with resolution

A final issue is whether our viscous simulations are adequately resolved. More specifically: above what critical Ra is a grid of points inadequate? We conducted simulations at Ra with and grid points and found generally good agreement between the two. The linear growth rates were almost identical and the ultimate nonlinear state statistically similar. The only noticeable difference was in the peak Reynolds stress, which was somewhat larger in the higher resolution run. Overall, we conclude that grid points are adequate to resolve a simulation with Ra.

Things start to deteriorate at a Rayleigh number of . In Figure 11 we plot the time-evolution of the -component of the volume-averaged Reynolds stress for two simulations run at (black) and (red). The lower resolution run possesses no extended period of negative , in contrast to the runs at Ra. We speculate that at this resolution physical viscosity is subdominant to the grid and non-axisymmetric disturbances are artificially enhanced, probably via aliasing. At , the stress is definitely negative in the linear phase and there is a strong negative peak shortly afterward. The physical viscosity is now permitted to work properly and appears to prohibit artificial non-axisymmetric disturbances. As a consequence, the axisymmetric modes preserve their control of the simulation for significantly longer and the simulation is in accord with those at lower Ra. More reassuring is the nonlinear phase a little later in which the two flows closely resemble each other. We conclude that in the linear phase is insufficient to describe a simulation of , but in the later nonlinear phase it is probably adequate.

## 5 Structures in forced compressible convection

In Section 4 we initialized our simulations with a convectively unstable vertical profile, but otherwise did not include any source of heating or cooling. Convection (and to a much lesser extent, conduction) transferred heat and mass vertically so as to zero the buoyancy frequency and send the box into a convectively stable equilibrium. As a result, we were only able to probe non-linear convection for a short period of time. Now we aim to continually sustain convection in order to explore this phase in greater depth. In the absence of self-sustaining convection, this means we have to maintain the convectively unstable profile artificially.

### 5.1 Set-up

SB96 perpetuated convection by forcing the temperature at the mid-plane to adjust to its value at initialization. This strategy, however, raises problems in conservative codes, such as PLUTO and ATHENA, because the energy injection at the midplane has no way to leave the box, except through the distant vertical boundaries (in ZEUS numerical losses on the grid supply a turbulence-dependent ‘cooling’, in addition to the cooling facilitated through the fixed-temperature boundaries). In practice, we found that the disk heats up to an inordinate level and, more importantly, it settles into a marginally unstable state, rather than a driven convective state.

Rather than forcing the code in this way, we mimic the effects of both heating and cooling through thermal relaxation. The idea is to add a source term to the energy equation such that the vertical internal energy profile relaxes to its value at equilibrium on a time-scale . Although this technique is artificial, it serves as a very basic tool for approximating the effects of realistic heating and cooling, as might be supplied by MRI turbulence and radiative transfer. The timescale then would be the characteristic time that the radiative MRI system achieves thermodynamical quasi-equilibrium. In our simulations, however, we take the relaxation timescale to be equal to the linear growth time of convection.

In addition, to mitigate the effects of mass outflows, we incorporate a mass source term to the simulations. At the end of the th time-step we calculated the total vertical mass loss through the upper and lower vertical boundaries. We then added this mass back into the box in the th time-step with the same vertical profile used to initialize the density (c.f. Equation (30)).

We implement the thermal relaxation term by making slight modifications to PLUTO’s built-in optically thin cooling module. The thermal energy equation is updated during a substep to take into account user-defined sources of heating and cooling. The resolution employed in the simulation described in this section is in a box of size , corresponding to about grid cells per scaleheight in the - and -directions and grid cells per scaleheight in the vertical direction. The numerical set-up, boundary conditions and initial conditions are the same as those describe in Section 4.1.1. The thermal relaxation time is taken to be equal to the inverse of the growth rate of the convective instability which we measured to be . Thus the thermal relaxation time is , i.e. the internal energy is relaxed back to its equilibrium profile on about 0.7 orbits. We have also included an explicit viscosity and thermal diffusivity of corresponding to a Rayleigh number of and a Prandtl number of .

Finally, we have repeated the simulations in a cube of resolution and size , as well as at a resolution of in a box of size , and also with periodic boundary conditions in the vertical direction, and observed similar results. We have also run a simulation with the HLL solver to partially explore any code dependence.

### 5.2 Large-scale oscillatory cells

In the left-hand panel of Figure 12 we show the time-evolution of the volume-averaged kinetic energy density associated with the vertical velocity. As in Section 4, exponential growth in the linear phase is followed by non-linear saturation. The forced simulations, however, do not subsequently decay. Instead the vertical kinetic energy increases at a slower rate until about orbit 20, at which point oscillations in the kinetic energy density begin to develop. The cycles increase in frequency and amplitude until about orbit 41 at which point the system settles into a quasi-equilibrium, with the volume-averaged vertical kinetic energy density oscillating about and the period remaining steady at orbits.

Associated with the oscillations are large fluctuations in the -component of the Reynolds stress tensor (shown in the right-hand panel of Figure 12). Instantaneous fluctuations in and in may be either positive or negative, but the time-averaged values over this cyclical phase (orbit 20 to the end of the simulation) are and , respectively. Furthermore, comparing the oscillations in the vertical kinetic energy density to the fluctuations in and also to the volume-averaged gas pressure (not shown here) it is apparent that peaks in the kinetic energy density are correlated with both and troughs in , while troughs in the kinetic energy density are correlated with and peaks in .

A clearer picture of the behaviour of the system emerges when we study the structure of the flow during the cyclical phase. Figure 13 shows snapshots of the -component of the velocity in the -plane. The first panel is taken just after the end of the linear phase (orbit 7.4); the subsequent panels are taken at five successive peaks and troughs in the kinetic energy. As the kinetic energy rises in the non-linear phase, the thin convective cells with radial wavelengths , where is the size of grid a cell in the -direction, slowly begin to merge into larger coherent structures which couple the two halves of the disc together. By orbit 20 (start of the cyclical phase), the radial wavelength of the convective cells has increased to , and by the time the quasi-steady equilibrium state sets in (at around orbit 40) the radial wavelength of the convective cells is .

Comparing the snapshots of to the peaks and troughs in the kinetic energy, it becomes apparent that peaks in the convective energy are associated with the large-scale axisymmetric convective cells, hence . Troughs in the kinetic energy are associated with destruction of those convective cells and with positive stress (outward angular momentum transport). Thus we are observing large-scale convective eddies that appear to be created and destroyed cyclically. Furthermore, it is evident from Figure 13 that after the eddies are destroyed, they are re-formed with the opposite rotation.

The reader may be alarmed that the large eddies extend all the way to the vertical boundaries of the domain, and indeed there is an uncomfortable level of mass loss during this phase. To check that the flows are not an artefact of our box size, we ran additional simulations with a vertical domain of . A density floor had to be activated in such runs, which unfortunately compromised our thermal equilibrium and complicated the onset of convection early in the simulation. However, they ultimately settled onto the cyclical state described above, but now the large-scale eddies peter out before reaching the vertical boundaries. As a result, the mass loss drops to negligible amounts. This confirms that these flows are physical and robust, though only marginally contained within our smaller boxes.

Next, to explore any code dependence of this final outcome, we have rerun the simulation with the more diffusive HLL solver. We observe a negative Reynolds stress during the linear phase and well into the non-linear phase. Associated with this is a remarkably axisymmetric flow field, as confirmed by viewing slices of the pressure perturbation in the -plane at different times. The simulation, nonetheless, enters the cyclical phase during which this axisymmetry is broken. As expected there is a flip in the sign of the Reynolds stress from negative to positive. The behavior thereafter mirrors that of the simulation of forced compressible convection run with the Roe solver: the cyclical creation and destruction of large scale convective cells and an oscillatory Reynolds stress. We conclude that for forced runs, the ultimate quasi-steady state depends negligibly on the algorithm.

### 5.3 Zonal Flows

In Figure 14 we plot, side-by-side, space-time diagrams in the -plane of the -averaged pressure perturbation and of the -averaged perturbation to the -component of the velocity . It is immediately evident from Figure 14 that around the onset of the cyclical phase (orbit 20) alternating streaks in and begin to develop. Respectively, these mark alternating bands (in ) of high () and low () pressure, and of super-Keplerian () and sub-Keplerian () flow.

From Figure 14 it is clear that the pressure and velocity perturbations are 90 degrees out of phase, which is characteristic of zonal flows. We hesitate, however, to claim that these structures are in geostrophic balance (i.e. that ) because, as is evident from the space-time diagrams, the flows are not stationary in time but appear to fluctuate over about one orbit – the same timescale over which the large convective cells are created and destroyed. Indeed the creation and destruction of the zonal flows tracks precisely that of the large-scale convective cells, showing clearly that the two phenomena are connected. The axisymmetry observed during the formation of large-scale convective cells and zonal flows is consistent with the inward transport of angular momentum observed while these structures remain coherent: axisymmetric convective modes dominate during the lifetime of these structures and transport angular momentum inwards.

In Figure 15 (right-hand panels) we plot snapshots in the -plane taken during the cyclical phase of the perturbation of the -component of the velocity and of the perturbed pressure . The axisymmetric structure of the zonal flows is clearly visible in panels (b) and (d), as is the phase difference between the pressure perturbation and the perturbed -component of the velocity .

### 5.4 Vortices

Because the zonal flows consist of alternating bands of sub- and super-Keplerian motion – with strong shear between the bands – we expect that they could give rise to vortices via the Kelvin-Helmholtz instability (modified by rotation and stratification). In Figure 15 we plot a snapshot in the -plane taken during the cyclical phase of the the -component of the vertical component of vorticity . Small but coherent anti-cyclonic blobs of vorticity are observed during the cyclical phase, and these occur precisely where the flow transitions from sub-Keplerian to super-Keplerian (i.e. at the edges of the zonal flows). Typically only anti-cyclonic vortices are observed, which is consistent with the fact that cyclonic vortices tend to be more unstable in Keplerian shear flows.

The vortices tend to be elongated at the start of the cyclical phase (around orbit 20), with an aspect ratio of about 0.4, and grow increasingly circular with time; in fact, once quasi-steady equilibrium has been reached at around orbit 41, many of the vortices have aspect ratios approaching unity. They appear to have limited three-dimensional extent, and we did not observe any vortices that remained coherent for depths exceeding half a scale-height.

### 5.5 Discussion

It is tempting to link our results to the large-scale emergent structures in recent simulations of rotating hydrodynamic convection with uniform (rather than Keplerian) rotation and in the Boussinesq (rather than compressible) regime (Favier et al., 2014; Guervilly et al., 2014). Both Favier et al. (2014) and Guervilly et al. (2014) observe growth in the vertical component of the kinetic energy at high Rayleigh numbers (), as we do. In addition, Favier et al. (2014) and Guervilly et al. (2014) notice the formation of depth-invariant large-scale vortices in the -plane of their simulations (i.e. in the plane perpendicular to the rotation axis).

There are both qualitative and quantitative differences between our results and those of Favier et al. (2014) and Guervilly et al. (2014). Although the Rayleigh number of our simulation with thermal relaxation was , which is consistent with the Rayleigh numbers at which large-scale vortices were observed in the simulations of Favier et al. (2014) and Guervilly et al. (2014), we observe large-scale convective cells (in the -plane) rather than large-scale vortices (in the -plane). We, too, observe anti-cyclonic vortices in the -plane, but these are small and do not appear to have any three-dimensional extent compared to the depth invariant vortices observed in the Boussinesq simulations with uniform rotation. It is possible that they are prevented from merging and thus growing in size because of the turbulence associated with repeated destruction of the large-scale convective cells. The strong shear might also inhibit their growth.

Due to the cyclical nature of the large-scale convective cells, it is tempting to link our results to the intermittent convection reported in the MRI shearing box simulations of Hirose et al. (2014) and Coleman et al. (2018)). In our runs the forcing is due to explicit thermal relaxation, whereas in the runs of Hirose et al. (2014) the forcing is self-consistently provided by MRI heating and radiative cooling. Thus the forcing, and in particular the time-scales associated with the forcing are rather different. However, both thermal relaxation and the MRI limit cycles are similar in that they lead to a cyclical build up of heat and its subsequent purge through vertical convective transport, with the role of opacity in the simulations of Hirose et al. (2014) being simply to modulate this cycle. Thus our results demonstrate that strongly non-linear hydrodynamic turbulent convection has a cyclical nature that might be generic, and that might therefore be robust to the inclusion of more realistic thermodynamics.

## 6 Conclusions

Motivated by recent radiation magnetohydrodynamic shearing box simulations that indicate that an interaction between convection and the magnetorotational instability in dwarf novae can enhance angular momentum transport, we have studied the simpler case of purely hydrodynamic convection, both analytically and through three-dimensional, fully-compressible simulations in PLUTO.

For the linear phase of the instability, we find agreement between the growth rates of axisymmetric modes calculated theoretically and those measured in the simulations to within a percentage error of , thus providing a useful check on our PLUTO code. The linear eigenmodes are worth examining not only to help understand the physical nature of convection in disks, but also because they may appear in some form during the nonlinear phase of the evolution, especially on large-scales, and during intermittent or cyclical convection.

We then explored the nonlinear saturation of the instability, both when convection is continually forced and when it is allowed to reshape the background gradients so that it ultimately dies out. We focussed especially on the old problem of whether hydrodynamic convection in a disc leads to inward () or outward () angular momentum transport. In both forced and unforced convection we found in general in the nonlinear phase. These results were confirmed by a separate run using the code ATHENA, but contradict the classical simulations of SB96 who reported inward transport in both cases using the code ZEUS.

This discrepancy reveals a set of unfortunate numerical difficulties that complicate the simulation of convection in disks. These, in large part, issue from the fact that the inviscid linear modes of convection grow fastest on the shortest possible scales. Thus, no matter the resolution, the nature of the code’s grid dissipation will always impact on the system’s evolution, certainly in the linear phase and possibly afterwards. We argue that a more diffusive numerical set-up, such as supplied by ZEUS at low resolution or Riemann solvers such as HLL, impose an axisymmetry on the flow which leads to generally inward transport of angular momentum. But, on the other hand, we suspect that less diffusive solvers such as Roe and HLLC artificially alias shearing waves in the linear phase of the evolution, leading to spurious non-axisymmetric flow early in a simulation. Though concerning, we believe this is only a problem in the low amplitude linear phase because physical mode-mode interactions will dominate once the perturbations achieve larger amplitudes. Nonetheless, shearing wave aliasing certainly deserves a separate study.

To properly dispense with these numerical issues one must add explicit viscosisty (and thermal diffusivity), as this regularises the linear problem. We find that at a Richardson number of , onset of convection is observed for a critical Rayleigh number . Just above this value convection is largely axisymmetric and . At a larger second critical Ra between and , the sign of switches and the flow becomes more turbulent and nonaxisymmetric. This sequence of states mirrors that simulated by Lesur and Ogilvie (2010). At large (resolved) Rayleigh numbers, viscous simulations are initially controlled by the axisymmetric modes; these are then attacked by secondary shear instabilities in both the and -planes, which break the axisymmetry and order of these structures, leading to a more chaotic state. At lower Ra, viscosity suppresses the non-axisymmetric shear instabilities and axisymmetry is never broken. (At even lower Ra, convection never begins, of course.)

In forced convective runs, rather than maintaining the convection by a fixed heating source at the mid-plane, we instead allowed the vertical equilibrium to relax to its initial, convectively unstable, state. Our thermal relaxation is artificially imposed, but its overall effect is to mimic the heating of the mid-plane and cooling of the corona due to physical mechanisms that maintain the convectively unstable entropy profile, such as the MRI and radiative losses present in the simulations of Hirose et al. (2014). We observed in the non-linear stage now the formation of large-scale convective cells (similar in some respects to elevator flow) that emerge and break down cyclically, in addition to zonal flows and vortices. Although further checks are required, it is tempting to link this cyclical convection to the convective limit cycle observed in the radiative magnetohydrodynamic simulations of Hirose et al. (2014), since both processes rely on the build-up and rapid evacuation of heat.

Despite our demonstration that hydrodynamic convection can lead to positive stress and outward transport of angular momentum, the fact remains that the stresses are small (typically we measured ). Having said that, the magnitude of is sensitive to the depth of the buoyancy frequency profile, and a deeper profile could increase by an order of magnitude or more.

Finally we have not observed self-sustaining hydrodynamic convection in any of our simulations. By self-sustaining convection we mean that (when ) energy extracted from the shear by convection might itself cause convective motions, which in turn extract energy energy from the shear, closing the loop. It is more likely that if convection is to occur in disks it will be as a byproduct of other processes, such as heating by density waves, emitted in the presence of a planet, or by dissipation of magnetorotational turbulence. We intend to investigate both mechanisms and their instigation of convection in future work.

## Acknowledgements

This work was funded by a Science and Technologies Facilities Council (STFC) studentship. The authors acknowledge useful input from John Papaloizou, Doug Lin, Jim Stone, and Steve Balbus. LEH would like to thank Antoine Riols and William Bethune for stimulating conversations and advice on using the PLUTO code.

## References

• Alexiades et al. (1996) Alexiades V., Amiez G., Gremaud P.-A., 1996, Communications in numerical methods in engineering, 12, 31
• Armitage (2011) Armitage P. J., 2011, Annual Review of Astronomy and Astrophysics, 49, 195
• Aurnou et al. (2015) Aurnou J., Calkins M., Cheng J., Julien K., King E., Nieves D., Soderlund K., Stellmach S., 2015, Physics of the Earth and Planetary Interiors, 246, 52
• Bell & Lin (1994) Bell K., Lin D., 1994, ApJ, 427, 987
• Bodo et al. (2013) Bodo G., Cattaneo F., Mignone A., Rossi P., 2013, The Astrophysical Journal Letters, 771, L23
• Boley & Durisen (2006) Boley A. C., Durisen R., 2006, The Astrophysical Journal, 641, 534
• Boyd (2001) Boyd J. P., 2001, Chebyshev and Fourier spectral methods. Courier Corporation
• Cabot (1996) Cabot W., 1996, The Astrophysical Journal, 465, 874
• Cannizzo (2011) Cannizzo J. K., 2011, The Limit Cycle Instability In Dwarf Nova Accretion Disks. World Scientific, pp 6–40, doi:10.1142/9789814350976_0002, http://www.worldscientific.com/doi/abs/10.1142/9789814350976_0002
• Chiang & Goldreich (1997) Chiang E., Goldreich P., 1997, The Astrophysical Journal, 490, 368
• Coleman et al. (2018) Coleman M. S., Blaes O., Hirose S., Hauschildt P. H., 2018, The Astrophysical Journal, 857, 52
• D’Alessio et al. (1998) D’Alessio P., Cantö J., Calvet N., Lizano S., 1998, The Astrophysical Journal, 500, 411
• Favier et al. (2014) Favier B., Silvers L. J., Proctor M. R. E., 2014, Physics of Fluids, 26, 096605
• Goldreich & Lynden-Bell (1965) Goldreich P., Lynden-Bell D., 1965, Monthly Notices of the Royal Astronomical Society, 130, 125
• Guervilly et al. (2014) Guervilly C., Hughes D. W., Jones C. A., 2014, Journal of Fluid Mechanics, 758, 407
• Hawley et al. (1995) Hawley J. F., Gammie C. F., Balbus S. A., 1995, The Astrophysical Journal, 440, 742
• Hirose (2015) Hirose S., 2015, Monthly Notices of the Royal Astronomical Society, 448, 3105
• Hirose et al. (2014) Hirose S., Blaes O., Krolik J. H., Coleman M. S., Sano T., 2014, The Astrophysical Journal, 787, 1
• Klahr et al. (1999) Klahr H., Henning T., Kley W., 1999, The Astrophysical Journal, 514, 325
• Kley et al. (1993) Kley W., Papaloizou J., Lin D., 1993, The Astrophysical Journal, 416, 679
• Latter & Papaloizou (2017) Latter H. N., Papaloizou J., 2017, Monthly Notices of the Royal Astronomical Society
• Lesur & Ogilvie (2010) Lesur G., Ogilvie G. I., 2010, Monthly Notices of the Royal Astronomical Society: Letters, 404, L64
• Lin & Papaloizou (1980) Lin D., Papaloizou J., 1980, Monthly Notices of the Royal Astronomical Society, 191, 37
• Lin et al. (1993) Lin D., Papaloizou J., Kley W., 1993, The Astrophysical Journal, 416, 689
• Lyra et al. (2016) Lyra W., Richert A. J., Boley A., Turner N., Mac Low M.-M., Okuzumi S., Flock M., 2016, The Astrophysical Journal, 817, 102
• Mignone et al. (2007) Mignone A., Bodo G., Massaglia S., Matsakos T., Tesileanu O., Zanni C., Ferrari A., 2007, The Astrophysical Journal Supplement Series, 170, 228
• Mignone et al. (2012) Mignone A., Flock M., Stute M., Kolb S., Muscianisi G., 2012, Astronomy & Astrophysics, 545, A152
• Ruden et al. (1988) Ruden S. P., Papaloizou J., Lin D., 1988, The Astrophysical Journal, 329, 739
• Ryu & Goodman (1992) Ryu D., Goodman J., 1992, The Astrophysical Journal, 388, 438
• Spruit et al. (1990) Spruit H. C., Nordlund A., Title A., 1990, Annual review of astronomy and astrophysics, 28, 263
• Stone & Balbus (1996) Stone J. M., Balbus S. A., 1996, The Astrophysical Journal, 464, 364
• Stone & Gardiner (2010) Stone J. M., Gardiner T. A., 2010, The Astrophysical Journal Supplement Series, 189, 142
• Stone et al. (2008) Stone J. M., Gardiner T. A., Teuben P., Hawley J. F., Simon J. B., 2008, The Astrophysical Journal Supplement Series, 178, 137
• Zingale et al. (2002) Zingale M., et al., 2002, The Astrophysical Journal Supplement Series, 143, 539

## Appendix A Convectively unstable vertical disc profiles

We describe the two convectively unstable vertical profiles that are used to initialize our simulations. The reader should note that these profiles may or may not correspond to those in real astrophysical discs, which will be determined by several sources of heating and cooling, none considered in this paper. Here we simply present convectively unstable disc profiles that satisfy hydrostatic equilibrium, are convectively unstable, and can conveniently initialize simulations.

### a.1 Stone and Balbus (1996) profile

The SB96 profile employs a power law profile in temperature, see  Eq. (9). For , the density can be obtained analytically

 ρ=ρ0(1−s3)−(1+g/3)(1−s)gexp{2g[s−1√3tan−1(√3ss+2))]}

where is the midplane density, , and , in which , and is the midplane scaleheight. The pressure is obtained from the ideal gas equation of state. Note that for we must have vacuum, which ties the maximum numerical domain to the ratio .

The buoyancy frequency for is given by

 N2BΩ2=11−(zH)3/2f[(1−1γ)(zH)2−32(zH)3/2f] (28)

which corresponds to a profile that is negative (convectively unstable) within some region about the mid plane and positive (convectively unstable) outside of this region. The width of the convectively unstable region is given by

 |Lc|=[32(1−1γ)−1H2AT0]2. (29)

Although convenient within a limited choice of parameters, the SB96 profile suffers from the drawback that the width of the convectively unstable region is sensitive to the size of the box through the ratio . Increasing the vertical box size necessarily decreases the size of the convectively unstable region.

### a.2 Gaussian temperature profile

The drawbacks of the SB96 profile motivated us to search for a more convenient unstable profile, one that would leave the size and depth of the convectively unstable region independent of the vertical extent of the box. Setting the temperature to a Gaussian, cf. Eq. (8) provided such a profile.

The associated density is

 ρ=ρ0ez2/βH2exp(−βH2Ω22T0ez2/βH2), (30)

where is midplane density and the midplane scale height. Pressure is obtained from the ideal gas equation of state, as above.

The buoyancy frequency is given by

 N2BΩ2=2βH2z2[((1−1γ)βH2Ω2μ2T0Rez2/βH2−1]. (31)

The boundary of the convectively unstable region is hence given by

 |Lc|≡ ⎷βH2ln[(1−1γ)−12T0RβμH2Ω2]. (32)