MNRAS 480, 4797–4816 (2018) doi:10.1093/mnras/sty2097
Advance Access publication 2018 August 2
Hydrodynamic convection in accretion discs
Loren E. Held‹ and Henrik N. Latter
Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Centre for Mathematical Sciences, Wilberforce Road,
Cambridge CB3 0WA, UK
Accepted 2018 July 30. Received 2018 July 30; in original form 2018 March 27
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, and 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.
Key words: accretion, accretion discs – convection – hydrodynamics – instabilities, turbu-
lence.
1 IN T RO D U C T I O N
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, Nordlund & Title
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 mid-plane, 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 nu-
clear 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 associ-
 E-mail: leh50@cam.ac.uk
ated 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 out-
burst phase of dwarf novae (for a succinct review, see Cannizzo
1993). More recently, simulations of magnetorotational turbulence
in stratified discs have shown that the MRI is capable of gener-
ating a convectively unstable entropy gradient. Moreover, three-
dimensional, shearing box simulations with zero-net vertical mag-
netic 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 proto-
planetary 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 magnetoro-
tational instability is unlikely to operate (Armitage 2011). Although
irradiation by the central star usually results in an isothermal ver-
tical 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
C© 2018 The Author(s)
Published by Oxford University Press on behalf of the Royal Astronomical Society
D
ow
nloaded from
 https://academ
