MNRAS 462, 4549–4554 (2016) doi:10.1093/mnras/stw2172
Advance Access publication 2016 August 4
On the survival of zombie vortices in protoplanetary discs
Geoffroy R. J. Lesur1,2‹ and Henrik Latter3
1CNRS, IPAG, F-38000 Grenoble, France
2Univ. Grenoble Alpes, IPAG, F-38000 Grenoble, France
3DAMTP, University of Cambridge, CMS, Wilberforce Road, Cambridge CB3 0WA, UK
Accepted 2016 August 2. Received 2016 July 31; in original form 2016 June 3
ABSTRACT
Recently it has been proposed that the zombie vortex instability (ZVI) could precipitate
hydrodynamical activity and angular momentum transport in unmagnetized regions of proto-
planetary discs, also known as ‘dead zones’. In this Letter we scrutinize, with high-resolution
3D spectral simulations, the onset and survival of this instability in the presence of viscous and
thermal physics. First, we find that the ZVI is strongly dependent on the nature of the viscous
operator. Although the ZVI is easily obtained with hyperdiffusion, it is difficult to sustain
with physical (second order) diffusion operators up to Reynolds numbers as high as 107. This
sensitivity is probably due to the ZVI’s reliance on critical layers, whose characteristic length-
scale, structure, and dynamics are controlled by viscous diffusion. Second, we observe that
the ZVI is sensitive to radiative processes, and indeed only operates when the Peclet number
is greater than a critical value ∼104, or when the cooling time is longer than ∼10−1. As
a consequence, the ZVI struggles to appear at R  0.3 au in standard 0.01 M T Tauri disc
models, though younger more massive discs provide a more hospitable environment. Together
these results question the prevalence of the ZVI in protoplanetary discs.
Key words: hydrodynamics – instabilities – protoplanetary discs.
1 IN T RO D U C T I O N
The origin of angular momentum transport in accretion discs, espe-
cially protoplanetary discs, is a long-standing issue in the astrophys-
ical community. Angular momentum transport is the mechanism
that governs the global dynamics of the gas, in particular its accre-
tion on to the central star. It is therefore especially important if one
is to predict the long-term evolution and structure of protoplanetary
discs.
The magnetorotational instability (MRI; Balbus & Hawley 1991)
is believed to be the main driver of angular momentum transport in
accretion discs. By sustaining three-dimensional magnetohydrody-
namics (MHD) turbulence, the MRI transports angular momentum
outwards and leads to mass accretion at rates compatible with ob-
servations. It is far from assured, however, that cold protoplanetary
discs are sufficiently ionized to sustain MHD turbulence. This has
led to the concept of ‘dead zones’ (Gammie 1996), internal regions
of the disc where the MRI is quenched. The question of angular
momentum transport in dead zones is highly debated, and in which
hydrodynamical instabilities are likely to be key (Turner et al. 2014).
The radial Keplerian rotation profile of astrophysical discs is
known to be hydrodynamically stable, both linearly and non-linearly
(Lesur & Longaretti 2005; Edlund & Ji 2014). However, additional
 E-mail: geoffroy.lesur@univ-grenoble-alpes.fr
