MNRAS 471, 317–336 (2017) doi:10.1093/mnras/stx1548
Advance Access publication 2017 June 21
Gravitoturbulence and the excitation of small-scale parametric
instability in astrophysical discs
A. Riols,1‹ H. Latter1 and S.-J. Paardekooper1,2
1Department of Applied Mathematics and Theoretical Physics, Centre for Mathematical Sciences, University of Cambridge,
Wilberforce Road, Cambridge CB3 0WA, UK
2Astronomy Unit, School of Physics and Astronomy, Queen Mary, University of London, Mile End Road, London E1 4NS, UK
Accepted 2017 June 19. Received 2017 June 2; in original form 2017 April 5
ABSTRACT
Young protoplanetary discs and the outer radii of active galactic nuclei may be subject to
gravitational instability and, as a consequence, fall into a ‘gravitoturbulent’ state. While in this
state, appreciable angular momentum can be transported; alternatively, the gas may collapse
into bound clumps, the progenitors of planets or stars. In this paper, we numerically characterize
the properties of 3D gravitoturbulence, focusing especially on its dependence on numerical
parameters (resolution, domain size) and its excitation of small-scale dynamics. Via a survey
of vertically stratified shearing-box simulations with PLUTO and RODEO, we find (a) evidence
that certain gravitoturbulent properties are independent of horizontal box size only when
the box is larger than 40H, where H is the scaleheight, (b) at high resolution, small-scale
isotropic turbulence appears off the mid-plane around z  0.5–1H and (c) this small-scale
dynamics results from a parametric instability, involving the coupling of inertial waves with
a large-scale axisymmetric epicyclic mode. This mode oscillates at a frequency close to 
and is naturally excited by gravitoturbulence via a non-linear process to be determined. The
small-scale turbulence we uncover has potential implications for a wide range of disc physics,
e.g. turbulent saturation levels, fragmentation, turbulent mixing and dust settling.
Key words: accretion, accretion discs – instabilities – turbulence – protoplanetary discs.
1 IN T RO D U C T I O N
Gravitational instability (GI) attacks gaseous accretion discs that are
sufficiently cold and massive, as might be the case in the early stages
of a protoplanetary (PP) disc’s life (e.g. Durisen et al. 2007), or be-
yond roughly 0.01 pc in active galactic nuclei (AGN; e.g. Shlosman
& Begelman 1987). Recent observations of wide-orbit planets and
spiral structure have fuelled ongoing interest in GI within the PP
disc context (Kalas et al. 2008; Fukagawa et al. 2013; Christiaens
et al. 2014), while indications of transonic turbulence from maser
emission, and the long-standing questions of disc truncation, star
formation and the sustenance of high accretion rates motivate its
study in AGN (e.g. Wallin, Watson & Wyld 1998; Goodman 2003).
The critical parameter that sets the onset of the GI is the Toomre
Q (Toomre 1964),
Q = csκ
πG0
, (1)
where cs is the sound speed, κ the epicyclic frequency and 0
the background surface density. In a razor thin disc, linear ax-
isymmetric disturbances are unstable when Q < 1, but non-linear
 E-mail: ar764@cam.ac.uk