ic.oup.com
/m
nras/article-abstract/480/4/4797/5064252 by U
niversity of C
am
bridge user on 14 D
ecem
ber 2018
4798 L. E. Held and H. N. Latter
unstable gradient (Boley & Durisen 2006; Lyra et al. 2016). So too
might the heating from the dissipation of the strong zonal mag-
netic fields associated with the Hall effect (Lesur, Kunz & Fromang
2014).
As mentioned, accretion may instigate convection. But convec-
tion 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
disc contracting towards the mid-plane. They showed that vertical
convection arises if the opacity increases sufficiently fast with tem-
perature, with the heating source initially provided by the gravita-
tional contraction. But if convection is modelled 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 ini-
tial gravitational contraction as a source of heat for maintaining
convection.
In order to go beyond mixing length theory, Ruden, Papaloizou &
Lin (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 α ∼ 10−3–
10−2. 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 cen-
tred not so much on the size of α but rather on its sign. Dissipation
of non-linear axisymmetric convective cells was investigated by
Kley, Papaloizou & Lin (1993), who performed quasi-global simu-
lations (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 viscos-
ity. Convective shearing waves were first investigated by Ryu &
Goodman (1992), who concluded that linear non-axisymmetric per-
turbations would result in a net inward angular momentum flux at
sufficiently large time. These results were questioned however by
Lin, Papaloizou & Kley (1993), who examined analytically and nu-
merically a set of localized linear non-axisymmetric disturbances in
global geometry. They demonstrated that these modes could trans-
port angular momentum outwards in some cases, and opined that the
inward transport of angular momentum reported by Ryu & Good-
man (1992) was an artefact of the shearing box approximation.
Interest in hydro convection waned after local non-linear 3D com-
pressible 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 & Bal-
bus (1996, hereafter SB96) initialized inviscid, fully compressible,
and vertically stratified shearing box simulations with a convec-
tively unstable vertical temperature profile, and measured a time-
averaged value of α ∼ −4.2 × 10−5. 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 in-
ward 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, Henning
& Kley (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 15 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 indi-
cated 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 suffi-
ciently large (Lesur & Ogilvie 2010). It would appear then that, in
certain circumstances, vertical convection can drive outward angu-
lar momentum transport after all.
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 momen-
tum transport. We also isolate and characterize other generic fea-
tures of convection that should be shared by multiple disc classes
(dwarf novae, protoplanetary discs, 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 approxi-
mation. We restrict ourselves to an idealized 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 behaviour 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 sus-
tained 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 ME T H O D S
2.1 Governing equations
We work in the shearing box approximation (Goldreich & Lynden-
Bell 1965; Hawley, Gammie & Balbus 1995; Latter & Papaloizou
2017), which treats a local region of a disc as a Cartesian box
located at some fiducial radius r = r0 and orbiting with the angular
frequency of the disc at that radius 0 ≡ (r0). A point in the box
has Cartesian coordinates (x, y, z) that are related to the cylindrical
coordinates (r, φ, z) through x = r − r0, y = r0(φ − φ0 − 0t), and
z = z. The equations of gas dynamics are now
∂tρ + ∇ · (ρu) = 0, (1)
∂tu + u · ∇u = − 1
ρ
∇P − 20ez × 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 P = (R/μ)ρT or, equivalently, P = e(γ − 1)ρ in
terms of the specific internal energy e. Here, P is the thermal pressure
MNRAS 480, 4797–4816 (2018)
D
ow
nloaded from
 https://academ
ic.oup.com
/m
nras/article-abstract/480/4/4797/5064252 by U
niversity of C
am
bridge user on 14 D
ecem
ber 2018
Hydrodynamic convection in accretion discs 4799
of the fluid, T is the temperature, ρ is mass density, R 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
γ = 5/3 in all our simulations.1 The effective gravitational potential
is embodied in the tidal acceleration geff = 2q20xxˆ − 20zzˆ, where
q is the shear rate. For Keplerian discs q = 3/2 a value we adopt
throughout the paper.
Other terms in the equations include the viscous stress tensor
T ≡ 2ρνS. Here, ν is the kinematic viscosity and S ≡ (1/2)[∇u +
(∇u)T] − (1/3)(∇ · u)I is the traceless shear tensor. The thermal
conductivity is ξ , but in our simulations we specify thermal diffu-
sivity χ rather than ξ . The former is related to the latter via χ =
ξ /(cpρ), where cp is the specific heat capacity at constant pressure.
In some of our simulations, we include heating and cooling
through a thermal relaxation term 
relax ≡ (ρe − ρ0e0)/τ relax.
This relaxes the thermal energy in each cell back to its initial state
e0(z) ≡ e(z, t = 0) on a time-scale given by τ relax. We refer to sim-
ulations run with 
relax = 0 as simulations of forced compressible
convection. Conversely, simulations with 
relax = 0 are referred to
as unforced. The thermal relaxation time-scale τ relax 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
γ
∂ ln P
∂z
− ∂ ln ρ
∂z
]
20. (4)
But convection is opposed by the effects of viscosity and thermal
diffusivity, with the ratio of destabilizing and stabilizing processes
quantified by the Rayleigh number
Ra ≡ N
2
BH
4
νχ
, (5)
where H is the disc’s scale height (formally defined later). On the
other hand, the ratio of buoyancy to shear is expressed through the
Richardson number
Ri ≡ N
2
B
q220
. (6)
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 ex-
pression for the scaled intensity of N2B . Finally, the Prandtl number
is defined as the ratio of kinematic viscosity to thermal diffusivity
Pr ≡ ν
χ
. (7)
When the explicit diffusion coefficients are neglected (ν, χ =
0), convective instability occurs when the buoyancy frequency is
negative (i.e. in regions where N2B < 0). When explicit viscosity
and thermal diffusivity are included, however, convective instabil-
ity requires both that N2B < 0 and that the Rayleigh number exceeds
some critical value Rac. Though in real astrophysical discs the vis-
cosity is negligible, the inclusion of magnetic fields complicates the
stability criterion, as does the presence of preexisting turbulence (as
1Technically γ = 5/3 is valid for ionized hydrogen. Since all our simulations
are hydrodynamic it might be more appropriate to use γ = 7/5, which is
valid for molecular hydrogen, but the differences are small and do not affect
our results.
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 convec-
tion resists the disordered background flow. In any case, the sign of
N2B 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 em-
ploy the Roe Riemann solver, 3rd order in space WENO interpola-
tion, and the 3rd order in time Runge–Kutta algorithm. Other con-
figurations are explored in Section 4. To allow for longer time-steps,
we employ the FARGO scheme (Mignone et al. 2012). When ex-
plicit viscosity ν and thermal diffusivity χ are included, we further
reduce the computational time by using the Super-Time-Stepping
(STS) scheme (Alexiades, Amiez & Gremaud 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, evolv-
ing the total energy equation rather than the thermal energy equation,
where the total energy density of the fluid is given by E = (1/2)ρu2
+ ρe (kinetic + internal). In PLUTO the thermal relaxation term 
relax
= (ρe − ρ0e0)/τ relax 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. The first exhibits a Gaussian temperature profile
T = T0exp
[
− z
2
βH 2
]
, (8)
where T0 is mid-plane temperature and β is a dimensionless tuning
parameter. See Fig. 1 for all associated thermal profiles. For com-
parison, we also use the profile introduced by SB96 in which the
temperature follows a power law
T = T0 − Azp, (9)
where A and p are parameters, with p = 3/2 usually (see Fig. 2).
Further details of both profiles are given in Appendix A. The equi-
libria are convectively unstable within a confined region (of size
Lc) 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 ampli-
tudes typically less than 10−3 at the vertical boundaries and less
than 10−7 at the mid-plane.
The background velocity is given by u = −(3/2)0x ey . At ini-
tialization, we usually perturb all the velocity components with
random noise exhibiting a flat power spectrum. The perturbations
MNRAS 480, 4797–4816 (2018)
D
ow
nloaded from
 https://academ
ic.oup.com
/m
nras/article-abstract/480/4/4797/5064252 by U
niversity of C
am
bridge user on 14 D
ecem
ber 2018
4800 L. E. Held and H. N. Latter
Figure 1. Vertical disc structure for a Gaussian temperature profile in a disc of height 2H about the mid-plane. Left: vertical profiles for density ρ, pressure P,
and temperature T. Right: vertical profile of the buoyancy frequency N2B/2. For clarity, only half of the vertical domain is shown. The profile parameters are
{T0 = 1.0, ρ0 = 1.0, β = 3.0} and the adiabatic index is γ = 5/3.
Figure 2. Vertical disc structure of Stone & Balbus (1996) in a disc of height 2H about the mid-plane. Left: vertical profiles for density ρ, pressure P, and
temperature T. Right: vertical profile of the buoyancy frequency N2B/2. For clarity, only half of the vertical domain is shown. The profile parameters are {T0
= 1.0, ρ0 = 1.0, g = 5.0} and the adiabatic index is γ = 5/3.
δu have maximum relative amplitude of about 5 × 10−5 and can
be either positive or negative. In order to investigate specifically the
nature of linear axisymmetric convective modes, we initialize sev-
eral PLUTO simulations with linear axisymmetric modes calculated
semi-analytically rather than with random noise (see Section 3.5).
Time units are selected so that 0 = 1. From now, the subscript
on the angular frequency is dropped if it appears. The length unit
is chosen so that the initial mid-plane sound speed cs = 1, which in
turn defines a constant reference scale height H ≡ cs/0 = 1. 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 scale height H, defined above. We
employ resolutions of 643, 1283, 2563, and 5123 in boxes of size
4H × 4H × 4H, which correspond to resolutions of 16, 32, 64,
and 128 grid cells per scale height H in all directions. However, in
simulations in which we force convection, we employ a ‘wide-box’
(6H × 6H × 4H) with a resolution of 256 × 256 × 256 which
corresponds to about 43 grid cells per H in the horizontal direction
and 64 grid cells per H in the vertical direction.
2.3.4 Boundary conditions
We use shear-periodic boundary conditions in the x-direction (see
Hawley et al. 1995) and periodic boundary conditions in the y-
direction. In the vertical direction, we keep the ghost zones asso-
ciated with the thermal variables in isothermal hydrostatic equilib-
rium (the temperature of the ghost zones being kept equal to the
temperature of the vertical boundaries at initialization) in the man-
ner 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.
MNRAS 480, 4797–4816 (2018)
D
ow
nloaded from
 https://academ
ic.oup.com
/m
nras/article-abstract/480/4/4797/5064252 by U
niversity of C
am
bridge user on 14 D
ecem
ber 2018
Hydrodynamic convection in accretion discs 4801
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 X,
the volume-average of that quantity is denoted 〈X〉 and is defined
as
〈X〉(t) ≡ 1
V
∫
V
X(x, y, z, t)dV , (10)
where V is the volume of the box.
If we are interested only in the vertical structure of a quantity X
then we average over the x- and y-directions, only. The horizontal
average of that quantity is denoted 〈X〉xy and is defined as
〈X〉xy(z, t) ≡ 1
A
∫
A
X(x, y, z, t)dA, (11)
where A is the horizontal area of the box. Horizontal averages over
different coordinate directions (e.g. over the y- and z-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 X is
denoted 〈X〉t and is defined as
〈X〉t (x, y, z) ≡ 1
t
∫ tf
ti
X(x, y, z, t)dt, (12)
where we integrate from some initial time ti to some final time tf
and t ≡ tf − ti.
2.4.2 Reynolds stress, alpha, and energy densities
In accretion discs, the radial transport of angular momentum is
related to the xy-component of the Reynolds stress, defined as
Rxy ≡ ρuxδuy, (13)
where δuy ≡ uy + qx is the fluctuating part of the y-component of
the total velocity uy. Rxy > 0 corresponds to positive (i.e. radially
outward) angular momentum transport, while Rxy < 0 corresponds
to negative (i.e. radially inward) angular momentum transport.
The Reynolds stress is related to the classic dimensionless pa-
rameter α, 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
α ≡ 〈Rxy〉
q〈P 〉 . (14)
The kinetic energy density of a fluid element is defined as
Ekin ≡ 12ρu
2, (15)
where u is the magnitude of the total velocity of a fluid element.
Often we will plot the vertical kinetic energy density, in which case
u in the above is replaced by uz.
The total energy (kinetic plus thermal) is not conserved in the
shearing box even when there are no (explicit or numerical) dis-
sipative effects: work done by the effective gravitational potential
energy on the fluid and flux of angular momentum through the x-
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 pres-
sure, 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
γ
d ln 〈P 〉xy
dz
− d ln 〈ρ〉xy
dz
]
. (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 L I N E A R TH E O RY
In this section, we investigate the linear phase of the axisymmet-
ric 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 eigen-
functions (density, velocity, and pressure perturbations) of the ax-
isymmetric convective instability in the shearing box as a function
of radial wavenumber kx. The convectively unstable background
vertical profile used in all calculations is shown in Fig. 1. The pro-
file is discussed in detail in Section A2. All calculations were carried
out with profile parameters T0 = 1.0, ρ0 = 1.0, β = 3.0 and an adi-
abatic index of γ = 5/3. For our chosen background equilibrium,
this corresponds to a Richardson number of Ri ∼ 0.05.
We proceed as follows: first, we produce linearized equations for
the perturbed variables; secondly, these are Fourier analysed so that
X′(x, z, t) = eσ t eikxxX′(z) for each of the perturbed fluid variables
ρ ′, u′, and p′ . 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 = 2u′y −
ikx
ρ0
P ′, (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 X0 = X0(z). For equilibrium
ρ0 and P0 see Appendix A.
3.1 Spectral methods
We discretize the fluid variables on a vertical grid containing Nz
points by expanding each eigenfunction as a linear superposition
of Whittaker cardinal (sinc) functions. These are appropriate as
all perturbations should decay to zero far from the convectively
unstable region. This choice also saves us from explicitly imposing
boundary conditions. The discretized equations are next gathered
up into a single algebraic eigenvalue problem (Boyd 2001). Finally,
we solve this matrix equation numerically using the QR algorithm
MNRAS 480, 4797–4816 (2018)
D
ow
nloaded from
 https://academ
ic.oup.com
/m
nras/article-abstract/480/4/4797/5064252 by U
niversity of C
am
bridge user on 14 D
ecem
ber 2018
4802 L. E. Held and H. N. Latter
Figure 3. Growth rates as a function of radial wavenumber kx (scaled by the width of the convectively unstable region Lc). The solid lines were calculated
analytically using WKBJ methods and the dashed lines were calculated semi-analytically using pseudo-spectral methods. The squares, triangles, and diamonds
correspond to measurements taken from PLUTO simulations run at a vertical resolution of Nz = 128. The blue circle was taken from a PLUTO simulation run
at Nz = 256.
to obtain the growth rates σ and eigenfunctions for a given radial
wavenumber kx. The code that 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
d2U
dz2
+ k2z (z)U = 0, (24)
where
U =
(
γP0
σ 2 + κ2 + k2xc20
)1/2
u′z, (25)
in which c0 is the equilibrium sound speed and the ‘vertical
wavenumber’ of the perturbations is given by
k2z (z) ≈ −k2x
N2B (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 kx large. Equation (24) has two
turning points at k2z (z) = 0, which occur when −N2B (z) = σ 2. Phys-
ically, one turning point occurs at the disc mid-plane (z = 0), while
the other turning point occurs at the boundary of the convectively
unstable region z = Lc. When −N2B (z) > σ 2 then k2z (z) > 0, 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 decay exponentially in space outside.
The standard WKBJ solution is given by
U (z) ∼ c1 exp
(
+i
∫ z
kz(z′)dz′
)
+ c2 exp
(
−i
∫ z
kz(z′)dz′
)
,
where c1 and c2 are constants. By matching the interior (oscillatory)
solution on to the exterior (exponentially decaying) solution, we
obtain the eigenvalue equation (dispersion relation):∫ z1
z0
kz(z)dz =
(
n + 1
2
)
π n = 0, 1, 2, 3, · · · , (27)
where lower and upper bounds of integration z0 and z1, respectively,
are related by the implicit equation N2B (z0) = N2B (z1) = −σ 2, and n
is the mode number. Substituting for k2z (z) using equation (26), we
obtain an implicit equation relating the radial wavenumber kx and
the growth rate σ which we solve numerically via a root-finding
algorithm.
3.3 Eigenvalues
For a given radial wavenumber kx, 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 kx, these corresponding to vertical quantum numbers n = 0, 1,
and 2. Note that in the figure kx has been normalized by the width
of the convectively unstable region Lc (given by equation A5 in
Appendix A2).
As the radial wavenumber kx 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 kx (small radial wavenumber λx), 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 be-
MNRAS 480, 4797–4816 (2018)
D
ow
nloaded from
 https://academ
ic.oup.com
/m
nras/article-abstract/480/4/4797/5064252 by U
niversity of C
am
bridge user on 14 D
ecem
ber 2018
Hydrodynamic convection in accretion discs 4803
Figure 4. Vertical profiles of the perturbations Re{ρ ′ (z)}, Re{u′z(z)},
Re{P′ (z)}, and of Im{u′x (z)} and Im{u′y (z)} in the upper half plane z ∈
[0, 2]H, for the n = 0 odd mode at kx = 20.0. The imaginary and real parts
not shown are effectively zero. The black dot marks the upper extent of
the convectively unstable region Lc, while the red dot marks the most con-
vectively unstable points (see equation A5). The eigenfunctions have been
rescaled so that the maximum of each eigenfunction is unity.
tween even or odd modes, nonetheless we find very close agree-
ment between the WKBJ and semi-analytical results for all radial
wavenumbers, even at low kx 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 kx, min
below which both the numerical and WKBJ growth rates are zero.
This feature is just visible in the bottom left-hand side of Fig. 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
non-linear phase of convection the short scale linear axisymmetric
modes may be subordinate and this less 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 (n = 0) odd mode with a radial
wavenumber kx = 20.
In Fig. 4, we plot the vertical profiles of the non-zero compo-
nents 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 Fig. 5, we plot their real parts in the xz-plane with the
correct amplitudes. We find that |u′z| > |u′y | ∼ |u′x |  |ρ ′|  |P ′|.
Thus, vertical velocity perturbations are greatest in magnitude,
while pressure perturbations are smaller than vertical velocity per-
turbations 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 Fig. 4, u′z (solid black curve) is antisymmetric about the mid-
plane, whereas the remaining perturbations are symmetric. Figs 4
and 5 clearly show that both odd and even modes consist of a ra-
dial chain of convective cells above and below the mid-plane and
localized near the most convectively unstable point (denoted by the
large red dot in the former figure). The peak in the amplitude of
u′z that occurs near this point tells us that the fluid elements reach
Figure 5. Two-dimensional profiles in the xz-plane of the real parts of the selected components of odd and even eigenfunctions for n = 0 and kx = 20.0.
MNRAS 480, 4797–4816 (2018)
D
ow
nloaded from
 https://academ
ic.oup.com
/m
nras/article-abstract/480/4/4797/5064252 by U
niversity of C
am
bridge user on 14 D
ecem
ber 2018
4804 L. E. Held and H. N. Latter
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. P′ > 0), and
where hot elements begin to sink (lower background pressure), the
pressure perturbations decrease (i.e. P′ < 0). The vertical velocity
perturbation (black solid line) is out of phase with the radial and az-
imuthal 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 non-zero, but relatively small, u′z at the
mid-plane, and thus weakly couple the two sides of the disc. Such
modes may permit some exchange of mass, momentum, and thermal
energy across the mid-plane when reaching large amplitudes (in
particular, see Section 5.2). In particular, falling fluid elements may
’overshoot’ the midplane and thus circulate from one side of the
disk to the other.
As the radial wavenumber is increased from kx = 20.0 to 200.0
(very small radial wavelengths), the perturbations become increas-
ingly localized about the most convectively unstable point(s). The
localization in z however is not as narrow as the mode’s radial
wavelength. For larger kx, in fact, each pair of even and odd modes
changes their character and becomes entirely localized to either
the upper or lower disc. Activity in the two halves of the disc is
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 λx × 0.625H ×
4H. We maintain a fixed number of 32 grid cells per scale height
in the z-direction, 16 grid cells per scale height in the y-direction,
and 64 grid cells per radial wavelength in the x-direction. Thus, the
radial wavelength of the modes is always well resolved in our sim-
ulations. No explicit diffusion coefficients (i.e. viscosity or thermal
diffusivity) are used.
The results are plotted in Fig. 3 as squares (and one circle).
Numerical values appear in Table 1. We find excellent agreement
between simulation and theory: the relative error is < 1 per cent for
the fundamental mode in the interval kx ∈ [20.0, 60.0] and for the
first harmonic at kx = 25.0. As the radial wavenumber is increased
the simulation growth rates begin to deviate from the theoretical
curves, although the percentage error remains less than 6 per cent
even at large wavenumbers. This behaviour 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 λx is always well resolved in
the simulations, the vertical wavelength becomes increasingly less
so.
The numerical diffusion could be reduced by increasing the ver-
tical resolution of the simulations. To check this we have rerun the
PLUTO simulation, initialized with eigenfunctions corresponding to
the first harmonic (n = 1) even mode with kx = 100, at twice the ver-
tical resolution, i.e. Nz = 256 instead of Nz = 128. At a resolution of
Nz = 128, we measured a growth rate of σPLUTO = 0.2840 , cor-
responding to a percentage error of 5.8 per cent with the growth
rate calculated semi-analytically (σEIG = 0.3016 ). At a verti-
cal resolution of Nz = 256, however, we measured a growth
Table 1. Comparison of eigensolver growth rates with growth rates mea-
sured in PLUTO simulations. The simulations were initialized with the eigen-
functions corresponding to each mode n and radial wavenumber kx. All
simulations except the one marked with a footnote were run at a resolution
of 64 × 10 × 128. The standard deviation on the growth rates σ PLUTO mea-
sured in the simulations is typically in the range from ±0.0002 to ±0.0004.
Mode Parity kx σEIG σPLUTO
Percentage
error
n = 0 even 20.0 0.2727 0.2719 0.3
odd 20.0 0.2724 0.2715 0.3
even 40.0 0.3076 0.3059 0.6
even 60.0 0.3186 0.3160 0.8
even 100.0 0.3273 0.3223 1.6
n = 1 even 25.0 0.1694 0.1678 0.9
odd 25.0 0.1643 0.1624 0.9
even 40.0 0.2392 0.2343 1.1
even 60.0 0.2747 0.2655 3.3
even 100.0 0.3016 0.2840 5.8
– – – 0.2974 1.4a
n = 2 even 40.0 0.1593 0.1525 4.2
odd 40.0 0.1580 0.1517 3.9
even 60.0 0.2272 0.2251 0.9
even 100.0 0.2749 0.2637 4
Note: a Simulation run with a resolution of Nz = 256 grid points in the
vertical direction.
rate of σPLUTO = 0.2974 , corresponding to a percentage error
of just 1.4 per cent. Thus, the growth rates measured in the sim-
ulations 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 ob-
served inward angular momentum transport. This is in agreement
with the analytical argument presented in Stone & Balbus (1996)
regarding axisymmetric flow.
4 SI M U L AT I O N S O F U N F O R C E D
COMPRESSI BLE CONVECTI ON
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 ini-
tialized 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 de-
pends 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
MNRAS 480, 4797–4816 (2018)
D
ow
nloaded from
 https://academ
ic.oup.com
/m
nras/article-abstract/480/4/4797/5064252 by U
niversity of C
am
bridge user on 14 D
ecem
ber 2018
Hydrodynamic convection in accretion discs 4805
Figure 6. Left: semilog plot of the time evolution of the volume-averaged vertical kinetic energy density in PLUTO simulations at resolutions of 643 (blue
line), 1283 (red line), and 2563 (black line). Right: time evolution of volume-averaged xy-component of the Reynolds stress tensor.
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
643, 1283, and 2563 in boxes of size 4H × 4H × 4H, where H is the
scale height.2 All the simulations were initialized with the convec-
tively unstable vertical profiles for density and pressure described in
Appendix A2 with profile parameters T0 = 1.0, ρ0 = 1.0, β = 3.0
and an adiabatic index of γ = 5/3 (see Fig. 1). The Stone & Bal-
bus vertical profile SB96 was also trialled yielding very similar
results, most of which have been omitted in the interests of space.
Vertical outflow conditions were employed for our fiducial simu-
lations but periodic and free-slip boundary conditions in the verti-
cal direction produced the same behaviour. Random perturbations
to the velocity components were seeded at initialization, with a
maximum amplitude |δu| ∼ 10−5. Finally, we adopt a very small
but finite thermal diffusivity of χ = 2 × 10−6 to facilitate con-
duction of thermal energy through the vertical boundaries and to
aid code stability. No explicit non-adiabatic heating, cooling, or
thermal relaxation is 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-hand panel of Fig. 6, we track the time evolution of the
volume-averaged kinetic energy density associated with the vertical
velocity for simulations run at resolutions of 643, 1283, and 2563,
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 con-
vective instability begins in earnest, characterized by exponential
2A simulation run at a resolution 5123 is included in the list of fiducial
simulations in Table B1 (see Appendix B) but is not described in this section
for brevity.
growth of the perturbation amplitudes. During this phase, inter-
nal 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 sat-
uration, the kinetic energy decreases gradually: the convective cells
redistribute thermal energy and mass, thus shifting the thermal pro-
file 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 en-
ergy successively smaller over the same period of time in the 643,
1283, and 2563 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 xy-component of the Reynolds
stress is plotted in the right-hand 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 axisymmet-
ric 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 sug-
gests 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’
kx = 0 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 effec-
tively axisymmetric. (Even then the stress is extremely oscillatory.)
MNRAS 480, 4797–4816 (2018)
D
ow
nloaded from
 https://academ
ic.oup.com
/m
nras/article-abstract/480/4/4797/5064252 by U
niversity of C
am
bridge user on 14 D
ecem
ber 2018
4806 L. E. Held and H. N. Latter
Figure 7. From left to right: snapshots of uz in the xz-plane taken at non-linear saturation (just after the linear phase) at resolutions of 643, 1283, and 2563,
respectively. Blue denotes uz < 0 (sinking fluid) and red denotes uz > 0 (rising fluid).
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 non-linear regime is achieved before
the dominant shearing waves become strongly trailing and begin to
send angular momentum inward. Given that the typical time-scale
for shearing out is roughly a few orbits at best, this is marginal
but not impossible. An alternative, more plausible, and more trou-
bling hypothesis is that our inviscid numerical code misrepresents
trailing shearing waves: more specifically, the code is artificially
reseeding fresh leading waves from strongly trailing waves (‘alias-
ing’; 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 concern-
ing. But we hasten to add that this phenomenon should only be
problematic in the low-amplitude linear phase; once the perturba-
tions 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 Rxy from orbit 5 to orbit 15, a
period spanning non-linear saturation, show that the Reynolds stress
decreases as the resolution increases (see Table B1 in Appendix B).
The time- and volume-averaged values of Rxy over a period spanning
non-linear saturation (orbits 5–15) are +4.8 × 10−6, +2.7 × 10−6,
and +1.7 × 10−6 in the 643, 1283, and 2563 simulations, respectively.
In terms of the turbulent alpha-viscosity, this corresponds (from
equation 14) to α ∼ +2.3 × 10−5, +1.5 × 10−5, and 9.6 × 10−6,
respectively. The dependence on the numerical viscosity can be
explained by appealing to secondary shear instabilities (see the 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 con-
fined to the region |z| < Lc ∼ 1.11H, 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
−1.4 , −3.4 , and −8.8 per cent in the 643, 1283 and 2563, simula-
tions, respectively.
4.1.3 Structure of the flow
The development of convective instability and the associated con-
vective cells are best observed through snapshots of the vertical
component of the velocity in the xz-plane, shown at resolutions of
643, 1283, and 2563 just after non-linear saturation in Fig. 7. A full
set of convective cells is clearly visible at all resolutions. They are
thin, filamentary structures several grid cells wide comprising alter-
nating negative and positive velocities (updrafts and downdrafts).
The maximum vertical Mach number of the flow around non-
linear saturation in the 2563 run is about Mz ∼ 0.15, 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 de-
velop 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 ampli-
tudes 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 the right-hand 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
Fig. 8(a) shows the vertical profiles of horizontally averaged heat
and mass fluxes taken from snapshots just after non-linear satura-
tion from the simulation with resolution 2563. 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 sat-
uration (see Fig. 6).
Negative (positive) values for the fluxes for z < 0 and positive
(negative) values for the fluxes for z > 0 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 Fig. 8). In addition, there is some mass
flux towards the mid-plane within −H < z < H. 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
MNRAS 480, 4797–4816 (2018)
D
ow
nloaded from
 https://academ
ic.oup.com
/m
nras/article-abstract/480/4/4797/5064252 by U
niversity of C
am
bridge user on 14 D
ecem
ber 2018
Hydrodynamic convection in accretion discs 4807
Figure 8. (a) Vertical profiles of the vertical heat (solid red line) and mass
flux (dashed black line). The outer vertical dashed (blue) lines mark the
boundaries of the convectively unstable region at initialization, while the
inner vertical dashed (red) lines mark the most convectively unstable points
at initialization. (b) Vertical profile of the xy-component of the Reynolds
stress tensor Rxy. In all cases, the results have been time-averaged between
orbits 5 and 10, spanning non-linear saturation.
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. Pthermal > Pmagnetic; cf. the middle-panel of Fig. 5 in
their paper).
Beyond |z| > H, the mass flux is directed away from the mid-
plane, peaking just beyond the point that 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.
Finally, in Fig. 8(b) we show the vertical profile of the hori-
zontally averaged Rxy taken from the 2563 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.3 Thus, the bulk of outward an-
gular 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.
3Note however that at any given instant in time, we generally do observe
regions where the Reynolds stress is negative.
Figure 9. Comparison of time evolution of volume-averaged Reynolds
stress taken from a simulation run in ATHENA (black) with the same quantity
taken from a simulation run in PLUTO (red).
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 vis-
cosity) 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 & Balbus’s convectively unsta-
ble profile (see Appendix A1). Explicit diffusion coefficients were
omitted and vertical periodic boundary conditions implemented.
Both simulations were run in a shearing box of size 4H × 4H ×
4H and at a resolution of 64 × 64 × 64. 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 em-
ployed 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 initial-
ization in both simulations were set to  = 10−3 and cs = 10−3,
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. Fig. 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 the 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
MNRAS 480, 4797–4816 (2018)
D
ow
nloaded from
 https://academ
ic.oup.com
/m
nras/article-abstract/480/4/4797/5064252 by U
niversity of C
am
bridge user on 14 D
ecem
ber 2018
4808 L. E. Held and H. N. Latter
Figure 10. Left: time evolution of volume-averaged Reynolds stress tensor taken from a PLUTO simulation run with the less diffusive HLLC solver. Right:
the same, but taken from a simulation run with the more diffusive HLL solver. The change in the sign of the stress tensor provides compelling evidence that
whether convection can transport angular momentum inwards or outwards depends on the diffusivity of the underlying numerical scheme.
solver, PPM reconstruction, and CTU integration, the Reynolds
stress is highly oscillatory in time in both codes, indicative, perhaps,
of numerical instability. The different behaviour (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 sensitivity of our results to the
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 Rxy > 0, i.e. out-
ward angular momentum transport. This is in disagreement with
the ZEUS results of SB96, who reported a Reynolds stress of Rxy
< 0. 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 Fig. 10: the
left-hand panel shows the time evolution of the volume-averaged
Rxy 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 sim-
ulations 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 2563, but with no change in the
sign of Rxy. 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 are summa-
rized in Table B2. They indicate that the sign of Rxy 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 Rxy > 0, while the more dif-
fusive Riemann solvers (HLL and a simple Lax–Friedrichs solver)
gave Rxy < 0. Altogether we explored 12 different configurations of
Riemann solver, interpolation scheme, and time-stepping method.
These ranged from the least diffusive set-up, which was identical
to that employed in the simulations described in the previous sec-
tion (a Roe Riemann solver, third-order-in-space WENO interpo-
lation, and third-order-in-time Runge–Kutta time-stepping), to the
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 compiled evidence that, as suspected, angular mo-
mentum transport due to convection can be sensitive to the diffusiv-
ity of the numerical set-up. Oversmoothing of the flow by diffusive
Riemann solvers such as HLL or by codes such as ZEUS, which
employs artificial viscosity and finite-differencing of the pressure
term, imposes a spurious axisymmetry on the flow. This enforced
axisymmetry ensures that the flow transports 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 (while retaining 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 & Ogilvie
(2010).
MNRAS 480, 4797–4816 (2018)
D
ow
nloaded from
 https://academ
ic.oup.com
/m
nras/article-abstract/480/4/4797/5064252 by U
niversity of C
am
bridge user on 14 D
ecem
ber 2018
Hydrodynamic convection in accretion discs 4809
We carry out a suite of simulations at a resolution mainly of
2563 investigating the effects on hydrodynamic convection when
the Rayleigh number Ra is increased, but the Prandtl number Pr is
fixed at 2.5.4 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 Ra = 105, 106, 107, 108, 109, and 1010. (Note
that the simulations undertaken at the two highest Ra are probably
underresolved, as explored later.) The Richardson number at ini-
tialization is fixed by the initial vertical profile, which is described
in detail in Section 4.1.1 (and shown in Fig. 1). For the profile pa-
rameters chosen, the Richardson number at initialization is Ri ∼
0.05. Further details of the simulations are given in Table B3 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 & Ogilvie (2010). We observe no instability for Ra = 105. At
Ra = 106 instability occurs but the flow appears relatively laminar
and axisymmetric; in particular, the Reynolds stress is negative
throughout the linear and non-linear phases. We conclude that the
critical Rayleigh number for the onset of convection lies in the range
105 < Rac < 106. At Ra = 107 the instability is more vigorous and
the flow field significantly more chaotic and non-axisymmetric in
the non-linear phase, at which point the Reynolds stress has become
positive. We conclude that the critical Ra at which the sign of Rxy
flips lies between 106 and 107. At higher Ra the flow appears even
more turbulent and non-axisymmetric. It is a relief that in the non-
linear 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 con-
siderably higher than in Lesur & Ogilvie (2010), who report a value
of Rac = 6900, but this is perhaps explained by their much larger
Richardson number (=0.2) that remains constant throughout their
box. On the other hand, the critical Ra at which Rxy 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 unambiguously out-
wards after non-linear saturation when Ra ≥ 107, it is almost al-
ways 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 dis-
agrees with our inviscid fiducial simulations (see the 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 of how the system should evolve: (a) in the linear
phase growing axisymmetric modes outcompete shearing waves,
(b) at sufficiently large amplitudes the modes are subject to sec-
ondary shear instabilities in the xz-plane that buckle the rising and
falling plumes, (c) perhaps concurrently or on a short time later
(and viscosity permitting), secondary non-axisymmetric instabili-
ties also attack the plumes because they exhibit significant shear in
the xy-plane as well, and (d) at this point, the flow degenerates into
something more disordered, and importantly, non-axisymmetric and
4We have also repeated some of the simulations at a Prandtl number of unity.
The differences with the Pr = 2.5 simulations are nominal.
Figure 11. Time evolution of volume-averaged xy-component of Reynolds
stress tensor from two simulations run at resolutions of 2563 (black curve)
and 5123 (red curve) at a Rayleigh number of Ra = 109. The vertical dashed
lines mark the end of the linear phase in the 2563 (black) and 5123 (red)
simulations, respectively.
the Reynolds stress flips sign. We conclude that viscosity prefer-
entially damps shearing waves vis-a-vis axisymmetric modes or
effectively kills off the artificial aliasing of shearing waves in invis-
cid 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 2563
points inadequate? We conducted simulations at Ra=108 with 2563
and 5123 grid points and found generally good agreement between
the two. The linear growth rates were almost identical and the
ultimate non-linear 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 2563
grid points are adequate to resolve a simulation with Ra=108.
Things start to deteriorate at a Rayleigh number of Ra = 109.
In Fig. 11, we plot the time evolution of the xy-component of the
volume-averaged Reynolds stress for two simulations run at 2563
(black) and 5123 (red). The lower resolution run possesses no ex-
tended period of negative Rxy, in contrast to the runs at Ra=108.
We speculate that at this resolution physical viscosity is subordi-
nate to the grid and non-axisymmetric disturbances are artificially
enhanced, probably via aliasing. At 5123, the stress is definitely
negative in the linear phase and there is a strong negative peak
shortly afterwards. The physical viscosity is now permitted to work
properly and appears to prohibit artificial non-axisymmetric distur-
bances. 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 non-
linear phase a little later in which the two flows closely resemble
each other. We conclude that in the linear phase 2563 is insufficient
to describe a simulation of Ra = 109, but in the later non-linear
phase it is probably adequate.
5 STRUCTURES IN FORCED C OMPRESS IBLE
C O N V E C T I O N
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
MNRAS 480, 4797–4816 (2018)
D
ow
nloaded from
 https://academ
ic.oup.com
/m
nras/article-abstract/480/4/4797/5064252 by U
niversity of C
am
bridge user on 14 D
ecem
ber 2018
4810 L. E. Held and H. N. Latter
Figure 12. Left: semilog plot of time evolution of volume-averaged vertical kinetic energy density in a simulation with thermal relaxation. Right: time
evolution of volume-averaged xy-component of the Reynolds stress tensor. The thermal relaxation time was set to τrelax = 4.2 −1.
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 mid-plane has no way to leave
the box, except through the distant vertical boundaries (in ZEUS nu-
merical 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 disc heats up to an in-
ordinate 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 τ relax. 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 time-scale τ relax then would be the characteristic time
that the radiative MRI system achieves thermodynamical quasi-
equilibrium. In our simulations, however, we take the relaxation
time-scale to be equal to the linear growth time of convection.
In addition, to mitigate the effects of mass outflows, we incor-
porate a mass source term to the simulations. At the end of the
nth 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 (n + 1)th time-step with the same vertical profile
used to initialize the density (cf. equation A3).
We implement the thermal relaxation term by making slight mod-
ifications 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 256 × 256
× 256 in a box of size 6H × 6H × 4H, corresponding to about
43 grid cells per scale height in the x- and y-directions and 64 grid
cells per scale height 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 σ = 0.2404 . Thus, the thermal relax-
ation time is τrelax = 4.2 −1, 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 ν = 1.075 × 10−5
corresponding to a Rayleigh number of Ra = 109 and a Prandtl
number of Pr = 2.5.
Finally, we have repeated the simulations in a cube of resolution
643 and size 4H × 4H × 4H, as well as at a resolution of 128 × 128 ×
64 in a box of size 6H × 6H × 4H, 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 Fig. 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
〈Ekin, z〉 ∼ 2 × 10−3 and the period remaining steady at t = 0.88
orbits.
Associated with the oscillations are large fluctuations in the xy-
component of the Reynolds stress tensor (shown in the right-hand
panel of Fig. 12). Instantaneous fluctuations in 〈Rxy〉 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 Rxy
∼ +9.9 × 10−6 and α ∼ +3.9 × 10−5, respectively. Furthermore,
comparing the oscillations in the vertical kinetic energy density
to the fluctuations in 〈Rxy〉 and also to the volume-averaged gas
pressure 〈P 〉 (not shown here) it is apparent that peaks in the kinetic
energy density are correlated with both 〈Rxy〉 < 0 and troughs in
〈P 〉, while troughs in the kinetic energy density are correlated with
〈Rxy〉 > 0 and peaks in 〈P 〉.
A clearer picture of the behaviour of the system emerges when
we study the structure of the flow during the cyclical phase. Fig. 13
MNRAS 480, 4797–4816 (2018)
D
ow
nloaded from
 https://academ
ic.oup.com
/m
nras/article-abstract/480/4/4797/5064252 by U
niversity of C
am
bridge user on 14 D
ecem
ber 2018
Hydrodynamic convection in accretion discs 4811
Figure 13. Snapshots at different times of the z-component of the velocity
in the xz-plane taken from a simulation in which convection was sustained
using thermal relaxation. Narrow convective cells shown in (a) just after
non-linear saturation merge to form large scale structures shown in (b),
which are destroyed (c) and recreated (d) in a cyclical manner with the
opposite rotation.
shows snapshots of the z-component of the velocity in the xz-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 λx ∼
5 x ∼ 0.177H , wherex is the size of grid a cell in the x-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 λx ∼ H, and by the time the quasi-steady equilibrium state sets
in (at around orbit 40) the radial wavelength of the convective cells
is λx ∼ 3.4H.
Comparing the snapshots of uz 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 α < 0. 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 Fig. 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 ±3H. 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 on to 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 ax-
isymmetric flow field, as confirmed by viewing slices of the pressure
perturbation δP in the xy-plane at different times. The simulation,
none the less, enters the cyclical phase during which this axisymme-
try is broken. As expected there is a flip in the sign of the Reynolds
stress from negative to positive. The behaviour 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 negli-
gibly on the numerical algorithm.
5.3 Zonal Flows
In Fig. 14, we plot, side-by-side, space–time diagrams in the xt-
plane of the yz-averaged pressure perturbation δP(x, t) ≡ (〈P〉yz
−〈P〉)/〈P 〉 and of the yz-averaged perturbation to the y-component
of the velocity 〈δuy〉yz(x, t) ≡ 〈uy − 1.5 x〉yz. It is immediately
evident from Fig. 14 that around the onset of the cyclical phase (orbit
20) alternating streaks in δP and δuy begin to develop. Respectively,
these mark alternating bands (in x) of high (δP > 0) and low (δP
< 0) pressure, and of super-Keplerian (δuy > 0) and sub-Keplerian
(δuy < 0) flow.
From Fig. 14, it is clear that the pressure and velocity perturba-
tions 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 ∂xP ′ ∼ ρ0 δuy) because, as is ev-
ident from the space–time diagrams, the flows are not stationary
in time but appear to fluctuate over about one orbit – the same
MNRAS 480, 4797–4816 (2018)
D
ow
nloaded from
 https://academ
ic.oup.com
/m
nras/article-abstract/480/4/4797/5064252 by U
niversity of C
am
bridge user on 14 D
ecem
ber 2018
4812 L. E. Held and H. N. Latter
Figure 14. Space–time diagrams of yz-averaged y-component of the perturbed velocity δuy (left) and of the pressure perturbation δP (right). (The perturbed
pressure is defined in the first paragraph of Section 5.3).
time-scale over which the large convective cells are created and de-
stroyed. Indeed, the creation and destruction of the zonal flows track
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 ob-
served while these structures remain coherent: axisymmetric con-
vective modes dominate during the lifetime of these structures and
transport angular momentum inwards.
In Fig. 15 (right-hand panels), we plot snapshots in the xy-
plane taken during the cyclical phase of the perturbation of the
y-component of the velocity δuy and of the perturbed pressure δP.
The axisymmetric structure of the zonal flows is clearly visible in
panels (b) and (d), as is the π /2 phase difference between the pres-
sure perturbation δP and the perturbed y-component of the velocity
δuy.
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 Fig. 16, we
plot a snapshot in the xy-plane taken during the cyclical phase of the
the z-component of the vertical component of vorticity ωz. Small
but coherent anticyclonic 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 anticyclonic 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 vor-
tices have aspect ratios approaching unity. They appear to have
limited three-dimensional extent, and we did not observe any vor-
tices 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 struc-
tures in recent simulations of rotating hydrodynamic convection
with uniform (rather than Keplerian) rotation and in the Boussi-
nesq (rather than compressible) regime (Julien et al. 2012, Rubio
et al. 2014, Favier, Silvers & Proctor 2014; Guervilly, Hughes &
Jones 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 (Ra ∼ 107–109), as we do. In addition,
Favier et al. (2014) and Guervilly et al. (2014) notice the forma-
tion of depth-invariant large-scale vortices in the xy-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 derived from uniformly rotating, Boussinesq
convection. Although the Rayleigh number of our simulation with
thermal relaxation was Ra = 109, which is consistent with the
Rayleigh numbers at which large-scale vortices were observed in
the Boussinesq simulations, we observe large-scale convective cells
(in the xz-plane) rather than large-scale vortices (in the xy-plane). We
do observe anticyclonic vortices in the xy-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 simula-
tions 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. Alternatively,
our large-scale axisymmetric cells could be interpreted as vortical
structures sheared out by the disk’s differential rotation.
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 asso-
ciated 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 verti-
cal convective transport, with the role of opacity in the simulations
of Hirose et al. (2014) being simply to modulate this cycle. Thus,
MNRAS 480, 4797–4816 (2018)
D
ow
nloaded from
 https://academ
ic.oup.com
/m
nras/article-abstract/480/4/4797/5064252 by U
niversity of C
am
bridge user on 14 D
ecem
ber 2018
Hydrodynamic convection in accretion discs 4813
Figure 15. Two-dimensional slices in the xy-plane of the perturbation to the y-component of the velocity δuy, and of the fractional pressure perturbation δP
from a simulation of forced compressible convection at z ≈ 0.5H. Left column: slices taken from a snapshot just after non-linear saturation (orbit 7.4). Right
column: slices taken from a snapshot generated during the cyclical phase (orbit 41.6).
our results demonstrate that strongly non-linear hydrodynamic tur-
bulent convection has a cyclical nature that might be generic, and
that might therefore be robust to the inclusion of more realistic
thermodynamics.
6 C O N C L U S I O N S
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
< 1 per cent, 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 discs, but also because they
may appear in some form during the non-linear phase of the evolu-
tion, especially on large scales, and during intermittent or cyclical
convection.
We then explored the non-linear 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
MNRAS 480, 4797–4816 (2018)
D
ow
nloaded from
 https://academ
ic.oup.com
/m
nras/article-abstract/480/4/4797/5064252 by U
niversity of C
am
bridge user on 14 D
ecem
ber 2018
4814 L. E. Held and H. N. Latter
Figure 16. Two-dimensional slice in the xy-plane of the z-component of
the vorticity ωz taken from a snapshot generated at the start of the cyclical
phase (orbit 20) and at z ≈ 0.5H. For clarity, we have zoomed in on the
upper-right quadrant of the xy-plane.
convection in a disc leads to inward (α < 0) or outward (α > 0)
angular momentum transport. In both forced and unforced convec-
tion, we found α > 0 in general in the non-linear phase. These
results were confirmed by a separate run using the code ATHENA,
but contradict the classical simulations of SB96, who reported in-
ward transport in both cases (forced and unforced) using the code
ZEUS.
This discrepancy reveals a set of unfortunate numerical difficul-
ties that complicate the simulation of convection in discs. These, in
large part, issue from the fact that the inviscid linear modes of con-
vection 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, imposes 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 viscosity (and thermal diffusivity), as this regularizes the
linear problem. We find that at a Richardson number of Ri ∼ 0.05,
onset of convection is observed for a critical Rayleigh number 105 <
Rac ≤ 106. Just above this value, convection is largely axisymmetric
and α < 0. At a larger second critical Ra between 106 and 107, the
sign of α switches and the flow becomes more turbulent and non-
axisymmetric. This sequence of states mirrors that simulated by
Lesur & 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 xz- and xy-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 convec-
tion 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 over-
all effect is to mimic the heating of the mid-plane and cooling of
the corona due to physical mechanisms that maintain the convec-
tively 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. Al-
though 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
α ∼ 10−6–10−5). 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 con-
vection in any of our simulations. By self-sustaining convection we
mean that (when α > 0) energy extracted from the shear by con-
vection might itself cause convective motions, which in turn extract
energy from the shear, closing the loop. It is more likely that if
convection is to occur in discs it will be as a by-product of other
processes, such as heating by density waves, emitted in the pres-
ence of a planet, or by dissipation of magnetorotational turbulence.
We intend to investigate both mechanisms and their instigation of
convection in future work.
AC K N OW L E D G E M E N T S
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 V., Amiez G., Gremaud P.-A., 1996, Commun. Numer. Methods
Eng., 12, 31
Armitage P. J., 2011, ARA&A, 49, 195
Aurnou J., Calkins M., Cheng J., Julien K., King E., Nieves D., Soderlund
K., Stellmach S., 2015, Phys. Earth Planet. Inter., 246, 52
Bell K. R., Lin D. N. C., 1994, ApJ, 427, 987
Bodo G., Cattaneo F., Mignone A., Rossi P., 2013, ApJ, 771, L23
Boley A. C., Durisen R. H., 2006, ApJ, 641, 534
Boyd J. P., 2001, Chebyshev and Fourier Spectral Methods. Dover, New
York USA
Cabot W., 1996, ApJ, 465, 874
Cannizzo J. K., 1993, in Wheeler J. C., ed., The Limit Cycle Instability In
Dwarf Nova Accretion Disks. World Scientific, Singapore, p. 6
Chiang E., Goldreich P., 1997, ApJ, 490, 368
Coleman M. S., Blaes O., Hirose S., Hauschildt P. H., 2018, ApJ, 857, 52
D’Alessio P., Canto¨ J., Calvet N., Lizano S., 1998, ApJ, 500, 411
Favier B., Silvers L. J., Proctor M. R. E., 2014, Physics of Fluids, 26, 096605
Goldreich P., Lynden-Bell D., 1965, MNRAS, 130, 125
MNRAS 480, 4797–4816 (2018)
D
ow
nloaded from
 https://academ
ic.oup.com
/m
nras/article-abstract/480/4/4797/5064252 by U
niversity of C
am
bridge user on 14 D
ecem
ber 2018
Hydrodynamic convection in accretion discs 4815
Guervilly C., Hughes D. W., Jones C. A., 2014, J. Fluid Mech., 758, 407
Hawley J. F., Gammie C. F., Balbus S. A., 1995, ApJ, 440, 742
Hirose S., 2015, MNRAS, 448, 3105
Hirose S., Blaes O., Krolik J. H., Coleman M. S., Sano T., 2014, ApJ, 787,
1
Julien K., Knobloch E., Rubio A. M., Vasil G. M., 2012, PRL, 109, 254503
Klahr H., Henning T., Kley W., 1999, ApJ, 514, 325
Kley W., Papaloizou J., Lin D., 1993, ApJ, 416, 679
Latter H. N., Papaloizou J., 2017, MNRAS, 472, 1432
Lesur G., Kunz M. W., Fromang S., 2014, A&A, 56
Lesur G., Ogilvie G. I., 2010, MNRAS, 404, L64
Lin D., Papaloizou J., 1980, MNRAS, 191, 37
Lin D., Papaloizou J., Kley W., 1993, ApJ, 416, 689
Lyra W., Richert A. J., Boley A., Turner N., Mac Low M.-M., Okuzumi S.,
Flock M., 2016, ApJ, 817, 102
Mignone A., Bodo G., Massaglia S., Matsakos T., Tesileanu O., Zanni C.,
Ferrari A., 2007, ApJS, 170, 228
Mignone A., Flock M., Stute M., Kolb S., Muscianisi G., 2012, A&A, 545,
A152
Rubio A. M., Julien K., Knobloch E., Weiss J. B., 2014, PRL, 112, 144501
Ruden S. P., Papaloizou J., Lin D., 1988, ApJ, 329, 739
Ryu D., Goodman J., 1992, ApJ, 388, 438
Spruit H. C., Nordlund A., Title A., 1990, ARA&A, 28, 263
Stone J. M., Balbus S. A., 1996, ApJ, 464, 364
Stone J. M., Gardiner T. A., 2010, ApJS, 189, 142
Stone J. M., Gardiner T. A., Teuben P., Hawley J. F., Simon J. B., 2008,
ApJS, 178, 137
Zingale M. et al., 2002, ApJS, 143, 539
APPEN D IX A : C ONVECTIVELY UNSTA BL E
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 con-
vectively unstable disc profiles that satisfy hydrostatic equilibrium,
are convectively unstable, and can conveniently initialize simula-
tions.
A1 Stone and Balbus (1996) profile
The SB96 profile employs a power-law profile in temperature, see
equation (9). For p = 3/2, the density can be obtained analytically
ρ = ρ0(1 − s3)−(1+g/3)(1 − s)g exp
{
2g
[
s − 1√
3
tan−1
( √
3s
s + 2)
)]}
,
where ρ0 is the mid-plane density, s = (z/H)1/2f1/3, and g = f−4/3, in
which f = H3/2A/T0, and H = cs/ is the mid-plane scale height. The
pressure is obtained from the ideal gas equation of state. Note that
for |z| > (T0/A)1/p we must have vacuum, which ties the maximum
numerical domain to the ratio T0/A.
The buoyancy frequency for p = 3/2 is given by
N2B
2
= 1
1 − ( z
H
)3/2
f
[(
1 − 1
γ
)( z
H
)2
− 3
2
( z
H
)3/2
f
]
, (A1)
which corresponds to a profile that is negative (convectively unsta-
ble) within some region |Lc| < 0 about the mid-plane and positive
(convectively unstable) outside of this region. The width of the
convectively unstable region is given by
|Lc| =
[
3
2
(
1 − 1
γ
)−1
H 2
A
T0
]2
. (A2)
Although convenient within a limited choice of parameters, the
SB96 profile suffers from the drawback that the width of the con-
vectively unstable region is sensitive to the size of the box through
the ratio T0/A. Increasing the vertical box size necessarily decreases
the size of the convectively unstable region.
A2 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.
equation (8), provided such a profile.
The associated density is
ρ = ρ0ez2/βH 2 exp
(
−βH
22
2T0
ez
2/βH 2
)
, (A3)
where ρ0 is mid-plane density and H is the mid-plane scale height.
Pressure is obtained from the ideal gas equation of state, as above.
The buoyancy frequency is given by
N2B
2
= 2
βH 2
z2
[(
(1 − 1
γ
)
βH 22μ
2T0R e
z2/βH 2 − 1
]
. (A4)
The boundary of the convectively unstable region is hence given by
|Lc| ≡
√√√√βH 2 ln
[(
1 − 1
γ
)−1 2T0R
βμH 22
]
. (A5)
APPENDI X B: TABLES OF SI MULATI ONS
MNRAS 480, 4797–4816 (2018)
D
ow
nloaded from
 https://academ
ic.oup.com
/m
nras/article-abstract/480/4/4797/5064252 by U
niversity of C
am
bridge user on 14 D
ecem
ber 2018
4816 L. E. Held and H. N. Latter
Table B1. Simulations of unforced compressible hydrodynamic convection in the shearing box carried out in PLUTO. No explicit viscosity or heat-
ing/cooling/thermal relaxation was included. The box size in all simulations is 4H × 4H × 4H. The simulations were initialized either with a Gaussian
temperature convectively unstable vertical disc profile (see Appendix A2) or with the convectively unstable vertical disc profile used in Stone & Balbus (1996)
(SB96; see Appendix A1). 〈〈.〉〉t denotes a time- and volume-averaged quantity. Ekin, z is the vertical kinetic energy density, which we have measured both at
non-linear saturation (where it is maximum) and at the end of the simulation. ρ is the percentage change in density over the course of the simulation, Rxy is
the Reynolds stress, α is the alpha viscosity parameter, and σ is the linear phase growth rate given in units of . All simulations were run for 38 orbits.
Run Resolution Profile max(〈Ekin, z〉) 〈Ekin, z〉|end  (〈ρ〉) (%) 〈〈Rxy〉〉t 〈〈α〉〉t σ ()
NSTR22c01 643 Gaussian Temp. 1.7 × 10−4 2.0 × 10−6 − 1.4 +4.8 × 10−6 +2.3 × 10−5 0.2404
NSTR22c02 1283 Gaussian Temp. 1.1 × 10−4 3.7 × 10−6 − 3.4 +2.7 × 10−6 +1.5 × 10−5 0.2707
NSTR22c03 2563 Gaussian Temp. 7.3 × 10−5 9.8 × 10−6 − 8.8 +1.7 × 10−6 +9.6 × 10−6 0.2853
NSTR22c04 5123 Gaussian Temp. 6.4 × 10−5 1.9 × 10−5 − 12.7 +1.3 × 10−6 +7.1 × 10−6 0.2811
NSTR21c01 643 SB96 9.4 × 10−4 9.0 × 10−6 − 2.1 +1.9 × 10−5 +2.4 × 10−5 0.2059
NSTR21c03 2563 SB96 3.6 × 10−6 3.9 × 10−5 − 10.3 +6.0 × 10−6 +7.7 × 10−6 0.2424
Table B2. Simulations of unforced compressible hydrodynamic convection in the shearing box with different numerical schemes in PLUTO. No explicit
viscosity or heating/cooling/thermal relaxation was included. All simulations were run in a box of size 4H × 4H × 4H, and initialized with the Gaussian
temperature profile discussed in Appendix A2. The time-averages of the volume-averaged Reynolds stress 〈Rxy〉 and of the volume-averaged angular momentum
transport parameter 〈α〉 have been taken over an interval spanning non-linear saturation. All simulations were run for 38 orbits.
Run Resolution Solver Interpolation
Time-
stepping 〈〈Rxy〉〉t 〈〈α〉〉t Comments
NSTR22e01a 643 HLLC WENO3 RK3 +4.5 × 10−6 +2.5 × 10−5 –
NSTR22e01c 2563 HLLC WENO3 RK3 +1.8 × 10−6 +1.0 × 10−5 –
NSTR22e02a 643 HLL WENO3 RK3 −1.3 × 10−7 −7.4 × 10−7 |δu| = O(10−3)
NSTR22e02c 2563 HLL WENO3 RK3 −1.1 × 10−5 −6.1 × 10−5 –
NSTR22e03a 643 TVDLF WENO3 RK3 – – No instability
NSTR22e11a 643 HLLC Linear TVD Hancock +3.4 × 10−6 +2.2 × 10−5 –
NSTR22e11i 2563 HLLC Linear TVD Hancock +1.5 × 10−6 +9.2 × 10−6 UMIST limiter
NSTR22e12a 643 HLL Linear TVD Hancock −2.2 × 10−5 −1.4 × 10−4 –
NSTR22e12b 1283 HLL Linear TVD Hancock −1.6 × 10−5 −1.0 × 10−4 –
NSTR22e12c 2563 HLL Linear TVD Hancock −9.5 × 10−6 −6.2 × 10−5 –
NSTR22e12i 2563 HLL Linear TVD Hancock −1.2 × 10−5 −8.2 × 10−5 UMIST limiter
NSTR22e13a 643 TVDLF Linear TVD Hancock −3.5 × 10−5 −2.3 × 10−4 –
NSTR22e13b 1283 TVDLF Linear TVD Hancock −1.4 × 10−5 −9.4 × 10−5 –
NSTR22e21a 643 HLLC Linear TVD RK2 +4.0 × 10−6 +2.6 × 10−5 –
NSTR22e22a 643 HLL Linear TVD RK2 −2.2 × 10−5 −1.3 × 10−4 –
NSTR22e22b 1283 HLL Linear TVD RK2 −2.0 × 10−5 −1.2 × 10−4 –
NSTR22e23a 643 TVDLF Linear TVD RK2 −1.9 × 10−5 −1.2 × 10−4 –
NSTR22e23b 1283 TVDLF Linear TVD RK2 −1.3 × 10−5 −7.6 × 10−4 –
Table B3. Simulations of unforced compressible hydrodynamic convection in the shearing box with explicit kinematic viscosity ν and thermal diffusivity χ
included. No thermal relaxation was included. All simulations were run in a box of size 4H × 4H × 4H. 〈〈α〉〉t|linear is the volume-averaged alpha viscosity
parameter time-averaged over the linear phase, 〈〈α〉〉t|NL is the same quantity time-averaged over the non-linear phase, and min(〈$Rxy〉) is the minimum value
of the xy-component of the volume-averaged Reynolds stress.
Run Resolution Instability? ν χ Ra Ri Pr 〈〈α〉〉t|linear 〈〈α〉〉t|NL min(〈Rxy〉)
NSTR22Ra1 2563 N 1.075 × 10−3 4.300 × 10−4 105 0.05 2.5 – – –
NSTR22Ra2 2563 Y 3.400 × 10−4 1.360 × 10−4 106 0.05 2.5 −5.1 × 10−6 −3.9 × 10−5 −2.1 × 10−5
NSTR22Ra3 2563 Y 1.075 × 10−4 4.300 × 10−5 107 0.05 2.5 −3.5 × 10−6 +3.9 × 10−6 −2.7 × 10−5
NSTR22Ra4a 2563 Y 3.400 × 10−5 1.360 × 10−5 108 0.05 2.5 −1.0 × 10−6 +1.4 × 10−5 −1.1 × 10−5
NSTR22Ra4b 5123 Y 3.400 × 10−5 1.360 × 10−5 108 0.05 2.5 −2.8 × 10−6 +7.1 × 10−6 −1.4 × 10−5
NSTR22Ra5a 2563 Y 1.075 × 10−5 4.300 × 10−6 109 0.05 2.5 −3.5 × 10−8 +1.5 × 10−5 −9.9 × 10−8
NSTR22Ra5b 5123 Y 1.075 × 10−5 4.300 × 10−6 109 0.05 2.5 −1.6 × 10−7 +1.3 × 10−5 −3.4 × 10−6
NSTR22Ra6a 2563 Y 3.400 × 10−6 1.360 × 10−6 1010 0.05 2.5 +7.6 × 10−8 +1.5 × 10−5 +2.2 × 10−13
NSTR22Ra6b 5123 Y 3.400 × 10−6 1.360 × 10−6 1010 0.05 2.5 −1.1 × 10−8 +1.6 × 10−5 −4.9 × 10−8
This paper has been typeset from a TEX/LATEX file prepared by the author.
MNRAS 480, 4797–4816 (2018)
D
ow
nloaded from
 https://academ
ic.oup.com
/m
nras/article-abstract/480/4/4797/5064252 by U
niversity of C
am
bridge user on 14 D
ecem
ber 2018