physics, such as cooling, heating, and stratification, could unleash
new hydrodynamical instabilities. In recent years, several have been
identified, including the subcritical baroclinic instability (SBI; Pe-
tersen, Julien & Stewart 2007; Lesur & Papaloizou 2010), the ver-
tical shear instability (VSI; Nelson, Gressel & Umurhan 2013), the
convective overstability (Klahr & Hubbard 2014), and more recently
the zombie vortex instability (ZVI) which appears in rotating shear
flows exhibiting a stable vertical stratification.
The ZVI was first observed (but not clearly identified as such)
in the anelastic simulations of Barranco & Marcus (2005) and was
subsequently isolated by Marcus et al. (2013) using Boussinesq
spectral simulations. This instability, of non-linear nature, produces
‘self-replicating’ vortices thanks to the excitation of very thin criti-
cal layers. The ZVI also appears in compressible simulations with
various initial conditions. It has been proposed, but not yet demon-
strated, that the excitation of spiral density waves by zombie vortices
could lead to significant angular momentum transport in dead zones
(Marcus et al. 2015), thereby solving the angular momentum trans-
port problem in these regions. However, the physical mechanism
driving the instability remains mysterious. The existence and ex-
citation of critical layers by a perturbation is a well-known linear
mechanism in shear flows (Drazin & Reid 1981), but their non-
linear saturation and spontaneous transformation into new vortices
is largely unexplained. Recently, based on a linear and a quasi-linear
analysis, Umurhan, Shariff & Cuzzi (2016) proposed that critical
layers could be subject to a stratified variant of the Rossby wave
C© 2016 The Authors
Published by Oxford University Press on behalf of the Royal Astronomical Society
 at U
niversity of Cam
bridge on Septem
ber 28, 2016
http://m
nras.oxfordjournals.org/
D
ow
nloaded from
 