non-axisymmetric instability can occur for Q  1. If the cooling is
relatively inefficient, simulations show that the disc falls into a grav-
itoturbulent state comprised of a disordered field of spiral density
waves that can transport significant angular momentum (Gammie
2001; Rice et al. 2014). On the other hand, if the cooling is overly
efficient, the disc fragments into dense clumps that may be the
precursors of gaseous giant planets or stars (Cameron 1978; Boss
1997).
Most numerical work is undertaken in a global setting, and indeed
the characteristic length-scales of the gravitoturbulent spiral waves
are generally large. But owing to inadequate resolution, global mod-
els struggle to describe dynamics on scales of the order of the disc
thickness H or less, which are important for turbulent mixing, dust
settling and (possibly) fragmentation. Local 3D shearing-box mod-
els are better suited for exploring these processes, though this is a
research direction that has been neglected, for the most part, in the
modelling of GI in discs. Recently, Shi & Chiang (2014) presented
3D local box simulations that showed that the gravitoturbulent state
failed to exhibit appreciable structure on scales less than H (see also
Hirose & Shi 2017). However, their resolution (a maximum eight
points per H) was probably too low to properly characterize any
small-scale features. In particular, it was unclear if their results had
converged with resolution. While the density waves are expected to
C© 2017 The Authors
Published by Oxford University Press on behalf of the Royal Astronomical Society
318 A. Riols, H. Latter and S.-J. Paardekooper
exhibit little vertical structure (being essentially f modes; Ogilvie
1998), they are potentially vulnerable to parasitic instabilities rely-
ing on a parametric resonance between small-scale inertial modes
(Fromang & Papaloizou 2007; Bae, Nelson & Hartmann 2016b).
An open question is whether this secondary instability also operates
in gravitoturbulence, and what influence it has on the properties of
the spiral waves.
To address this issue, we performed 3D local shearing-box sim-
ulations of gravitoturbulent discs with vertical stratification and an
ideal gas law. Two different codes were used, PLUTO and RODEO,
and these were subject to a number of numerical checks, so as
to confirm the trustworthiness of our results. Our main finding is
that vigorous and isotropic small-scale turbulence develops slightly
off the disc mid-plane and disturbs the spiral waves generated by
GI. The small-scale turbulence is obtained only in high-resolution
runs (at least 20 cells per disc scaleheight), and thus has not been
captured in previous under-resolved simulations: either smoothed
particle hydrodynamics or grid-based, local or global.
We are persuaded that the small-scale activity is excited via a
3D parametric instability that couples pairs of inertial waves and
a large-scale axisymmetric epicyclic oscillation. Unstable inertial
waves develop in the mid-plane but break non-linearly at a higher
altitude, resulting in helical, incompressible and non-axisymmetric
motions around z = 0.5–1 H. We also consider the merits of the
Kelvin–Helmholtz instability, vertical splashing caused by colliding
density wakes and the preferential breaking of wakes in the disc
atmosphere, but find each in conflict with our numerical diagnostics.
The provenance of the large-scale epicyclic oscillations (upon which
the parametric instability feeds) is discussed, and we speculate on
their importance to both the subcritical transition to GI turbulence
and its continued sustenance.
The structure of the paper is as follows. First, in Section 2, we
introduce the basic equations of the problem and discuss our nu-
merical setup. In Section 3, we study the dependence of the 3D
gravitoturbulent state on the numerical parameters, in particular
the grid resolution and box size. We show evidence that turbu-
lent quantities are converged with horizontal box sizes Lx = Ly if
Lx  40H0. In Section 4, we characterize the small-scale dynamics
and the large-scale axisymmetric oscillation in Fourier space. In
Section 5, we analyse the nature of the small-scale parasitic tur-
bulence and show how it may result from a parametric instability.
Finally, in Section 6, we discuss the possible implications of our
results for disc physics, in particular for turbulent mixing, dust dy-
namics and the evolution of magnetic fields.
2 N U M E R I C A L M O D E L
2.1 Governing equations
We assume that the gas orbiting around the central object is ideal,
its pressure P and density ρ related by γP = ρc2s , where cs is
the sound speed and γ the ratio of specific heats. The pressure is
related to internal energy U by P = (γ − 1)U. We neglect molecular
viscosity. We treat a local Cartesian model of an accretion disc (the
shearing sheet; Goldreich & Lynden-Bell 1965). In this model,
the axisymmetric differential rotation is approximated locally by
a linear shear flow and a uniform rotation  =  ez, with S =
(3/2) for a Keplerian equilibrium. We denote by (x, y, z) the
shearwise, streamwise and spanwise directions, corresponding to
radius, azimuth and the vertical. The shearing sheet approximation
works best when the disc’s angular thickness is small (see for e.g.
Balbus & Papaloizou 1999), and thus PP discs are only marginally
covered by the model. In fact, for our largest domains, it is difficult
to justify the shearing box. AGN discs, being much thinner, are
better described. We persist with the shearing sheet, however, as it
provides a well-defined platform and the superior resolution needed
to probe the smaller scales.
The evolution of density ρ, total velocity field v and internal
energy U follows:
∂ρ
∂t
+ ∇ · (ρv) = 0, (2)
∂v
∂t
+ v · ∇v + 2 × v = −∇ − ∇P
ρ
, (3)
∂U
∂t
+ ∇ · (Uv) = −P∇ · v − U
τc
. (4)
The total velocity field may be decomposed into
v = −Sx ey + u, (5)
where u is the perturbed velocity field. The potential  is the sum
of the tidal potential induced by the central object in the local frame,
c = 122z2 − 322 x2, and the gravitational potential induced by
the disc itself, s, which obeys the Poisson equation
∇2s = 4πGρ. (6)
In the energy equation (4), the cooling varies linearly with U with
a typical time-scale τ c referred to as the ‘cooling time’. This pre-
scription is not very realistic but allows us to control the rate of
energy loss via a single parameter. We neglect thermal conductiv-
ity. Finally, −1 defines our unit of time and H0, the initial disc
scaleheight (defined below), our reference unit of length.
2.2 Disc background equilibrium
When τ c → ∞, the governing equations admit an equilibrium
state characterized by u = 0, and vertical hydrostatic balance. We
consider homentropic equilibria, for which the vertical profile is
polytropic P = Kργ . Here K = c2s0/(γρ
γ−1
0 ), with cs0 and ρ0 de-
noting the mid-plane sound speed and density, respectively. The
equilibrium equations are
K
[
1
ρ
dργ
dz
]
+ z2 + ds
dz
= 0, (7)
d2s
dz2
= 4πGρ, (8)
which reduce to an inhomogeneous form of the Emden–Fowler
equation. Appendix A gives details on how to solve these equations
numerically. These solutions form part of our simulations’ initial
condition.
2.3 Numerical methods
2.3.1 Codes
Direct numerical simulations of the three-dimensional flow are per-
formed in the shearing box. The box has a finite domain of size
(Lx, Ly, Lz), discretized on a mesh of (NX, NY, NZ) grid points. For
most of the numerical runs, we use the Godunov-based PLUTO code
(Mignone et al. 2007), which is well adapted to highly compress-
ible flow and shocks. This scheme uses a conservative finite-volume
MNRAS 471, 317–336 (2017)
GI and turbulence in astrophysical discs 319
method that solves the approximate Riemann problem at each inter-
cell boundary. It is known to successfully reproduce the behaviour
of conserved quantities like mass, momentum and total energy. The
Riemann problem is handled by the HLLC solver which properly
describes contact discontinuities and has the advantage of being
robust and positivity preserving.
To allow longer time-steps and eliminate numerical artefacts at
the boundaries, where the background shear flow is often very
strong, we use the orbital advection algorithm of PLUTO. It is based
on splitting the equation of motion into two parts, the first containing
the linear advection operator due to the background Keplerian shear
and the second the standard fluxes and source terms. Finally, PLUTO
conserves the total energy and so the heat equation is not solved
directly as in equation (4). Hence, the code captures the irreversible
heat produced by shocks due to numerical diffusion, consistent with
the Rankine–Hugoniot conditions.
In Section 3.5, we compare our results with another code, RODEO,
that has previously simulated self-gravitating disc in the 2D shearing
box (Paardekooper 2012). Similar to PLUTO, it is a Godunov-based
method but one that relies on the Roe solver (Roe 1981) to calcu-
late interface fluxes. Second-order accuracy in space and time is
achieved through a wave limiting procedure as described, for exam-
ple, in LeVeque (2002). Most of the results were obtained with the
minmod wave limiter (Roe 1986), although we briefly explore dif-
ferent limiters. One other significant difference with PLUTO is that the
RODEO simulations were performed in shearing coordinates, where
the usual y coordinate is replaced by y′ = y + 3/2tx. In this way,
the use of an orbital advection algorithm is avoided, at the expense
of a periodic remap every tremap, which, by a clever choice of tremap,
can be done by shifting an integer number of grid cells in y. This
procedure therefore does not introduce any numerical diffusion.
2.3.2 Poisson solver
While in RODEO a pure Fourier Poisson solver was used (Koyama &
Ostriker 2009), a different approach was taken in the PLUTO simula-
tions in order to more robustly establish our numerical results. This
is explained in this subsection.
To compute the gravitational potential, we take advantage of the
shear-periodic boundary conditions. In a similar manner to Riols &
Latter (2016), we first shift back the density in y to the time it was
last periodic (t = tp). Then, for each plane of altitude zk, we perform
the direct 2D Fourier transform of the density ρ and obtain a ver-
tical profile for each coefficient ρˆkx ,ky (z) of its Fourier decomposi-
tion, where kx and ky are radial and azimuthal wavenumbers. Using
equation (6), it is straightforward to show that in Fourier space,
the coefficients of the gravitational potential satisfy the Helmholtz
equation[
d2
dz2
− k2
]
ˆkx,ky (z) = 4πGρˆkx ,ky (z) (9)
with k = k2x + k2y and ˆkx,ky the Fourier transform of the poten-
tial. This differential equation is solved in the complex plane by
means of a fourth-order finite-difference scheme and a direct in-
version method. The discretized system takes the form of a linear
problem AX = B, where X is a vector representing the discretized
potential’s z profile, A is a pentadiagonal matrix and B is a column
vector containing the right-hand side of equation (9) and extra co-
efficients setting the boundary conditions. We use a fast algorithm
involving O(NZ) flops to invert the matrix and obtain the discretized
coefficients ˆkx,ky (zk). For each altitude zk, we finally compute the
inverse Fourier transform of the potential and shift it back to the
initial sheared frame. Note that gravitational forces are obtained by
computing the derivative of the potential in each direction, using a
fourth-order finite-difference method.
In total, the computational cost is of order O(Nlog (NXNY)), where
N = NXNYNZ, and is hence less expensive than a full 3D Fourier
decomposition that would be O(Nlog (N)). Moreover, methods us-
ing vertical Fourier decomposition generally assume periodic or
vacuum boundary conditions for the potential. Hence, a correct
treatment of self-gravity in these methods requires the addition of
two buffer zones of size Lz/2 on either sides of the box, greatly
increasing the numerical workload (Koyama & Ostriker 2009; Shi
& Chiang 2014). In contrast, our setup can handle any boundary
condition without artificially augmenting the vertical domain. The
stratified disc equilibria as well as the linear stability of these equi-
libria has been tested to ensure that our implementation is correct
(see Appendices A and B).
2.3.3 Boundary conditions
The shearing-box framework implicitly assumes periodic boundary
conditions in y and shear-periodic boundary condition in x. The
most delicate part is to assign suitable vertical boundary conditions.
We use the standard outflow conditions for velocity and density
fields but compute a hydrostatic balance in the ghost cells for pres-
sure, taking into account the large-scale vertical component of the
self-gravity (averaged in x and y). In this way, we significantly re-
duce the excitation of waves near the boundary (see Fig. A2 in
Appendix A).
For the gravitational potential, we impose
d
dz
ˆkx,ky (±Lz/2) = ∓k ˆkx,ky (±Lz/2). (10)
This condition is an approximation of the Poisson equation in
the limit of low density. Finally, we impose a density floor of
10−4 /H0, which prevents the time-step getting too small because
of evacuated regions near the vertical boundaries.
2.4 Simulation setup
2.4.1 Box size and resolution
The axisymmetric linear theory tells us that the fastest growing
mode possesses a radial length-scale of the order of H Q when
Q < 1. Although our simulations focus on the regime Q  1, we
expect typical length-scales for spiral waves to be  H. So in order
to obtain good statistical averages of the turbulent properties and
ensure that the largest structures remain smaller than the domain
size, it is necessary to set Lx, Ly 	 H. As a compromise between
this constraint and numerical feasibility (for sufficient resolution),
we employ boxes of intermediate size Lx = Ly = 20 H0, though we
explore larger sizes in Section 3.2.
In most of the PLUTO simulations, we simulate a symmetric flow
with respect to the mid-plane. Thus, the vertical domain of the
box only extends from the mid-plane to 3 H0. We checked that anti-
symmetric modes do not affect at all the properties of the turbulence
and the results of the present paper.
According to previous shearing-box simulations (Gammie 2001;
Shi & Chiang 2014), a minimum resolution of four to five grid cells
per H0 in the horizontal directions is required to correctly capture
the onset of non-linear instability. However, to properly determine
MNRAS 471, 317–336 (2017)
320 A. Riols, H. Latter and S.-J. Paardekooper
certain properties of the ensuing turbulence, as well as the fragmen-
tation criterion, more may be required. For example, Paardekooper
(2012) showed that in 2D, the critical cooling time τ c for fragmen-
tation is still dependent on resolution when the latter exceeds 40
points per H0. In Section 3.2, we compare different gravitoturbulent
states obtained at different resolutions, ranging from 3.2 points per
H0 to almost 26 points per H0 in the horizontal directions. In the
vertical direction, we normally used 64 points in total over 3H0.
2.4.2 Initial conditions and simulations parameters
In all our simulations, we use a fixed heat capacity ratio γ = 5/3
and an initial mid-plane sound speed cs0 = H0/ = 1. The initial
conditions are similar to those used in Riols & Latter (2016), except
that we must stipulate vertical profiles for density and pressure.
First, for an initial Toomre parameter Q0 slightly larger than
1, we compute the vertical density and pressure profile associated
with a homentropic and self-gravitating disc equilibrium with no
cooling (see equations 7, 8 and Appendix A). Second, random
non-axisymmetric density and velocity perturbations are imposed
upon this equilibrium. They possess a finite amplitude decreasing
exponentially with altitude. These initial fluctuations are intensified
by the instability and break down into a turbulent flow after a short
period of time t  10 – 30−1.
In order to prevent the disc from fragmenting in the early stages,
cooling is only introduced once the average turbulent quantities
have converged to a fixed value. This also provides a relatively
‘soft landing’ upon the gravitoturbulent state, one that makes fewer
demands on the numerical scheme.
If mass is lost through the vertical boundary, it is replenished
near the mid-plane. The distribution of mass added to the disc at
each time-step exhibits a Gaussian profile ∝ exp[−z2/(2H 20 )]. The
mass-injection rate varies in time so that 0 = 1 remains constant
during the simulation. We checked that the amount of mass injected
at each orbital period is negligible compare to the total mass and
does not affect the results.
2.5 Diagnostics
2.5.1 Averaged quantities
To analyse the statistical behaviour of the turbulent flow, we define
two different volume averages of a quantity X. The first one is the
standard average
〈X〉 = 1
LxLyLz
∫
V
X dV . (11)
The second is the density-weighted average, defined by
〈X〉w =
∫
V
ρX dV∫
V
ρ dV
. (12)
An important quantity that characterizes self-gravitating discs is the
average 2D Toomre parameter defined by
Q2D = 〈cs〉w 
πG 〈〉 , (13)
where 〈〉 = Lz〈ρ〉 is the average surface density of the disc. To
simplify notation, we denote Q = Q2D.
Another quantity that characterizes the turbulent dynamics is the
coefficient α, which measures the angular momentum transport.
This quantity is related to the average Reynolds stress Hxy and
gravitational stress Gxy,
α = 2
3γ 〈P 〉
〈
Hxy + Gxy
〉
, (14)
where
Hxy = ρuxuy and Gxy = 14πG
∂
∂x
∂
∂y
.
The radial flux of angular momentum gives rise to the only source
of energy in the system that can balance the cooling. This energy,
initially in the form of kinetic energy, is irremediably converted into
heat by turbulent motions.
Lastly, in order to study the energy budget of the flow, we intro-
duce the average kinetic and internal energy denoted respectively
by
Ec = 12 〈ρu
2〉, U = (γ − 1)〈P 〉.
2.5.2 Fourier decomposition
In Section 4, we investigate the flow in Fourier space. Any field F
can be decomposed in the following way:
F =
NX/2∑
=−NX/2
NY /2∑
m=−NY /2
ˆF,m(z, t) exp
[
i(kx(t)x + kyy)
]
, (15)
with
kx(t) = kx0 +
3
2
mky0 t and ky = mky0 . (16)
The Lagrangian radial wavenumber kx0 allows us to define and
label a given ‘shearing wave’ , if ky is known. This quantity rep-
resents the wavenumber that a mode would have in a system of
coordinates following the shear. In the fixed coordinate system,
however, this wave has an Eulerian wavenumber kx(t) that increases
with time. The wave is referred to as ‘leading’ when kykx(t) < 0 and
‘trailing’ when kykx(t) > 0.
In practice, to compute the FFT transform of a given field, we first
have to express this field (in real space) in a system of coordinates
comoving with the shear, in a manner similar to our calculation
of the self-gravitating potential (see Section 2.3.2). The change of
variables is
y −→ y ′ − 3
2
x(t − tp), (17)
with tp the time corresponding to the nearest periodic point. This
process enforces strict periodicity upon the field. We then compute
the standard FFT and obtain, for each altitude z, a spectral map in 
and m. If we consider an interval of time [ntp, (n + 1)tp], a point in
this map represents the amplitude of a given wave  that has initially
kx(ntp) = kx0 . However, this representation is not valid for larger
time, because the radial wavenumber of the wave has increased
during that time (especially for large ky) and one needs to remap
the 2D Fourier spectrum on to a grid in Eulerian wavenumbers
(kx(t), ky).
3 M E A N T U R BU L E N T PRO P E RT I E S A N D
N U M E R I C A L C O N V E R G E N C E
In this section, we analyse the properties of 3D gravitoturbulence
and its dependence on resolution and box size for a fixed cooling
time τ c = 20−1. We then move on to explore the effects of different
initial conditions, numerical methods and cooling times.
MNRAS 471, 317–336 (2017)
GI and turbulence in astrophysical discs 321
Table 1. Parameters and properties of runs in a box of Lx = Ly = 20H0. The third column indicates the time over which quantities have been averaged
(not including the transient phase). Q is the average Toomre parameter, Ec and U are the box- and time-averaged kinetic and internal energy, Hxy, Gxy
and α are respectively the averaged Reynolds, gravitational stress and transport efficiency.
Run Resolution Time (in −1) τ c Q Ec U Hxy Gxy α
PL20-64 64 × 64 × 32 500 20 1.14 0.047 0.332 0.001 98 0.007 51 0.0176
PL20-128 128 × 128 × 64 3000 20 1.24 0.105 0.396 0.004 11 0.009 15 0.020 05
PL20-256 256 × 256 × 64 100 20 1.25 0.115 0.427 0.005 0.0109 0.022
PL20-512 512 × 512 × 64 200 20 1.26 0.109 0.412 0.005 94 0.009 63 0.0226
PL20-128-b 128 × 128 × 64 400 20 1.22 0.0762 0.377 0.003 76 0.009 77 0.0212
PL20-256-RK3 256 × 256 × 64 200 20 1.23 0.0741 0.384 0.005 09 0.009 55 0.0229
PL20-512-β10 512 × 512 × 64 200 10 1.27 0.126 0.382 0.0091 0.0157 0.039
RO20-128 128 × 128 × 128 400 20 1.18 0.0777 0.472 0.004 19 0.0112 0.0198
RO20-264 264 × 264 × 128 200 20 1.20 0.08 0.489 0.004 86 0.0124 0.0212
RO20-264-FL2 264 × 264 × 128 200 20 1.19 0.081 0.475 0.005 39 0.0118 0.0217
RO20-512 512 × 512 × 128 90 20 1.18 0.0754 0.464 0.004 88 0.0116 0.0213
3.1 Fiducial run: short- and long-term evolution
We first explore the gravitoturbulent flow on long time-scales. To
that end, we ran a simulation, labelled PL20-128, with horizontal
resolution similar to that of Shi & Chiang (2014) (≈6.5 cells per H0)
in a box of size Lx = Ly = 20 for almost 3000 −1. This simulation
will be considered as our reference test run. Some average turbulent
properties are summarized in Table 1 and their evolution in time
shown in Fig. 1 (blue/dark curves).
A quasi-steady turbulent state as well as a thermodynamic equi-
librium is obtained within a few tens of orbits, characterized by
Q  1.24, very similar to the value calculated by Shi & Chiang
(2014). We verified that the average transport efficiency α, defined
in equation (14), matches the prediction of Gammie (2001),
α = 1
q2τcγ (γ − 1) , (18)
which is based on total energy conservation. The average kinetic
energy associated with the fluctuations remains smaller than the
average internal energy, as in classical 2D razor thin disc simu-
lations. The gravitational stress is always positive and contributes
to most of the angular momentum transport, while the Reynolds
stress Hxy is subdominant, in agreement with Shi & Chiang (2014).
We also checked that these results do not depend on initial condi-
tions. For that purpose, we ran a simulation PL20-128-b, initialized
with a different noise realization (i.e. different Fourier amplitudes).
This simulation exhibits similar Q and mean average properties as
PL20-128.
While all mean quantities fluctuate on a time-scale of several or-
bits, Fig. 1 shows that Hxy and kinetic energy Ec undergo significant
variations on a much longer time-scale, of the order of ∼1000−1.
They exhibit bursts of activity, for example between t = 500 and
1700 −1, followed by quiescent phases in which the kinetic en-
ergy can be three times smaller (such as between t = 1700 and 2200
−1). The last panel of Fig. 1 shows that bursts of activity are gen-
erally associated with the formation of transient clumps where the
local density can exceed 100 times the background density. These
transient clumps do not collapse during the simulation but could be
important in the process of stochastic fragmentation (Paardekooper
2012). We conclude that this long-time behaviour is not driven by
the correlation length of the turbulence reaching the box size, and
is hence physical and not an artefact of the shearing-box model.
However, it does make it difficult to obtain meaningful averages for
Ec and Hxy because simulations must be run for extremely long.
Although the long-term evolution of the Reynolds stress Hxy is
stochastic and bursty, on short time it oscillates quasi-periodically
between negative and positive values with a frequency close to the
orbital frequency (on average it has a positive contribution). These
fast and regular oscillations are also visible in the simulations of Shi
& Chiang (2014). In Sections 4 and 5, we inspect these oscillations
in more detail and show that they control several aspects of the
dynamics.
Finally, we analysed the time-averaged rms fluctuations as a func-
tion of altitude z. Fig 2 shows that for our test reference simulation
(right-hand panel), the radial turbulent velocity is always larger
than the other components, and in fact is roughly twice the az-
imuthal component, a signature of large-scale density waves. As
z increases, the ratio between vertical and radial rms fluctuations
increases, and approaches 1 as z ∼ H0. This trend also appears in
Shi & Chiang (2014). The non-negligible value of vz in the atmo-
sphere is possibly due to the ‘vertical breathing’ motion of spiral
waves, characteristic of f modes in polytropic or self-gravitating
isothermal gas (Korycansky & Pringle 1995; Ogilvie 1998,
Appendix B). The vertical velocity may also be enhanced by ‘ver-
tical splashing’ as spiral density waves collide. We note that the
average sound speed (or temperature) increases slowly with z and is
always greater than each velocity component taken individually, and
thus the mean motions are slightly subsonic. We discuss in the next
subsection how these different results depend on grid resolution and
box size.
3.2 Dependence on resolution
For a fixed box size Lx = Ly = 20 H0, we compare four differ-
ent setups, detailed in Table 1, with 64, 128, 256, 512 points in
the horizontal directions (ranging from 3.2 points to 25.6 points
per H0.).
Table 1 shows that some quantities like the gravitational stress
and Toomre parameter Q do not vary significantly as resolution is
increased, provided that NX ≥ 128. In contrast, the kinetic energy
and, in particular, the Reynolds stress differ from one resolution to
another, because the higher resolution simulations have not been
run sufficiently long. Recall that in the previous subsection these
quantities were shown to vary on periods of thousands −1, but such
run times are inaccessible at high resolution (the highest resolved
simulations have been run for ∼200−1). The values listed for
the kinetic energy and Reynolds stress are probably statistically
insignificant.
MNRAS 471, 317–336 (2017)
322 A. Riols, H. Latter and S.-J. Paardekooper
Figure 1. Time evolution of various quantities, averaged over a box whose size is Lx = 20, Ly = 20 and Lz = 3 H0. From top to bottom, density-weighted
average Toomre parameter Q, box-averaged kinetic energy, Reynolds stress, gravitational stress and density maximum. The resolution is 128 × 128 × 64 (run
PL20-128).
Fig. 2 shows the vertical profile of rms velocity fluctuations for
PL20-128 and PL20-512 (averaged over 40−1). Increasing resolu-
tion does not seem to affect the balance between each spatial com-
ponent, although the radial and azimuthal fluctuations are 1.5 times
stronger with a resolution of 512 points in Lx and Ly. We note that
vertical motions remain important at z ≥ H0 whatever the resolution
used (even for the lowest one, NX = 64).
Finally, we visually examined the density snapshots at differ-
ent altitudes and resolutions. Fig. 3 shows the density structures
in the mid-plane and one scaleheight above. At z = 0, the turbu-
lence is comprised of large-scale and distinct spiral waves, of radial
length-scale ≈4–5H0, which break non-linearly. However, at larger
z  H0 and in higher resolution runs, the coherence of the spi-
ral waves is disturbed by a form of small-scale turbulence. This
does not appear in our low-resolution runs PL20-128, PL40-256
or PL80-512, nor in previous simulations of 3D gravitoturbulence.
Fig. 4 shows, in particular, a comparison between different resolu-
tions at z = H0 (and for a box of size Lx = 40 H0). Visually these
small-scale motions take the form of wispy deformations of the spi-
ral wavefronts. The fluctuations are relatively strong (of the order
of the host density wave), highly non-axisymmetric, and mainly
located between z = 0.5 H0 and H0. The density at that location is
still relatively high, varying between 0.1 and 1. These small-scale
features remain visible throughout the simulation and are observed
independently of the code, boundary conditions (see Section 3.5)
and box size (see Figs 3 and 4 for comparison).
In conclusion, there is evidence that we achieve convergence for
certain mean gravitoturbulent quantities (Q and Gxy notably) with
resolution, provided that we use more than five grid cells per H0 in x
and y. We cannot decide on the convergence of the kinetic energy Ec
and Reynolds stress Hxy, because they exhibit strong fluctuations on
times of several hundreds of orbits, of the order of or longer than our
MNRAS 471, 317–336 (2017)
GI and turbulence in astrophysical discs 323
Figure 2. Vertical profiles of turbulent rms velocities, time averaged over
40−1 for the run PL512 at high resolution (left) and for the run PL128 at
low resolution (right).
highest resolution simulations. Finally, at resolution larger than 20
cells per H0, we discover a form of small-scale turbulence attacking
the large-scale spiral waves. The behaviour of this new dynamical
feature with resolution is probably still unconverged. And though its
influence on the mean turbulent quantities is probably mild (possibly
appearing in the vertical profiles of ux and uy), it is readily identified
visually in snapshots of the density field, and in later section we
show how it can be quantified through various diagnostics.
3.3 Dependence on box size
Another important test, especially for local models, is the box-size
dependence of the turbulent properties. Tables 1, 2 and 3 show,
respectively, a set of simulations computed for Lx = 20, 40 and 80
H0. To aid comparison, fiducial runs at different box sizes but the
same resolution per H0 are highlighted in red.
First, doubling the horizontal size from 20 to 40 H0 while keeping
the resolution per scaleheight fixed modifies the average Q only
moderately. It is almost 15 per cent larger in PL40-256 compared
with PL20-128. This result seems to be robust since it is obtained
at different resolutions and with different codes (see Section 3.5).
The internal energy also increases. But more strikingly, the ratio
between the gravitational and the Reynolds stress is multiplied by a
factor 2 or even 3. Unlike the variations of Hxy with resolution, this
result is statistically robust, since quantities have been recorded for
at least 900 −1 (run RO40-256) and the same trend is obtained for
all simulations with Lx = 40 H0. Oscillations of Hxy at a frequency of
 are still clearly discernible, although their amplitudes are weaker