4550 G. R. J. Lesur and H. Latter
instability (Lovelace et al. 1999). However, this scenario needs to
be supplemented with a mechanism to self-sustain the original per-
turbation, and further tested with more realistic non-linear models
of the layer, including viscous effects. Indeed, because diffusive
physics determines the layer’s structure and evolution, the nature
of viscosity (and whether it is physical or numerical) should be a
fundamental ingredient in any theory.
In this Letter, we scrutinize the foundations of the ZVI with
our focus squarely on the role of viscosity and cooling. We first
present our physical model and numerical methods, which are very
similar to Marcus et al. (2013). We then move to the question of
the physical convergence of the ZVI as a function of the viscous
operator. We also explore the dependence of the ZVI on cooling,
and compare our results to realistic protoplanetary disc models. We
finally summarize our results and propose future routes of research
into the ZVI.
2 M E T H O D S
2.1 Physical model
We represent the local dynamics of the disc with the shearing box
approximation. To further simplify the dynamics, we employ incom-
pressibility but include vertical buoyancy effects via the Boussinesq
approximation. In this framework, the equations of motion read
∂u
∂t
+ u · ∇u = −∇ − 2 × u + 2Sxex
−N2θez + νu + (ν6∇2)3u, (1)
∂θ
∂t
+ u · ∇θ = uz + χθ + (χ6∇2)3θ − θ
tc
, (2)
∇ · u = 0. (3)
In the above formulation we have defined the local rotation fre-
quency , the shear rate S, the Brunt–Vaissala frequency N, and the
generalized pressure , which allows us to satisfy the incompress-
ible condition (3). In addition, we have introduced several explicit
diffusion operators: the usual second-order viscosity ν and thermal
diffusivity χ are supplemented with sixth-order ‘hyperdiffusion’
operators with coefficients ν6 and χ6. Hyperdiffusion has no real
physical motivation but can be useful numerically to reduce diffu-
sion on large scales without accumulating energy at the grid scale.
Note that a similar hyperdiffusion operator was used by Marcus
et al. (2015). Finally, we have added a Newtonian cooling in the
form of a constant thermal relaxation time tc.
The above set of equations admits a simple solution of pure shear
flow u0 = −Sxey . In the following, we define perturbations (not
necessarily small) to this global shear flow v = u − u0.
The equations of motions are supplemented by a set of periodic
boundary conditions in the y and z directions. In the x direction, we
use shear-periodic boundary conditions, following Hawley, Gam-
mie & Balbus (1995).
2.2 Dimensionless numbers and units
The set of equations above includes several dynamical time-scales
that can be usefully compared via appropriate dimensionless num-
bers. We define and use the following ones.
(i) The Rossby number q = S/. In Keplerian accretion discs q
= 3/2, which we will be the framework of this Letter.
(ii) The Froude number Fr = N/. In this work, we will always
assume Fr = 2 which corresponds to the fiducial case studied by
Marcus et al. (2013) of a moderately stratified flow. Note however
that expected Froude numbers in protoplanetary discs are somewhat
lower than this value, with Fr  0.3 (Dubrulle et al. 2005). Our
set-up therefore represents an upper bound on the amplitude of
stratification effects.
(iii) The Reynolds number Re = L2/ν compares the ampli-
tude of non-linear advection terms to viscous diffusion. Equiva-
lently, we define a Reynolds number based on hyperviscosity Re6
= 1/3L2/ν6.
(iv) The Peclet number Pe = L2/χ compares non-linear ad-
vection to thermal diffusion. As for the Reynolds number, we also
define a hyperdiffusion Peclet number Pe6.
(v) The dimenionsless cooling time τ = tc.
Unless mentioned otherwise, we use −1 as our time unit and
the box size L as our length unit.
2.3 Numerical technique
We employ SNOOPY to integrate the equations of motion. SNOOPY is a
spectral code using a Fourier decomposition of the flow to compute
spatial derivatives. Time integration is performed using a low stor-
age third-order Runge–Kutta scheme. Diffusive operators are solved
by an implicit operator which maintains the third-order accuracy of
the scheme. To avoid spectral aliasing due to the quadratic non-
linearities, we use a standard 2/3 anti-aliasing rule when computing
each non-linear term. The code and ZVI set-up is freely available
on the author’s web site.
In this Letter, we use two sets of initial condition: single vortex
initial conditions (runs label ending with ‘v’) and Kolmogorov-like
noise (runs label ending with ‘k’).
Our single vortex initial condition is similar to Marcus et al.
(2013) with an isolated Gaussian vortex centred at the origin of
the box with a size σ and a velocity amplitude v0. The initial
perturbation reads
v0x(x) = yv0/σ exp
[−(x2 + y2 + z2)/σ 2] ,
v0y(x) = −xv0/σ exp
[−(x2 + y2 + z2)/σ 2] .
Our simulations start with σ = 0.07 and v0 = 0.03 in order to get
results close to Marcus et al. (2013). This perturbation corresponds
to a stratified anticyclonic vortex with a vertical vorticity ωz =
∂xvy − ∂yvx  −0.8.
When using Kolmogorov-like noise, we randomly excite each
velocity wavenumber isotropically in phase and amplitude and set
the energy spectrum to E(k) ∝ k−5/3. We normalize our initial con-
ditions so that
√
〈v2〉 = 4 × 10−2 at t = 0.
3 R ESULTS
3.1 Fiducial case and hyperdiffusion
We first introduce our fiducial model h-256 (see Table 1) which es-
sentially reproduces the results of Marcus et al. (2013). We choose a
resolution of 2563 Fourier modes in a cubic box representing a Kep-
lerian disc with q = 3/2 and Fr = 2. Neither viscosity nor diffusion
is imposed, ν = χ = 0. We instead use sixth-order hyperdiffusion to
dissipate energy at small scales that would otherwise accumulate,
spectral codes being inherently energy conserving schemes. We set
MNRAS 462, 4549–4554 (2016)
 at U
niversity of Cam
bridge on Septem
ber 28, 2016
http://m
nras.oxfordjournals.org/
D
ow
nloaded from
 