and rarely reach negative values.
Finally, for a box even larger, of size Lx = Ly = 80 H0, we show
that turbulence properties are comparable to those obtained for Lx
= Ly = 40 H0. In particular, it appears that Q has settled down to
a value 1.35–1.4. Discrepancies in kinetic energy and Reynolds
stress are probably due to insufficient statistics. In conclusion, we
have some evidence that there exists some critical Lx = Ly between
20 and 40 H0 above which simulations are converged with respect
to the box size. However, further long-term simulations are needed
to nail this down.
Figure 3. Density snapshots of the 3D gravitoturbulent state obtained in
a box of Lx = Ly = 20H0 and resolution of ≈25 cells per H0 with τ c =
20−1. The top panel is a view of the disc mid-plane z = 0, the central
panel is for z = H0 and the bottom panel is a vertical cut, with ρ represented
in a logarithmic scale.
3.4 Dependence on cooling time
We ran a small number of simulations with different cooling times,
not enough to ascertain a critical rate at which fragmentation occurs,
but to fix some basic ideas. A high-resolution simulation, PL20-512-
β10, using τ c = 10−1 and run for 200−1 shows that Q does
not vary significantly compared to when τ c = 20−1 (with same
MNRAS 471, 317–336 (2017)
324 A. Riols, H. Latter and S.-J. Paardekooper
Figure 4. Density snapshots of the 3D gravitoturbulent state obtained in
a box of Lx = Ly = 40 H0 at altitude z = H0. The top panel comes from
a low-resolution run of ≈6.5 cells per H0 (PL40-256), whereas the bottom
panel comes from a high-resolution run, ≈25 cells per H0 (PL40-1024).
resolution). However, the average stress is multiplied by a factor
2, as expected from equation (18). Note that, despite the shorter
cooling time, this simulation failed to exhibit collapsing or even
transient fragments. We did check however that when τ c = 3−1
the system fragments in less than 20 orbits, which supports the
results of Shi & Chiang (2014).
3.5 Code comparison and dependence on numerical details
To verify that our results are robust and that the simulated gravitotur-
bulence does not depend strongly on our numerical implementation,
we compare our results with a different code, RODEO. It employs a
different Riemann solver, Poisson solver and boundary conditions
(see Sections 2.3.1 and 2.3.2).
Tables 1, 2 and 3 show relatively good agreement between the two
codes (within the constraints of run time), in particular for the value
of mean Q. There are nevertheless some differences: the internal
energy is systematically larger in RODEO than in PLUTO. Moreover,
in the Lx = 80 H0 runs, Gxy is also larger. These difference remain
relatively small and could be due to insufficiently long integrations.
We also checked that simulations run with RODEO depend only
mildly on the flux limiter used (see comparison in Table 1 between
RO20-264 and RO20-264-FL2). Finally, we show that results are
independent of the order of the Runge–Kutta method in the time
integrations (see runs PL20-256 and PL20-256-RK3).
4 SM A L L - S C A L E DY NA M I C S A N D
AXI SYMMETRI C OSCI LLATI ONS
We showed previously that gravitoturbulent states are characterized
by regular oscillations of the Reynolds stress and, in the best re-
solved simulations, vigorous small-scale dynamics around z  H0.
In this section, we better understand these structures by extracting
their signatures in Fourier space.
4.1 Turbulent power spectra and small-scale structures
We denote by ûw(kx, ky, z, t) the horizontal 2D Fourier transform
of the weighted turbulent velocity field ρ1/3u (see Section 2.5.2
for definitions). The density-weighted, time-averaged, 2D kinetic
energy power spectrum may then be defined through
E(kx, ky, z) = 
1/3H
−1/3
0
2T
∫ t0+T
t0
∣∣ûw(kx, ky, z, t)∣∣2 dt, (19)
where the time average begins at t = t0 and is carried out over an
interval of duration T.
Table 2. Simulations in larger boxes, Lx = Ly = 40H0.
Run Resolution Time (in −1) τ c Q Ec U Hxy Gxy α
PL40-256 256 × 256 × 64 100 20 1.38 0.0753 0.4 0.002 49 0.0114 0.021 15
PL40-512 512 × 512 × 64 70 20 1.40 0.132 0.524 0.0034 0.0186 0.025
PL40-1024 1024 × 1024 × 64 50 20 1.36 0.076 0.4 0.0028 0.012 0.023
RO40-256 256 × 256 × 128 900 20 1.36 0.117 0.625 0.002 81 0.0213 0.0233
Table 3. Simulations in large boxes, Lx = Ly = 80H0.
Run Resolution Time (in −1) τ c Q Ec U Hxy Gxy α
PL80-512 512 × 512 × 64 100 20 1.41 0.131 0.523 0.0032 0.0182 0.0244
RO80-512 512 × 512 × 128 200 20 1.39 0.149 0.664 1e−6 0.0272 0.0246
MNRAS 471, 317–336 (2017)
GI and turbulence in astrophysical discs 325
Figure 5. 2D horizontal kinetic energy power spectrum E(kx, ky, z) as a
function of Eulerian wavenumbers in simulation PL20-512. The spectrum
is averaged over 40−1. The top panel is for z = 0 and the bottom panel is
for z = H0.
The weighted spectrum, with the one-third scaling in density,
is often used in the study of supersonic and compressible turbu-
lence (Lighthill 1955; Kritsuk et al. 2007). Using simple scaling
arguments, one can show that E(k) ∼ k−5/3 and hence follows the
incompressible Kolmogorov power law. Note, however, that this
scaling fails for highly compressible forcing and cannot be consid-
ered universal (Federrath 2013).
We define next the time- and kx-averaged 1D spectrum
E1D(ky, z) = 1
kx0
∫
E(kx, ky, z) dkx, (20)
and the time-averaged isotropic 1D spectrum
Eiso(k, z) =
∫ ∫
E(kx, ky, z) δ(|k′| − k) dkxdky, (21)
where |k′| =
√
(k′x)2 + (k′y)2 is the radius in horizontal wavenum-
ber space. The kx-averaged spectrum is particularly useful in distin-
guishing the GI wakes (which possess a low ky) from small-scale
turbulence (for which ky is bigger). Note that both features may
possess comparable kx.
Fig. 5 shows the time-averaged 2D power spectrum of gravi-
toturbulence E(kx, ky, z) as a function of the horizontal Eulerian
wavenumbers at two different altitudes z = 0 and H0. Both spectra
are computed from high-resolution simulation data (PL20-512) and
averaged over 40−1. For z = 0, the energy is contained in an
inclined elliptical band along the ky = 0 axis. Turbulent structures
are weakly non-axisymmetric and elongated along the azimuthal
direction with a small pitch angle. In contrast, at z = H0, the signal
appears far more isotropic; energy is clearly spread on to small-
scale non-axisymmetric modes. This broadening of the spectrum
Figure 6. Left: 1D kinetic energy spectrum E1D. Right: 1D isotropic spec-
trum Eiso. Two different altitudes have been considered: blue curves are for z
= 0 and green curves are for z = H0. Data are computed from the simulation
PL20-512 and averaged in time over 40−1.
is undoubtedly related to the small-scale disturbances afflicting the
spiral waves described earlier in Section 3.2.
To better distinguish the small-scale motions, we plotted E1D(ky,
z), the kx-averaged 1D spectrum, in the left-hand panel of Fig. 6
evaluated at z = 0 and H0. Note that the ratio of large-scale energy
(ky ≤ 2π/Ly) to small-scale energy (ky = 50 for instance) strongly
varies with altitude z. At z = H0, this ratio is 15 times smaller than at
z = 0, indicating a significant distribution of energy to smaller scales
higher up in the disc. The spectrum in ky appears also less steep for
azimuthal scales larger than 0.5H0, which suggests that the energy
of non-axisymmetric modes is injected and cascades differently at
those two altitudes. Lower resolution runs do not exhibit nearly the
same amount of power on the resolved small scales at z = H0.
Lastly, out of interest we plot the isotropic spectrum Eiso(k, z),
even though it is difficult to distinguish the two features in it because
both share similar kx. Note, however, that the slope of E(k) follows
a Kolmogorov law for k < 10H−10 but then appears much steeper
(well before approaching the grid scale).
In conclusion, the small-scale parasitic turbulence attacking at z
= H0 is quantifiable and possesses a clear signature in the Fourier
decomposition of the flow. This supports the conclusions of the pre-
vious section obtained by visual inspection of the density maps. In
light of Figs 3, 4 and the present results, we assume that the mech-
anism producing those small-scale features could be more efficient
than a simple inertial turbulent cascade and could be potentially due
to a secondary instability of the large-scale structure.
4.2 Large-scale modes: spiral waves and axisymmetric
epicyclic oscillations
The turbulent power spectra shown in Fig. 6 reveal that kinetic
energy is not concentrated at wavelengths ≈4–5H0, the typical size
of GI density wakes. Apart from the short scales, energy keeps
increasing to very large scales ≈Lx. In this section, we investigate
in more detail the nature and origin of these large-scale features.
Fig. 7 shows the time evolution of kinetic energy associated with
large-scale m = 0 axisymmetric modes (top panel) and m = 1
non-axisymmetric shearing waves (bottom two panels), in the mid-
plane region z = 0 and for the run PL20-512. Surprisingly, we found
that the largest scale axisymmetric mode (kx ≡ kx0 = 2π/Lx, ky =
0) dominates the energy content. This mode oscillates at a very
well defined angular frequency ωF = 0.83 κ , close to the orbital
or epicyclic frequency. This is indicative of a large-scale inertial
oscillation modified by self-gravity. The mode is clearly responsible
MNRAS 471, 317–336 (2017)
326 A. Riols, H. Latter and S.-J. Paardekooper
Figure 7. Top: time evolution of the kinetic energy E(kx0 , 0, 0),
E(2kx0 , 0, 0) and E(3kx0 , 0, 0) associated with the first few axisymmet-
ric modes. Bottom: kinetic energy of several non-axisymmetric shearing
waves, with ky = 2π/Ly , labelled by their Lagrangian wavenumber  (to
improve readability, we separate the even  in the centre panel from the odd
 in the bottom panel). All the spectral quantities have been computed in the