ZVI in protoplanetary discs 4551
Table 1. List of simulations discussed in this Letter.
Run Re Pe Re6 Pe6 τ Resolution ZVI
h-256-v ∞ ∞ 5 × 105 5 × 105 ∞ 2563 Yes
h-1024-v ∞ ∞ 5 × 105 5 × 105 ∞ 10243 Yes
v6-1024-v 106 106 ∞ ∞ ∞ 10243 No
v7-1024-v 107 107 ∞ ∞ ∞ 10243 No
d5-256-v ∞ 6.4 × 105 5 × 105 ∞ ∞ 2563 Yes
d4-256-v ∞ 1.0 × 105 5 × 105 ∞ ∞ 2563 Yes
d4-256-k ∞ 8.0 × 104 5 × 105 ∞ ∞ 2563 Yes
d3-256-v ∞ 5.0 × 104 5 × 105 ∞ ∞ 2563 Yes
d3-256-k ∞ 4.0 × 104 5 × 105 ∞ ∞ 2563 Yes
d2-256-v ∞ 2.0 × 104 5 × 105 ∞ ∞ 2563 No
d2-256-k ∞ 2.0 × 104 5 × 105 ∞ ∞ 2563 No
d1-256-v ∞ 1.0 × 104 5 × 105 ∞ ∞ 2563 No
d1-256-k ∞ 1.0 × 104 5 × 105 ∞ ∞ 2563 No
t5-256-v ∞ ∞ 5 × 105 ∞ 128 2563 Yes
t4-256-v ∞ ∞ 5 × 105 ∞ 64 2563 Yes
t4-256-k ∞ ∞ 5 × 105 ∞ 64 2563 Yes
t3-256-v ∞ ∞ 5 × 105 ∞ 32 2563 Yes
t3-256-k ∞ ∞ 5 × 105 ∞ 32 2563 Yes
t2-256-v ∞ ∞ 5 × 105 ∞ 16 2563 No
t2-256-k ∞ ∞ 5 × 105 ∞ 16 2563 No
Figure 1. Vertical vorticity ωz in a x–z cut of our fiducial simulation with
Re6 = Pe6 = 5 × 105 at t = 500. Similarly to Marcus et al. (2013), we
observe the formation and replication of anticyclonic vortices on a fixed
lattice.
Re6 = Pe6 = 5 × 105. This simulation allows us to reproduce the
main results of Marcus et al. (2013): self-replicating vortices on a
fixed lattice (Fig. 1), and a growth in kinetic energy EK ≡ 〈v2/2〉
associated with these vortices (Fig. 2, blue line). As expected, new
vortices appear at critical layers defined, from the initial vortex, by
xc = ±FrL/(2πmq)  ±0.21/m, where m is the azimuthal mode
number.
More interesting is the behaviour of these simulations when res-
olution and dissipation processes are modified. To illustrate this,
let us consider higher resolution simulations with 10243 Fourier
modes. We first perform a resolution test (h-1024) to reproduce our
fiducial run with hyperdiffusion which confirms that our 2563 run
is numerically converged, at least with respect to EK (Fig. 2, green
line). We then restart this high-resolution simulation at t = 400 but
revert to classical dissipation coefficients. We consider two cases:
Re = Pe = 106 (v6-1024) and Re = Pe = 107 (v7-1024). The Re =
106 shows a clear and steep decay indicating that the ZVI disappears
for this Reynolds number. If we move up to Re = 107, a decline is
still seen, but we cannot say for sure that the ZVI is deactivated. A
Figure 2. Volume averaged kinetic energy as a function of time for several
simulations at Fr = 2 and q = 3/2. The blue curve corresponds to our
fiducial case. High-resolution runs including diffusion are restarted from
the hyperdiffusion run at t = 400.
careful examination of the critical layers at t = 600 in the Re = 107
case shows that they are resolved by only four to five collocation
points (Fig. 3). We therefore conclude that the critical Reynolds
number Rec for the ZVI (if it exists) is certainly larger than 106, and
possibly larger than 107. Simulations with at least 20483 (or even
40963) points will be required to confirm the existence of the ZVI
with second-order dissipation operators.
3.2 Cooling and heating processes
The sensitivity to the Reynolds number indicates that the ZVI mech-
anism is highly dependent on dissipation and diffusion. In proto-
planetary discs, Reynolds numbers are huge, so a high Rec is not
physically a problem (although it is definitely a problem for numer-
ical simulations). Cooling in these discs, on the other hand, is far
from negligible. It is therefore desirable to test the existence of the
ZVI in the presence of cooling and heating.
In protoplanetary discs, cooling and heating are dominated
by radiative transfer, thermal conductivity being unimportant.
MNRAS 462, 4549–4554 (2016)
 at U