mid-plane z = 0.
for the pseudo-periodic variation of the Reynolds stress observed in
Section 3. Smaller scale axisymmetric modes such as kx = 2kx0 or
kx = 3kx0 are negligible compared to the fundamental mode (their
specific energy decreases with their radial scale) and their evolution
appears more chaotic.
Fig. 7 (bottom panels) shows that, energetically, individual spi-
ral or non-axisymmetric waves grow and decay aperiodically, but
their amplitude remains on average weaker than the largest scale
axisymmetric mode (in particular  = 2 at t = 42/ or the orange
and blue ones at t = 57/ and 60/). The energy of the strongly
amplified waves evolves quasi-exponentially during their growth
phase before soon decaying. Overall, the weaker strength of spiral
waves is surprising, given that we expect gravitoturbulence to be
primarily supported by non-axisymmetric structures (axisymmetric
modes being stable for Qturb ≥ 1). It is important to note, however,
that we only have shown the waves’ kinetic energy. The associated
density perturbations associated with spiral waves in fact can be
larger than for axisymmetric modes. Indeed, we find that the den-
sity of the fundamental axisymmetric mode is rather small and does
not oscillate regularly in time.
Are these axisymmetric modes physical or artificially enhanced
by the box size, or even the shearing periodic boundary conditions?
To answer this question, we determined the energetics of the ax-
isymmetric structures in larger boxes. Fig. 8 shows that when Lx
= 40 and 80 H0, the kinetic energy of each individual axisymmet-
ric mode is smaller than when Lx = 20 H0 but in total represents
still ∼20 per cent of the energy. Also, for Lx ≥ 40 H0, the kinetic
Figure 8. Ratio between the kinetic energy of individual axisymmetric
modes and the total kinetic energy (averaged over at least 100−1) as a
function of horizontal box size Lx = Ly (but same resolution per H0). Black,
purple, blue, green and red account for modes kx = 2π/80j , with j = 1,
. . . , 5 in this order. The diamond marker indicates the fundamental mode
with wavelength equal to the box size.
energy associated with each axisymmetric mode appears to con-
verge to a similar value. Most importantly, for a sufficiently large
box, the fundamental mode (on the box size) ceases to be domi-
nant, its place taken by a shorter scale mode. Larger boxes can host
multiple axisymmetric inertial oscillations with frequencies close
to κ . These results indicate that the large-scale axisymmetric dy-
namics is not controlled by the box size, which is consistent with
Section 3.2 that showed some turbulent quantities had converged
by Lx = Ly  40 H0.
Finally, we checked that for large boxes, axisymmetric modes
still exhibit an oscillatory behaviour with a frequency close to 0.8κ
(this is in particular true for the harmonics kx = 2kx0 or kx = 3kx0 ).
We checked also that the dynamics and amplitude of these modes
are similar in both PLUTO and RODEO simulations. In sum, oscillatory
axisymmetric waves are an important component of the GI dynam-
ics, regardless of the box size. Although their amplitude might be
enhanced artificially if the box size is too small, their existence is
likely to have a physical origin.
4.3 Origin and role of axisymmetric modes
Though non-axisymmetric shearing waves form the backbone of the
non-linear GI dynamics, we have seen that large-scale axisymmetric
oscillations are conspicuous participants. Several questions may
then be asked.
(a) What is their origin? Though linearly stable, how can they
reach and sustain large amplitudes?
(b) What role do they play in the dynamics? Are they involved
in the subcritical transition to, and the maintenance of, turbulence,
or are they just passive modes?
In order to answer (a), we determine the force balance associated
with the fundamental axisymmetric mode kx = 2π/Lx in PL20-
512. It can be decomposed into a high- and low-amplitude standing
wave, uˆ1(t) cos(kxx) and uˆ2(t) sin(kxx), respectively, with uˆ1 	 uˆ2.
Fig. 9 shows that the first wave possesses an inertial character since
it is mainly driven by the Coriolis force and partially restored by
pressure. Self-gravity is always opposed to the latter and tends to
destabilize the wave (although Q is too large for the wave to be-
come unstable). Non-linear terms seem negligible in its dynamical
evolution, although they have a cumulative positive feedback on its
energy budget. We checked that the force balance is similar in the
azimuthal direction (not shown here) and for the smaller amplitude
second wave. In sum, we identify a coherent quasi-linear epicyclic
MNRAS 471, 317–336 (2017)
GI and turbulence in astrophysical discs 327
Figure 9. Radial force balance of the large-scale axisymmetric mode kx0 =
2π/Lx in the mid-plane region z = 0.
oscillation. We found that this oscillation is maintained for the
entirety of our simulations (which is 3000−1 for PL20-128b), and
the mode cannot be a vestige of our initial condition, since turbulent
viscosity would have destroyed it in a shorter time-scale (typically
tν  1/(αcsk2)  500−1). This suggests that large-scale ax-
isymmetric perturbations are non-linearly excited and regenerated
by gravitoturbulence.
We performed several additional simulations, shown in
Appendix C, that firmly indicate that the mode kx = kx0 is the result
of a non-linear baroclinic interaction between non-axisymmetric
density waves. This dynamics is reminiscent of the mode coupling
proposed by Laughlin, Korchagin & Adams (1997) for gravito-
turbulent discs and by Lee (2016) in the context of planet–disc
interactions. Laughlin et al. (1997) showed, in particular, that self-
interaction of m = 2 modes can lead to a non-linear growth of the
m = 0 mode in global discs. However, their initial density profile is
subject to a Rossby wave instability, which complicates interpreta-
tion of their results.
Regarding question (b), there are several reasons to think that
these modes are active participants in the disc dynamics. First,
they appear to be sufficiently coherent in time to enter in reso-
nance with a pair of inertial waves. Such a resonance, based on
a non-linear triadic interaction, is known to provide a path to-
wards parametric instabilities. We investigate this possibility in the
next section.
Secondly, the large-scale axisymmetric modes could be involved
in the regeneration and instability of leading shearing waves and
thus in the sustenance of the whole GI dynamics. A potentially
related regeneration process has been formulated for the magne-
torotational dynamo problem (Herault et al. 2011). Individual non-
axisymmetric waves are sheared out into strongly trailing decaying
structure. Long-lived dynamics can only be excited if new leading
waves are seeded each t = 23Ly/Lx−1. Herault et al. (2011)
found that axisymmetric modes play the role of ‘walls’ and thus
‘confine’ the non-axisymmetric waves. The dying trailing struc-
tures can then be reflected and give birth to a new leading shearing
wave. This mechanism is probably essential to any kind of sustained
non-axisymmetric turbulence in the presence of a background shear
flow. Analogously, axisymmetric zonal flows are known to support
a linear instability of self-gravitating spiral waves (Lithwick 2007;
Vanon & Ogilvie 2016). We suspect then that such an instability is
at work in our simulations of gravitoturbulence and might feature
prominently in the self-sustaining process. The dynamics involving
the non-linear regeneration of axisymmetric modes (see previous
paragraph and Appendix C) as well as the non-axisymmetric GI lin-
ear instability (supported by axisymmetric structures) might form
the basis of the fully developed turbulence.
5 E X C I TAT I O N O F A PA R A M E T R I C
INSTABILITY
Section 4 revealed that a large-scale axisymmetric epicyclic mode,
oscillating at a frequency ωF  0.8 κ , naturally emerges from grav-
itoturbulence. From the perspective of the smaller scales, this in-
ternal coherent oscillation may function equivalently to an external
periodic forcing or deformation of the disc. Rotating flows subject
to periodic forcing or deformations are known to excite parametric
instability (Goodman 1993), and in fact the standing epicycle we
witness in our simulations is similar in many respects to that treated
by Fromang & Papaloizou (2007) who showed that an axisymmet-
ric density wave is subject to such an instability. In this section, we
explore this idea and consider whether a similar instability occurs
in our problem and is responsible for the small-scale turbulence
appearing at z  H0.
5.1 Theoretical expectations
5.1.1 Onset of instability
Parametric instabilities can arise when an oscillator undergoes a
forcing at twice its natural frequency. In accretion discs, they occur
when a large-scale, time-periodic disturbance enters into resonance
with a pair of small-scale inertial waves (Goodman 1993). This
fundamentally 3D instability has been studied in the context of
eccentric (Ogilvie 2001; Papaloizou 2005; Barker & Ogilvie 2014)
and warped discs (Gammie, Goodman & Ogilvie 2000; Ogilvie &
Latter 2013), and in the case of mean-motion resonances between
a disc and its companion (Lubow 1991). More recently, it has been
shown to attack spiral density waves excited in the disc (Bae et al.
2016a,b).
Although relatively robust, parametric instabilities require cer-
tain conditions to work. If we denote by ωi1 and ωi2 the frequencies
of two inertial waves and ωF the frequency of the large-scale oscil-
lation, it is possible to show that growth is obtained when (Fromang
& Papaloizou 2007)
ωF = ωi1 + ωi2 . (22)
To obtain a constraint on the background oscillation ωF, we need
the dispersion relation for each inertial wave. This is available in
the WKBJ limit of short radial and vertical wavelengths (Goodman
1993). We consider only the coupling of ‘neighbouring’ inertial
wave branches, for which ωi1 ≈ ωi2 , because of the density of the
spectrum on short scales. The resonance condition immediately
yields
ωi1 = ωi2 =
1
2
ωF, (23)
and because ω2i1 < κ
2
, we deduce that ωF < 2κ . From the dispersion
relation again, we find that resonance implies ω2i1 > N
2
, and so the
condition for the excitation of parametric instability is
2N < ωF < 2 κ (24)
(Goodman 1993; Fromang & Papaloizou 2007). As discussed by
Bae et al. (2016a), the influence of buoyancy is considerable and
must stabilize the flow above a certain disc altitude zcrit for which
N > κ .
MNRAS 471, 317–336 (2017)
328 A. Riols, H. Latter and S.-J. Paardekooper
We should emphasize that these simple theoretical arguments
have their limits. First, condition (24) is derived for axisym-
metric disturbances, though for weakly non-axisymmetric waves
(ky  k) Goodman (1993) showed that condition (24) still holds.
Unlike axisymmetric waves, however, their lifetime is finite and they
can only be amplified transiently. Secondly, condition (24) only ap-
plies to localized wave packets, and not to vertically global modes.
In principle, one could solve the full linear eigenvalue problem with
vertical structure, similarly to Ogilvie & Latter (2013) and Barker
& Ogilvie (2014), but we leave this task for another time. In any
case, given the disordered background in which the inertial waves
find themselves, they may manifest more as localized packets rather
than larger scale organized structures. Finally, self-gravity has been
neglected, though for Q  1 and kx large, the dispersion relation of
linear inertial disturbance is almost unchanged when compared to
the case without self-gravity (see fig. 3 in Mamatsashvili & Rice
2010).
In the isothermal case, Fromang & Papaloizou (2007) estimated
the growth rate of parametric instability due to an axisymmetric den-
sity wave to be roughly proportional to the product of the fractional
amplitude of the density wave and ωF. Applying this formula to
the large-scale epicyclic oscillations in our simulations, this gives a
growth rate ∼0.1–0.2, and so we expect the parametric instability
to develop after few orbits in our simulations, unless it is impeded
by other participants in the gravitoturbulence.
5.1.2 Propagation and breaking of inertial waves
Suppose that the conditions for exciting inertial waves are met. How
do these waves propagate and dissipate in the disc? First note that,
as a localized wave packet travels upwards from the mid-plane, the
background state around it changes and, as a consequence, the wave
packet properties will also change. In fact, the wave will ‘refract’
and become ‘channelled’ by the disc’s vertical structure, which acts
similarly to a waveguide (Lubow & Pringle 1993; Bate et al. 2002).
The packet’s upward trajectory will bend more and more radially,
before ultimately settling into a normal eigenmode of the vertically
global problem, travelling solely in the radial direction.
Perhaps before that occurs, the wave will break non-linearly. In
an unstratified medium, inertial waves are known to degenerate into
turbulence via the action of secondary instabilities, driven by shear,
which occurs when displacements are comparable to wavelengths.
If the inertial wave has a form uˆ exp i(kξ ξ − ωt), where ξ is in the
direction of propagation, then breaking occurs when the particle
velocities exceed the phase velocity
uˆ ≈ ω/kξ . (25)
This relation is equivalent to saying the Rayleigh stability criterion
is locally violated (the angular momentum gradient becomes locally
negative; see Wienker 2016). This secondary instability is actually
then another parametric instability involving a primary inertial wave
and two daughter waves. The route to turbulence in such systems
may then be viewed as a cascade of parametric instabilities.
In a stratified disc, inertial waves naturally break as they prop-
agate towards the surface of the disc. Indeed, since the density
decreases above the mid-plane, the velocity components increase
with z (consider the linear eigenfunctions of a stratified disc, Ap-
pendix A). We expect inertial waves to break above a critical height
that depends on the mode structure. For instance, an n = 1 iner-
tial mode, for which uˆx ∝ z (see Appendix A), possesses a critical
height around z  H0. Higher order n modes break at altitudes some-
what lower, but all are consistent with our observation of small-scale
Figure 10. Power spectral density (PSD) of the averaged kinetic energy
as a function of temporal frequencies. The upper panel was obtained from
simulations PL20-512 and the lower panel from PL20-128. The red arrows
on the right indicate the frequency of the axisymmetric oscillation (note that
the kinetic energy oscillates at twice the natural frequency of the mode).
The left arrows show a secondary peak at 0.8 κ , which might correspond to
the resonant inertial modes.
activity above the mid-plane. Most of the kinetic energy transfers to
small scale at this critical height or above. Energy is then expected
to saturate at a value of the order of u2 = ω2/kξ , although a more
complete calculation, in the context of eccentric discs, shows that
it also increases with the amplitude of the parent wave (Wienker
2016).
5.2 Evidence of parametric instability from the numerical
data
5.2.1 Temporal spectra
In Section 4, we identified a large-scale axisymmetric mode that
oscillates at a frequency ωF ≈ 0.8κ . According to condition (24),
a parametric instability is then possible in the vicinity of the mid-
plane where N = 0 and should induce the growth of a pair of
inertial waves at a frequency ω = ωF/2. To assess whether or not
the instability occurs in our simulations, we compute the temporal
Fourier spectrum of the turbulent signal. Fig. 10 (top) shows the
power spectrum of the box-averaged kinetic energy as a function
of (temporal) frequencies for the PL20-512 run. A strong peak is
observed at a frequency 2ωF, which corresponds to the large-scale
axisymmetric mode (quadratic quantities like energy oscillate at
twice the frequency of the associated mode). From the right to
the left, we see a second peak, smaller in amplitude, at half the
frequency of the fundamental peak ω = ωF. If we look further to
the left, several distinct peaks are seen respectively at ω =ωF/2 and
ω = ωF/4. This cascade of subharmonics peaks at frequency ωF/2n
may be interpreted as a direct signature of parametric instability via
inertial mode couplings.
A similar result is obtained from the PSD of PL20-128, calculated
over a much longer time 3000−1, although the signal is more noisy
and the peaks less pronounced. Because of the rich assortment of
long-time variability (discussed in Section 3.1), it is difficult to
extract as clear a signature of the parasitic inertial waves, as in the
shorter time high-resolution run PL20-512.
In larger boxes with Lx = 40−1 and 80−1, the signal is spread
over more frequencies and the spectrum is more difficult to inter-
pret. Now multiple large-scale axisymmetric modes are hosted in the
box with similar frequencies (cf. Fig. 8), each potentially unstable
MNRAS 471, 317–336 (2017)
GI and turbulence in astrophysical discs 329
Figure 11. Helicity in the poloidal (x, z) plane weighted by the average density profile in z, for different time. All arrows indicate a direction of 45◦ with
respect to the vertical axis. The curly bracket shows some examples of intricate and tangled patterns probably associated with inertial waves. In any of these
figures, it is possible to see an alternation of inclined red and blue bands, with typical wavelength ≈0.1–0.2H0.
to different parametric resonances. The ensuing turbulent dynam-
ics is more complicated and the temporal spectrum a less useful
diagnostic in this case.
5.2.2 Helicity maps
Another way to trace the inertial waves and the associated paramet-
ric instability is to compute the helicity of the flow
H = u · (∇ × u). (26)
Because inertial waves satisfy ∇ × u  ±ku, where k is their
wavenumber, they are intrinsically helical, and the direction of
propagation of a given packet corresponds to a surface of constant
helicity. These different properties have been studied in particular
in the context of the solar and geodynamo (Singer & Olson 1984;
Davidson 2014; Ranjan & Davidson 2014; Wei 2015).
Fig. 11 shows the helicity H, weighted by the averaged vertical
density profile, projected on to the (x, z) plane at different y and
times. Helical structures have a typical size of  H, in both the z
and x directions. They form inclined and interleaved layers of posi-
tive and negative helicity, corresponding to the troughs and crests of
inertial waves. Thus, the associated wavevectors cut across the lay-
ers, while the group velocities point along the layers (Bordes et al.
2012). In Fig. 11, the inclined black arrows have been superimposed
to highlight the inertial layers.
The preferential angle of propagation is in a range between ±30◦
and ±50◦ measured from the vertical axis. However, using the local,
linear and axisymmetric dispersion for inertial waves (Goodman
1993), and setting ωi = 12F ≈ 0.42κ , we find that unstable inertial
waves should possess a group velocity pointing within a range of
angles30◦ with respect to the rotation axis. The slight discrepancy
with the simulation data we attribute to wave channelling by the disc
structure, as discussed earlier. In fact, the bottom-right panel of Fig.
11 presents a nice example of wave channelling.
5.2.3 Helmholtz decomposition and connection with the
small-scale turbulence
We showed in Sections 3.2 and 4.1 vigorous small-scale turbulent
activity taking place z ∼ H0, while in the previous subsection we
compiled numerical evidence for inertial waves propagating up-
wards at altitudes less than H0. In this subsection, we attempt to
connect the two.
First, we check that the small-scale motions are not directly
excited by large-k gravitational modes. We performed a simula-
tion, starting from a gravitoturbulent state, taken from PL20-512,
and filtered the gravitational potential on small scales by retaining
only modes with kx < 6kx0 and ky < 6ky0 . Despite the filtering, the
small-scale turbulence persisted for long times in the reconfigured
simulation.
This activity, however, might still consist of a field of small-scale
p modes, different in character to inertial waves. In fact, short-
wavelength linear p modes preferentially localize at the disc sur-
face (Korycansky & Pringle 1995; Ogilvie 1998), where the turbu-
lence is observed. They may be generated by disordered motions on
large scales, in particular by the non-linearly breaking of large-scale
wakes, which may favour the less dense upper layers.
One way to disentangle these different small-scale motions is to
decompose the velocity field into compressible and incompressible
parts. Inertial waves are generally incompressible (or rather anelas-
tic) if they vary on a vertical scale comparable to or less than the
background support (¯kz  1). On the other hand, p modes are asso-
ciated with compressible motions, though they are not necessarily
curl free and can also possess an incompressible part. To separate
the compressible part of a vector field F from its incompressible
part, we use the Helmholtz decomposition
F = Fc + Fic = ∇ϕ + ∇ × , (27)
whereϕ is a scalar field, defined up to a constant, and  a vector field
defined up to a gradient field. The first term is the compressible and
potential part of F. This is a curl-free component, i.e. ∇ × Fc = 0.
MNRAS 471, 317–336 (2017)
330 A. Riols, H. Latter and S.-J. Paardekooper
Figure 12. Velocity components ux (left) and uz (right) in the mid-plane z = 0, weighted by the average density at that location. The top panels represent the
total field, the centre and bottom panels represent respectively its incompressible and compressible parts.
The second term is the incompressible and solenoidal part, and it is
divergence free, i.e. ∇ · Fic = 0. Applying the divergence operator
on both sides of equation (27) and using the solenoidal properties,
we obtain
∇2ϕ = ∇ · F. (28)
Therefore, to find the potential ϕ, we have to solve a Poisson equa-
tion with a source term ∇ · F. In the shearing-box frame, the method
is very similar to the one used to obtain the disc’s potential (see Sec-
tion 2.3.2), so we simply adapted our Poisson solver to deal with a
different source term. Several tests on simple flows were completed
before applying our method to the full turbulent field. The calcula-
tion of  requires solving three Poisson equations but in practice
we obtain it by subtracting the compressible part from the initial
field F.
We applied the Helmholtz decomposition to the vector F = ρ0 u,
where ρ0 is the vertical profile of the horizontally averaged den-
sity. This decomposition permits us to analyse the incompressibil-
ity of the flow, remembering that the incompressible component
(ρ0 u)ic = ∇ ×  will preferentially exhibit an inertial character.
Fig. 12 shows the incompressible (centre panels) and compress-
ible parts (bottom panels) of the flow components ρ0 ux and ρ0 uz in
the mid-plane z = 0. For the radial component ρ0 ux, the two parts
have similar amplitudes and both exhibit large-scale spiral waves.
In the incompressible part, the waves manifest as large-scale fea-
tures, while in the compressible part they are visible as thin shocks.
In addition, we observe small-scale bundles, of typical size 0.2 H0,
in the incompressible part exclusively. The vertical velocity (right-
hand panels) is clearly dominated by the incompressible component,
which exhibits thin filament structures organized along the wakes.
We attribute both the small vertical velocity filaments and the hori-
zontal bundles with incompressible inertial waves developing at the
mid-plane.
Fig. 13 shows the same decomposition at z = H0. Here the small-
scale dynamics is far more developed, in both the radial and vertical
velocities. Most important, however, is that these small-scale fea-
tures are almost completely confined to the incompressible parts
of both velocity components. The compressible parts exhibit pre-
dominantly large-scale structure associated with the gravity wakes,
in the case of uz arising probably from the vertical breathing or
MNRAS 471, 317–336 (2017)
GI and turbulence in astrophysical discs 331
Figure 13. Velocity components ux (left) and uz (right) in the plane z = H0, weighted by the average density at that location. The top panels represent the total
field, the centre and bottom panels represent respectively its incompressible and compressible parts.
splashing of the wakes. There is no evidence of small-scale p modes
as might be generated by wave breaking at this altitude. In summary,
the marked incompressibility of small-scale turbulence at z = H0 is
strong evidence that it is indeed comprised of inertial waves excited
by parametric instability, breaking and mixing at these altitudes (as
described in Section 5.1.2).
6 C O N C L U S I O N S A N D A S T RO P H Y S I C A L
I M P L I C ATI O N S
We performed a set of 3D vertically stratified shearing-box sim-
ulations of gravitoturbulence, reproducing the local behaviour of
marginally gravitationally unstable discs. For a fixed cooling time
τ c = 20−1, we provided evidence showing that averaged turbulent
properties converge for box size larger than 40 H0 and for sufficient
resolution. The convergence of some quantities, however, is difficult
to ascertain, such as the kinetic energy and Reynolds stress, due to
their bursty and very long time variation (hundreds of orbits).
Our highest resolution runs exhibit small-scale disordered mo-
tions at around z ≈ H0. It is unclear whether the properties of
this small-scale turbulence have yet converged. We compile the-
oretical and numerical evidence that the origin of the turbulence
is a parametric instability involving a pair of inertial waves and
a large-scale axisymmetric epicycle emerging naturally from the
gravitoturbulence. It is likely that this large-scale mode is an im-
portant participant in the onset and sustenance of non-linear non-
axisymmetric GI.
We discount alternative causes of the turbulence. Kelvin–
Helmholtz instability arising from the gravity wakes’ horizontal
shearing should lead to activity at all altitudes, not just at z =
H0. Non-linear breaking of the gravity wakes should lead to com-
pressible motions, but we show that the small-scale turbulence is
primarily incompressible at z = H0. Finally, vertical splashing due
to wake collisions should give rise to larger scale motions dom-
inated by the vertical velocities, which is numerically observed
only in the compressible part of the flow. The small-scale tur-
bulence we find is incompressible and not dominated by vertical
motions.
It should be stated clearly that our simulations employ large
horizontal domains (>20H) for which it is challenging to
MNRAS 471, 317–336 (2017)
332 A. Riols, H. Latter and S.-J. Paardekooper
justify the locality of the shearing-box model. This is certainly
the case for the thicker PP discs, where H/r ∼ 0.05, but there
is more leeway with AGN. We persist with this numerical tool
because it affords us a well-defined and relatively clean platform
to probe shorter scale dynamics at high resolution. Previous grid-
based global simulations can compete on vertical resolution, but
suffer from poor azimuthal resolution especially (e.g. Michael et al.
2012; Steiman-Cameron et al. 2013). The radial boundary condi-
tions in global simulations impose further complications that muddy
their interpretation.
There are several astrophysical implications of this work. We
focus on those most relevant to PP discs. Small-scale turbulence
disturbs the coherence of the large-scale spiral waves, while not ex-
actly destroying them. This ‘scrambling’ may be possible to observe
in direct imaging of large-scale spiral arms, giving them a somewhat
‘flocculent’ appearance (see discussion in Bae et al. 2016b).
Perhaps of greatest importance is the impact of small-scale tur-
bulence on marginally coupled particles off the mid-plane. Two-
dimensional planar simulations of GI with dust indicate that this
class of particle can collect in filaments, sometimes achieving over-
densities two orders of magnitude over the background. Moreover,
their relative velocities are diminished to below 1 per cent of the
sound speed, facilitating gravitational collapse (Gibbons, Rice &
Mamatsashvili 2012; Gibbons, Mamatsashvili & Rice 2014, 2015;
Shi et al. 2016). If, however, within these filaments there is
an additional source of turbulent motions, on a similar scale
(as we find), then such overdensities may be difficult to at-
tain, and collapse more difficult. Moreover, because the eddy
sizes of the parasitic turbulence are H, it is more likely that
pre-collisional pairs are decorrelated and will thus have large
impact speeds. However, for these effects to be significant the ver-
tical thickness of the particle subdisk must be ∼H, and this may not
be a common situation, or one that lasts very long, in the coupled
evolution of the disk and dust.
More generally, parasitic turbulence will enhance whatever diffu-
sion is present. Fig. 2 shows that better resolved turbulence exhibits
larger rms velocities by factors of ∼1, and one would expect dif-
fusivities to be similarly enhanced. As mentioned earlier, even our
highest resolved simulations are still probably unconverged with
respect to this physics.
Though not explored in this paper, small-scale motions may make
fragmentation more difficult. It should be emphasized that the typi-
cal inertial wave motions may be too slow to be effective, and being
localized to z ∼ H will certainly be inefficient at the mid-plane where
most fragmentation begins. Work specifically targeting fragmenta-
tion in high-resolution 3D disc models will appear in a separate
paper.
Finally, one may ask how these features affect the magnetohy-
drodynamics in PP and AGN discs, as both may be moderately well
ionized in their outer regions (Menou & Quataert 2001; Turner, Lee
& Sano 2014). The preponderance of helical motions automatically
provides a means of dynamo action, but the magnetorotational in-
stability (in some non-ideal form) must also play a role here. One
may then imagine a rich and complex interplay between the GI,
the parasitic turbulence feeding off it and the magnetorotational
instability.
AC K N OW L E D G E M E N T S
The authors thank the reviewer for a set of helpful comments,
Richard Booth, Charles Gammie, Sebastien Fromang, and Ji-Ming
Shi for important feedback on an earlier draft of the paper, and
Kacper Kornet for advice on numerical issues. This research is
partially funded by STFC grant ST/L000636/1. Many of the simu-
lations were run on the DiRAC Complexity system, operated by the
University of Leicester IT Services, which forms part of the STFC
DiRAC HPC Facility (www.dirac.ac.uk). This equipment is funded
by BIS National E-Infrastructure capital grant ST/K000373/1 and
STFC DiRAC Operations grant ST/K0003259/1. DiRAC is part of
the UK National E-Infrastructure.
R E F E R E N C E S
Bae J., Nelson R. P., Hartmann L., Richard S., 2016a, ApJ, 829, 13
Bae J., Nelson R. P., Hartmann L., 2016b, ApJ, 833, 126
Balbus S. A., Papaloizou J. C. B., 1999, ApJ, 521, 650
Barker A. J., Ogilvie G. I., 2014, MNRAS, 445, 2637
Bate M. R., Ogilvie G. I., Lubow S. H., Pringle J. E., 2002, MNRAS, 332,
575
Bordes G., Moisy F., Dauxois T., Cortet P.-P., 2012, Phys. Fluids, 24, 014105
Boss A. P., 1997, Science, 276, 1836
Boyd J. P., 2001, Chebyshev and Fourier Spectral Methods, 2nd edn. Dover,
New York
Cameron A. G. W., 1978, Moon Planets, 18, 5
Christiaens V., Casassus S., Perez S., van der Plas G., Me´nard F., 2014, ApJ,
785, L12
Davidson P. A., 2014, Geophys. J., 198, 1832
Durisen R. H., Boss A. P., Mayer L., Nelson A. F., Quinn T., Rice W. K. M.,
2007, Protostars and Planets V. Univ. Arizona Press, Tucson, p. 607
Federrath C., 2013, MNRAS, 436, 1245
Fromang S., Papaloizou J., 2007, A&A, 468, 1
Fukagawa M. et al., 2013, PASJ, 65, L14
Gammie C. F., 2001, ApJ, 553, 174
Gammie C. F., Goodman J., Ogilvie G. I., 2000, MNRAS, 318, 1005
Gibbons P. G., Rice W. K. M., Mamatsashvili G. R., 2012, MNRAS, 426,
1444
Gibbons P. G., Mamatsashvili G. R., Rice W. K. M., 2014, MNRAS, 442,
361
Gibbons P. G., Mamatsashvili G. R., Rice W. K. M., 2015, MNRAS, 453,
4232
Goldreich P., Lynden-Bell D., 1965, MNRAS, 130, 125
Golub G., Van Loan C., 1996, Matrix Computations. Johns Hopkins
University Press, Baltimore, MD
Goodman J., 1993, ApJ, 406, 596
Goodman J., 2003, MNRAS, 339, 937
Herault J., Rincon F., Cossu C., Lesur G., Ogilvie G. I., Longaretti P.-Y.,
2011, Phys. Rev. E, 84, 036321
Hirose S., Shi J.-M., 2017, MNRAS, 469, 561
Kalas P. et al., 2008, Science, 322, 1345
Korycansky D. G., Pringle J. E., 1995, MNRAS, 272, 618
Koyama H., Ostriker E. C., 2009, ApJ, 693, 1316
Kritsuk A. G., Padoan P., Wagner R., Norman M. L., 2007, in Shaikh D.,
Zank G. P., eds, AIP Conf. Ser. Vol. 932, Turbulence and Nonlinear
Processes in Astrophysical Plasmas. Am. Inst. Phys., New York, p. 393
Laughlin G., Korchagin V., Adams F. C., 1997, ApJ, 477, 410
Lee W.-K., 2016, ApJ, 832, 166
LeVeque R., 2002, Finite Volume Methods for Hyperbolic Problems.
Cambridge Univ. Press, Cambridge
Lighthill M., 1955, IAU Symp. 2, ch. 22, p. 121
Lithwick Y., 2007, ApJ, 670, 789
Lubow S. H., 1991, ApJ, 381, 259
Lubow S. H., Pringle J. E., 1993, ApJ, 409, 360
Mamatsashvili G. R., Rice W. K. M., 2010, MNRAS, 406, 2050
Menou K., Quataert E., 2001, ApJ, 552, 204
Michael S., Steiman-Cameron T. Y., Durisen R. H., Boley A. C., 2012, ApJ,
746, 98
Mignone A., Bodo G., Massaglia S., Matsakos T., Tesileanu O., Zanni C.,
Ferrari A., 2007, ApJS, 170, 228
Ogilvie G. I., 1998, MNRAS, 297, 291
MNRAS 471, 317–336 (2017)
GI and turbulence in astrophysical discs 333
Ogilvie G. I., 2001, MNRAS, 325, 231
Ogilvie G. I., Latter H. N., 2013, MNRAS, 433, 2420
Paardekooper S.-J., 2012, MNRAS, 421, 3286
Papaloizou J. C. B., 2005, A&A, 432, 743
Ranjan A., Davidson P. A., 2014, J. Fluid Mech., 756, 488
Rice W. K. M., Paardekooper S.-J., Forgan D. H., Armitage P. J., 2014,
MNRAS, 438, 1593
Riols A., Latter H., 2016, MNRAS, 460, 2223
Roe P. L., 1981, J. Comput. Phys., 43, 357
Roe P. L., 1986, Annu. Rev. Fluid Mech., 18, 337
Shi J.-M., Chiang E., 2014, ApJ, 789, 34
Shi J.-M., Zhu Z., Stone J. M., Chiang E., 2016, MNRAS, 459, 982
Shlosman I., Begelman M. C., 1987, Nature, 329, 810
Singer H. A., Olson P. L., 1984, Geophys. J., 78, 371
Steiman-Cameron T. Y., Durisen R. H., Boley A. C., Michael S., McConnell
C. R., 2013, ApJ, 768, 192
Toomre A., 1964, ApJ, 139, 1217
Turner N. J., Lee M. H., Sano T., 2014, ApJ, 783, 14
Vanon R., Ogilvie G. I., 2016, MNRAS, 463, 3725
Wallin B. K., Watson W. D., Wyld H. W., 1998, ApJ, 495, 774
Wei X., 2015, Geophys. Astrophys. Fluid Dyn., 109, 159
Wienker A., 2016, Master of Philosophy, Univ. Cambridge
A P P E N D I X A : BAC K G RO U N D E QU I L I B R I U M
We present in this appendix several tests to verify that our self-
gravity module in PLUTO and boundary conditions are correctly im-
plemented in 3D. Note that the 2D linear problem (both axisym-
metric and non-axisymmetric) has been already checked in Riols
& Latter (2016). In the presence of self-gravity, the vertical back-
ground equilibrium is governed by equations (7) and (8), where cs0
and ρ0 denote the sound speed and the density in the mid-plane. By
introducing the ‘isothermal’ Toomre parameter
Q2D0 =
cs0
πG
,
and the dimensionless z¯ = z/H0, ρ¯ = ρ/ρ0 as well as the ratio 
= 4H0ρ0/, the set of equations (7) and (8) reduce to a single
dimensionless equation
1
γ
d
dz¯
[
1
ρ¯
dρ¯γ
dz¯
]
+ 1 + 
Q2D0
ρ¯ = 0. (A1)
We numerically solve this equation by use of a finite-difference
method. We impose  = 1 and a certain value for Q2D0 . Imposing
these two quantities makes the calculation somewhat tricky because
the mid-plane density ρ0, and hence , is not known in advance.
It is then necessary to begin with a guess for ρ0 and to refine
the solution iteratively, using a shooting method, until the surface
density matches our prescription.
Note that the isothermal Toomre parameter defined here is not
exactly equal to the Q in the simulations, as the density-weighted
average cs differs from the mid-plane one. However, the sound speed
profile associated with the equilibrium scales as ρ1/3 with γ = 5/3,
so that the difference between Q2D0 and Q defined in Section 2.5 is
small.
Fig. A1 shows different density (and pressure) profiles for an
isothermal gas (γ = 1) and a polytropic gas with γ = 5/3, obtained
by numerical integration of equation (A1). When self-gravity is
included and Q is of the order of 1, the disc becomes more com-
pressed and its effective height, defined such that ρ = e−0.5ρ0, is
Heff = 0.52H0 for the isothermal case and Heff = 0.41H0 for the
polytrope .
For Q2D0 = 1, we use these semi-analytical density profiles (and
the associated pressure) as an initial condition in PLUTO. We ran
Figure A1. Density profiles of disc equilibria, for an isothermal gas (top)
and polytrope with γ = 5/3 (bottom; fixed  = 1). The solid blue line
corresponds to the theoretical profile computed from equation (A1) without
self-gravity (Q2D0 = ∞) while the dashed thick blue line is computed for
Q2D0 = 1. In each case, the thin red line is the density profile obtained with
PLUTO after 100−1 and is almost indiscernible from the theoretical one.
Figure A2. Vertical profile of the rms vertical velocity fluctuations at t =
100−1 excited by vertical boundary conditions.
the code for 100−1 and checked that the density has not evolved
during that time (see the red solid lines in Fig. A1). We also plot in
Fig. A2 the vertical profile of velocity fluctuations at t = 100−1.
For stable configurations, the fluctuation amplitude remains smaller
than 10−2cs0 . For a polytrope, these fluctuations can be larger than
the local sound speed, but do not influence the equilibrium in the
mid-plane. This shows that the boundary conditions are correctly
implemented and do not introduce too much noise, even for the
extreme case of a polytrope for which the density is floored to a
certain value in the upper atmosphere.
MNRAS 471, 317–336 (2017)
334 A. Riols, H. Latter and S.-J. Paardekooper
A PPENDIX B: LINEAR AXISYMMETRIC
M O D E S
To further test our Poisson solver in PLUTO, we simulate the lin-
ear solutions and compare them with semi-analytic calculations.
Unlike the 2D case, it is not possible to analytically derive a disper-
sion relation and we need to solve numerically the full eigenvalue
problem. We restrict our analysis to the isothermal case and as-
sume axisymmetric density and velocity perturbations of the form
[ρˆ(z), uˆ(z)] exp(ikxx − iωt). We denote by ρ = ρe(z) and e(z)
the density and potential equilibrium, respectively. We normalize
densities (background and perturbations) by the mid-plane density
ρ0, velocities by the uniform sound speed cs0 and pressure by ρc2s0 .
Finally, by introducing the dimensionless wave frequency ω¯ = ω/
and wavelength ¯kx = kxH0, the linearized and normalized system
of equations reads
− iω¯ρˆ + i¯kxρeuˆx + ddz¯ (ρeuˆz) = 0 (B1)
ρe(−iω¯uˆx − 2uˆy) = −i¯kx ˆP − i¯kxρe ˆ (B2)
ρe
(
−iω¯uˆy + 12 uˆx
)
= 0 (B3)
ρe(−iω¯uˆz) = −d
ˆP
dz¯
− ρˆ
(
z¯ + de
dz¯
)
− ρe d
ˆ
dz¯
(B4)
(
d2
dz¯2
− ¯k2x
)
ˆ = ρˆ
Q2D0
. (B5)
Solutions to this problem without self-gravity have been produced
by Lubow & Pringle (1993) for an isothermal stratified disc and
by Korycansky & Pringle (1995) and Ogilvie (1998) in the case of
a polytropic disc. More recently, the self-gravitating problem was
tackled by Mamatsashvili & Rice (2010).
The eigenfunctions in an isothermal non-self-gravitating disc are
the Hermite polynomials and the dispersion relation is
(−ω¯2 + n)(−ω¯2 + 1) − ¯k2xω¯2 = 0, (B6)
where n is the order of the Hermite polynomial and is related to the
vertical structure of the mode (the number of nodes in the eigen-
functions). There exist two branches of solutions for each n > 0, a
low-frequency branch with ω¯ < 1 that has an inertial character and
a high-frequency branch with ω¯ > 1 that has an acoustic character.
Note that the case n = 0 is often considered separately. At frequen-
cies larger than the epicyclic, it corresponds to a two-dimensional
acoustic–inertial wave, which can be associated with an f mode
(free surface oscillation) for the polytropic case.
We restrict our study to the isothermal case and a fixed ¯kx = 2π/5.
The surface density is set to  = 1. For a given value of the Toomre
parameter Q2D0 , we first compute the corresponding background
density and gravitational potential, using the methods described
in Appendix A (the equilibria fix the value of , which depends
on Q2D0 ). We then compute the eigenvalues ω¯i and corresponding
eigenmodes of the linearized system, using a Chebyshev collocation
method on a Gauss–Lobatto grid. This method results in a matrix
eigenvalue problem that can be solved using a QZ method (Golub
& Van Loan 1996; Boyd 2001). Numerical convergence is guar-
anteed by comparing eigenvalues for different grid resolutions and
eliminating the spurious ones. The vertical domain has a size 6H0
and is resolved by a maximum resolution of 600 points. Boundary
conditions for perturbed modes are ρeuˆ = 0 (momentum tends to 0)
and d ˆ/dzˆ = ± ˆkx ˆ (justified if the background density decreases
quasi-exponentially).
Figure B1. Properties of linear axisymmetric modes in the local isothermal
thin disc approximation. Top panel: dependence of the wave frequency ω
on Q2D0 for the inertial branch modes n = 1, 2, . . . , 7. Bottom panel:
growth rates of the first unstable n = 0 mode. In all cases, ¯kx = 2π/5.
Solid curves are theoretical growth rates or frequencies, calculated with the
linear eigensolver. Red and green stars are the values obtained from PLUTO
simulations, with resolution of 256 and 512 points per azimuthal wavelength,
respectively. The dotted line is Q2D = 1 while the dashed line represents the
critical Q2D for which the n = 0 mode becomes unstable.
When self-gravity is negligible, Q2D0 = ∞, we checked that the
eigenvalues match the dispersion relation (B6) and that for n > 0
there exist two different types of solutions corresponding to either
ω¯ > 1 (p modes) or ω¯ < 1 (r modes). We then vary Q2D0 and focus
first on the modes with inertial character. For each branch labelled n
= 1, 2, . . . , 7 (n defined as the order of the Hermite polynomial in the
limitQ2D0 = ∞), we plot in the top panel of Fig. B1 the evolution of
the wave frequency as a function of Q2D0 . We find that ω¯ decreases
with increasing Q2D0 but never becomes complex (unstable), at
least for the range of Q2D0 studied (0.25 < Q2D0 < ∞). Instead, the
gravitational unstable mode is associated with the inertial–acoustic
branch n = 0 (f modes). The corresponding eigenfunction has no
structure in z for Q2D0 = ∞ but appears to develop some structure
as Q2D0 is decreased. Fig. B2 shows in particular that the vertical
velocity component of the unstable mode develops an odd symmetry
as for n = 1 .
Having solved the linear problem directly, we check that our
version of PLUTO reproduces its main results in simulations of stable
MNRAS 471, 317–336 (2017)
GI and turbulence in astrophysical discs 335
Figure B2. Vertical structure of the n = 0 eigenmode for three different
Q2D0 = {∞, 1, 0.4}. Left to right: density and vertical velocity perturba-
tions. Note that the mode is normalized such that ρˆ(z = 0) = 1. The noise
at the boundary (in particular for Q2D0 = 0.4) results from the fact that we
are solving the linear problem for the momentum and not for the velocity
components directly.
Figure B3. The top two panels describe the time evolution of the density
and radial velocity perturbation of the n = 1 inertial mode for different
values of Q2D0 [amplitudes of perturbations are expressed in terms of the
Hermite norm, defined through equation (B7)]. The last panel shows the
evolution of the average kinetic energy for the n = 0 mode for four values
of Q2D0 .
and unstable axisymmetric waves. For that purpose, we used a
box of size 5H × 5H × 6H (full disc in z) and introduced at t =
0 a perturbation with kx = 2π/Lx , ky = 0 whose vertical shape
is determined by the eigensolver (in addition to the background
density equilibrium). Fig. B3 shows the density and radial velocity
amplitude of the n = 1 inertial mode, expressed in terms of the
Fourier–Hermite norm,
An[X] =
∫ ∫
ˆX(z) cos(kxx)Hen(z)e−z2/2dxdz (B7)
for Q2D0 = 1 and 10. In the above Hen is the Hermite polynomial
of order n. At a resolution of 256 × 1 × 256, the solution remains
periodic and conserves its vertical shape for more than 200−1. The
wave frequencies, inferred from these plots, are superimposed on
the dispersion diagram of Fig. B1 and compared with the theoretical
frequencies, for different Q2D0 (red stars denote a resolution of 256
× 1× 256 whereas green stars are for a resolution of 512× 1× 512).
We find a very good agreement between the frequencies computed
with the eigensolver and those measured from the simulations. The
relative error remains always smaller than 1 per cent (typically the
mean error is 0.3 per cent). We have undertaken a similar test for
the n = 2 branch, although only two points have been computed in
that case.
Finally, we simulated the n = 0 branch that becomes unstable for
Q2D0 = 0.63. Fig. B3 (bottom) shows the evolution of the averaged
kinetic energy in the box for different Q2D0 . Growth rates are easily
measured by computing the slope of these curves. In Fig. B1 (bottom
panel), we superimpose the numerical growth rate obtained (green
stars) with the predicted growth rate as a function ofQ2D0 . We found
again that numerical and theoretical growth rates match remarkably
well, even for small values of Q2D0 .
APPENDI X C : N ON-LI NEAR ‘I NSTABI LITY ’
O F A X I S Y M M E T R I C M O D E S
Simulations of gravitoturbulence in shear flows exhibit strong large-
scale axisymmetric motions (see Section 4). To better understand
how these motions are generated by the turbulent flow, we perform
a simple experiment: we considered a gravitoturbulent state in a
small box Lx = 20H0 and we used the symmetry of the shearing
box to construct a stacked version of this state in a box of twice the
original size. We then examined the mid-plane time evolution of
the fundamental axisymmetric mode kx = kx0 , which possesses a
wavelength with the size of the new box. The results are summarized
in Fig. C1. At t = 0, the energy associated with this mode is zero,
since we started from a state computed in box of half size. The
mode is forced by non-linearities and grows with a typical growth
rate γ a = 0.3. Its short-term evolution is very stochastic, but its
long-term evolution is more or less exponential .
To understand the nature of this instability, we investigated the
kinetic energy budget of this fundamental mode:
1
2
duF 2
dt
= 3
2
uxF uyF + uF · ∇ − uF ·
(∇P
ρ
)
F
− uF · (u · ∇u)F . (C1)
The subscript ‘F’ denotes here the fundamental axisymmetric mode
kx = kx0 or its related projection for non-linear terms. On the right-
hand side, the two first linear terms are associated with the back-
ground shear and self-gravity. The third and fourth terms are related
to the pressure gradient and non-linear advection. The pressure gra-
dient term can be decomposed into a linear and non-linear part:(∇P
ρ
)
F
=
(∇P
ρ
)NL
F
+
(∇PF
ρ0
)
, (C2)
where ρ0 is the mid-plane background density.
Fig. C1 (bottom) shows the cumulative contribution of each term
involved in the mode growth. First, we found that the work done
by self-gravity is mainly negative and hence does not participate
in the instability. This is expected since the axisymmetric mode
is gravitationally linearly stable. The shear term −Sux0uy0 is also
stabilizing. Note that the Coriolis force does no work since it only
MNRAS 471, 317–336 (2017)
336 A. Riols, H. Latter and S.-J. Paardekooper
Figure C1. Top: time evolution of the specific kinetic energy in the mid-
plane z = 0 of the fundamental axisymmetric mode kx = kx0 in the box Lx =
40H0. Bottom: energy budget of the fundamental axisymmetric perturbation.
Each curve represents the cumulative work associated with a term (or a force)
in the Navier–Stokes equation.
redistributes energy from the y component to the x component or
vice versa.
The main contribution to the kinetic energy, at least for the first
orbits, is clearly the non-linear term associated with the pressure
work. Physically, this term corresponds to a baroclinic effect, due
to the misalignment of pressure and density gradients. Indeed, if the
pressure term is assumed to be dominant in the energy evolution of
the axisymmetric mode, then the vorticity equation reduces to
∂ω
∂t
= ∇P ×∇ρ
ρ2
. (C3)
The formation of axisymmetric cells associated with vorticity
aligned in the y direction is possible if baroclinic non-axisymmetric
perturbations (or spiral waves) interact with each other. In other
words, if the average in y of (∇P ×∇ρ)F · yˆ has a positive feed-
back.
A possible origin for the global process is a triadic interaction
between the m = 1 spiral density waves and axisymmetric modes.
The scenario can be described as follows: non-axisymmetric waves
tap into two reservoirs of energy that are the shear and the self-
gravity. The non-linear interaction between a swinging leading wave
and a dying trailing wave produces an axisymmetric modulation
(through the baroclinic term), which in turn reinforces the non-
axisymmetric waves, via a mechanism similar to that described in
Lithwick (2007). This scenario provides a potential explanation for
why axisymmetric modes are growing exponentially, despite the
absence of a classical ‘linear’ instability.
This paper has been typeset from a TEX/LATEX file prepared by the author.
MNRAS 471, 317–336 (2017)