niversity of Cam
bridge on Septem
ber 28, 2016
http://m
nras.oxfordjournals.org/
D
ow
nloaded from
 
4552 G. R. J. Lesur and H. Latter
Figure 3. Vertical velocity as a function of x measured at z = 0.1 in run
v7-1024-v at t = 600. The inset zooms on the critical layer at x  0.21 with
dots representing the spectral collocation points. Note the sharpness of the
critical layers despite the large resolution.
In optically thick regions, if the length-scales λ under consider-
ation are larger than the photon mean free path 
ph = 1/κρ, where
κ is the opacity and ρ is the gas density (i.e. the disc is optically thick
on scale λ), cooling and heating can be approximated by thermal
diffusion with a diffusion coefficient
χ = 16σT
3
3κρ2cv
, (4)
whereσ is the Stefan–Boltzmann constant and cv is the heat capacity
of the gas, which we will assume to be diatomic. In the opposite
limit λ < 
ph, radiative cooling acts like a scale-free Newtonian
cooling with a characteristic time-scale tc,
tc =

2ph
3χ
. (5)
Note that this cooling time-scale is not the same as the global cool-
ing time-scale of a vertically integrated disc (subject to external
heating and radiative cooling), since here we look at small-scale
thermal perturbations embedded in an optically thick medium.
Similar estimates for the cooling function were obtained by Lin
& Youdin (2015) for the VSI.
To illustrate the typical Peclet number and cooling times in pro-
toplanetary discs, we have considered a typical T Tauri disc model
 = 140R−1au g cm−2 and T = 280 R−1/2au K which corresponds to a
0.01 M mass disc extending to 100 au. We assume the disc to be
vertically isothermal as we do not solve the full radiative transfer
equations. Rosseland opacities including gas and dust contributions
are obtained from Semenov et al. (2003) assuming spherical homo-
geneous dust grains of solar composition. To compute the resulting
Peclet number, we have identified the box scale L to the disc pres-
sure scale height1 H ≡ cs/, where cs is the local sound speed
and  the local Keplerian frequency. The resulting map for thermal
diffusion (Peclet number) is shown in Fig. 4. As mentioned above,
the thermal diffusion approximation is valid only for scales λ >

ph. The typical photon mean free path 
ph is shown in Fig. 5. The
smallest 
ph/H is found close to the midplane in the inner parts of
the disc. These are the regions expected to be well described by the
1 In principle L can be arbitrarily smaller than H since we work in the
incompressible limit. We have not considered this case since it leads to
critical disc Pe even larger than the one discussed here, leading to a smaller
domain of existence for the ZVI.
Figure 4. Peclet number as a function of position in a 0.01 M disc model.
The magenta contour defines the τ = 1 surface below which the disc is
optically thick: 
ph < H.
Figure 5. Photon mean free path 
ph compared to the disc scale height H
in a 0.01 M disc model. Shortest mean free paths are found close to the
midplane in the innermost parts of the disc.
thermal diffusion approximation on most relevant scales.2 On the
contrary, the outer regions R > 10 au have 
ph  H. In these regions,
cooling is best described by a constant cooling time characterized
by the dimensionless parameter τ ≡ tc. Since the cooling time
(equation 5) does not depend on density, τ is only a function of
radius, shown in Fig. 6. From these three figures, we deduce that
for R  1 au, Pe < 103, and τ < 10−2. Note that in this discussion,
we have only considered optically thick regions of the disc. Opti-
cally thin regions will likely have longer cooling time. However,
they will also be affected by ionizing radiation which will allow
MHD processes such as the MRI or winds to occur, making the ZVI
irrelevant in these regions.
Our last task is to test in which parameter regime the ZVI lives.
To this end, we have performed a set of simulations identical to
our fiducial simulation, except that thermal hyperdiffusion is now
replaced by a classical thermal diffusion operator with Pe ∈ [104,
106] (runs dxxxx) or by a fixed cooling parameter τ ∈ [10, 200] (runs
txxxx). We have used either the Gaussian vortex initial conditions
(runs ending with ‘v’) or Kolmogorov noise initial conditions (runs
ending with ‘k’). The energy evolution of the simulations starting
2 Note however that the thickness of the critical layers involved in the ZVI
can be several orders of magnitude smaller than the disc scale as it is set
by gas molecular viscosity. It is therefore possible that critical layers are
always in the optically thin regime for realistic Reynolds numbers.
MNRAS 462, 4549–4554 (2016)
 at U
niversity of Cam
bridge on Septem
ber 28, 2016
http://m
nras.oxfordjournals.org/
D
ow
nloaded from
 
ZVI in protoplanetary discs 4553
Figure 6. Dimensionless cooling time as a function of radius in a 0.01 M
disc model. Jumps are due to the condensation of molecules on to dust grains
which abruptly change opacities (Semenov et al. 2003).
Figure 7. Volume averaged kinetic energy for runs dx-256-v varying ther-
mal diffusivities. Pe > 2 × 104 is needed to sustain the ZVI.
Figure 8. Volume averaged kinetic energy for runs tx-256-v varying ther-
mal relaxation time-scales. τ > 16 is needed to sustain the ZVI.
with a Gaussian vortex (Figs 7 and 8) clearly indicates that the
ZVI requires Pe > 2 × 104 and τ > 16. Runs with τ ≤ 16 or
Pe ≤ 104 becomes axisymmetric at t ∼ 2000 which ensure that
the ZVI is definitely switched off for this range of parameters.
Very similar limits are obtained when using Kolmogorov noise as
an initial condition (Table 1). Our limits of existence for the ZVI
therefore do not depend strongly on the chosen initial condition.
These dimensionless numbers are clearly excluded in our typical
disc model presented above, except maybe in the diffusive regime
in the innermost regions (R  0.1 au) which are also likely to be
unstable to the MRI due to their proximity to the central star (Latter
& Balbus 2012).
4 C O N C L U S I O N S
In this Letter, we have explored the sensitivity of the ZVI to dif-
fusive and thermal processes. We find that one can easily produce
this instability with hyperdiffusion operators, but not with classical
viscous operators. We conjecture that a resolution of at least 20483
collocation spectral points and a Reynolds number higher than 107
are required to ascertain the presence of the ZVI with physical
dissipation. This should not come as a surprise since the instabil-
ity mechanism relies on the physics of buoyancy critical layers,
which are themselves controlled by diffusion (the process that sets
their characteristic length-scale). It is therefore essential to prop-
erly resolve these structures with realistic dissipation operators (i.e.
neither hyperdiffusion nor numerical dissipation). Note that finding
the ZVI with finite volume codes does not solve this issue since
these codes are also strongly affected by numerical diffusion.
We have also explored the sensitivity of the ZVI to cooling. If
radiative diffusion or Newtonian cooling is too efficient then the
action of buoyancy is diminished, as expected, and the instability
switches off. The critical Peclet number below which ZVI fails is
∼104, while the critical cooling time is ∼10−1. This critical Pe
have been obtained with a fixed hyperdiffusivity so that the viscous
scale is always much smaller than the thermal diffusion scale, as
in a real protoplanetary disc. However, the ZVI may also show a
dependence on Pr = ν/χ or other combinations of dimensionless
parameters. These dependencies have not been explored in this
work.
Using a typical T Tauri disc model of 0.01 M mass, we find
that the ZVI may struggle to survive except in the densest and
innermost regions of the disc (R ∼ 0.1 au) which are in any case
MRI unstable. This is true whether the characteristic length-scale of
the ZVI falls in the diffusive or Newtonian cooling regimes. Taken
on face value, these results cast doubt on the ZVI as a potential
source of turbulent transport and vortices in late-type objects (class
II). We note, however, that younger discs (M ∼ 0.1 M) or discs
with steeper density profiles such as the minimum mass solar neb-
ula (Hayashi 1981) could reach Pe ∼ 105 at R  1 au thanks to
the increase in gas density. More massive discs would be subject
to gravitational instabilities in their outer part but could be ZVI
unstable in their inner part. Nevertheless, these scenarios must be
confirmed by (a) demonstrating the existence and convergence of
the ZVI with explicit viscous dissipation and (b) including a proper
radiative transfer modelling to compute cooling accurately.
Note finally that this work has been performed in a local approxi-
mation (constant stratification, incompressibility, constant cooling).
The ZVI being a local instability (Marcus et al. 2013, 2015), it is
well captured and described by this model. Our study does not ex-
clude the possibility of a global instability which would be due
to the vertical structure of the disc. However, such a hypothetical
instability would be driven by a different physical mechanism than
that of the ZVI. Note also that global simulations will inherit (and
probably exacerbate) the ZVI’s numerical convergence problem (cf.
Section 3.1).
AC K N OW L E D G E M E N T S
The computations presented here were performed using the
Froggy platform of the CIMENT infrastructure (https://ciment.
ujf-grenoble.fr). HL acknowledges funding from STFC grant
MNRAS 462, 4549–4554 (2016)
 at U
niversity of Cam
bridge on Septem
ber 28, 2016
http://m
nras.oxfordjournals.org/
D
ow
nloaded from
 
4554 G. R. J. Lesur and H. Latter
ST/L000636/1 and helpful advice from John Papaloizou, Michael
McIntyre, and Steve Lubow.
R E F E R E N C E S
Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214
Barranco J. A., Marcus P. S., 2005, ApJ, 623, 1157
Drazin P. G., Reid W. H., 1981, NASA STI/Recon Technical Report A, 82
Dubrulle B., Marie´ L., Normand C., Richard D., Hersant F., Zahn J.-P., 2005,
A&A, 429, 1
Edlund E. M., Ji H., 2014, Phys. Rev. E, 89, 021004
Gammie C. F., 1996, ApJ, 457, 355
Hawley J. F., Gammie C. F., Balbus S. A., 1995, ApJ, 440, 742
Hayashi C., 1981, Progress Theor. Phys. Suppl., 70, 35
Klahr H., Hubbard A., 2014, ApJ, 788, 21
Latter H. N., Balbus S., 2012, MNRAS, 424, 1977
Lesur G., Longaretti P.-Y., 2005, A&A, 444, 25
Lesur G., Papaloizou J. C. B., 2010, A&A, 513, A60
Lin M.-K., Youdin A. N., 2015, ApJ, 811, 17
Lovelace R. V. E., Li H., Colgate S. A., Nelson A. F., 1999, ApJ, 513, 805
Marcus P. S., Pei S., Jiang C.-H., Hassanzadeh P., 2013, Phys. Rev. Lett.,
111, 084501
Marcus P. S., Pei S., Jiang C.-H., Barranco J. A., Hassanzadeh P., Lecoanet
D., 2015, ApJ, 808, 87
Nelson R. P., Gressel O., Umurhan O. M., 2013, MNRAS, 435, 2610
Petersen M. R., Julien K., Stewart G. R., 2007, ApJ, 658, 1236
Semenov D., Henning T., Helling C., Ilgner M., Sedlmayr E., 2003, A&A,
410, 611
Turner N. J., Fromang S., Gammie C., Klahr H., Lesur G., Wardle M., Bai
X.-N., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T.,
eds, Protostars and Planets VI. Univ. Arizona Press, Tuscan, AZ, p. 411
Umurhan O. M., Shariff K., Cuzzi J. N., 2016, ApJ, in press, preprint
(arXiv:1601.00382)
This paper has been typeset from a TEX/LATEX file prepared by the author.
MNRAS 462, 4549–4554 (2016)
 at U
niversity of Cam
bridge on Septem
ber 28, 2016
http://m
nras.oxfordjournals.org/
D
ow
nloaded from