# Convective overstability in radially global protoplanetary disks I - Pure
gas dynamics

Marius Lehmann<EMAIL_ADDRESS>Institute of Astronomy and
Astrophysics, Academia Sinica, Taipei 10617, Taiwan  Min-Kai Lin Institute of
Astronomy and Astrophysics, Academia Sinica, Taipei 10617, Taiwan  Physics
Division, National Center for Theoretical Sciences, Taipei 10617, Taiwan

###### Abstract

Protoplanetary disks are prone to several hydrodynamic instabilities. One
candidate, Convective Overstability (COS), can drive radial semi-convection
that may influence dust dynamics and planetesimal formation. However, the COS
has primarily been studied in local models. This paper investigates the COS
near the mid-plane of radially global disk models. We first conduct a global
linear stability analysis, which shows that linear COS modes exist only
radially inward of their Lindblad resonance (LR). The fastest-growing modes
have LRs near the inner radial domain boundary with effective radial
wavelengths that can be a substantial fraction of the disk radius. We then
perform axisymmetric global simulations and find that the COS’s nonlinear
saturation is similar to previous incompressible shearing box simulations. In
particular, we observe the onset of persistent zonal and elevator flows for
sufficiently steep radial entropy gradients. In full 3D, non-axisymmetric
global simulations, we find the COS produces large-scale, long-lived vortices,
which induce outward radial transport of angular momentum via the excitation
of spiral density waves. The corresponding $\alpha$-viscosity values of order
$10^{-3}$ agree well with those found in previous 3D compressible shearing box
simulations. However, in global disks, significant modifications to their
radial structure are found, including the formation of pressure bumps.
Interestingly, the COS typically generates an outward radial mass transport,
i.e. decretion. We briefly discuss the possible implications of our results
for planetesimal formation and for interpreting dust rings and asymmetries
observed in protoplanetary disks.

Protoplanetary disks (1300); Hydrodynamics (1963); Astrophysical fluid
dynamics (101); Planet formation (1241); Planetesimals (1259)

##  1 Introduction

The proposed formation of km-sized planetesimals in protoplanetary disks
(PPDs) via the core accretion model (Safronov, 1972) faces several obstacles,
such as the so-called bouncing and drift barriers, inhibiting efficient
coagulation of dust grains (e.g. Blum 2018; Drazkowska et al. 2022). In order
to overcome these, mechanisms to effectively concentrate dust grains in the
mid-plane of PPDs are required. The most popular such mechanism is the
streaming instability (SI: Youdin & Goodman 2005; Youdin & Johansen 2007;
Johansen & Youdin 2007), which in its nonlinear saturation can produce dense
clumps of dust particles, which can subsequently collapse under the action of
selfgravity to form planetesimals (Johansen & Youdin, 2007; Bai & Stone, 2010;
Simon et al., 2016; Carrera et al., 2021). However, for small dust grains the
clumping by the SI is strong enough only if the metallicity (or dust
abundance) takes at least a few times the solar value (Li & Youdin 2021).

In the recent years, the possible occurrence of purely hydrodynamic
instabilities in PPDs has received an increased amount of attention. The main
reason is that the magneto-rotational instability (MRI: Balbus & Hawley 1991),
which had long been considered to be the primary source for mass accretion in
PPDs, is likely weak or absent in large parts of PPDs, due to low ionization
rates (Gammie, 1996; Turner & Drake, 2009; Armitage, 2011; Turner et al.,
2014). While purely hydrodynamic instabilities, most prominently the vertical
shear instability (VSI: Urpin & Brandenburg 1998; Urpin 2003; Nelson et al.
2013; Barker & Latter 2015; Lin & Youdin 2015), the convective overstability
(COS: Klahr & Hubbard 2014; Lyra 2014; Latter 2016) and its nonlinear
counterpart — baroclinic vortex amplification (Petersen et al., 2007a, b;
Lesur & Papaloizou, 2010; Lyra & Klahr, 2011; Raettig et al., 2013) are likely
too weak to explain global mass accretion rates in PPDs, they can strongly
affect the concentration of dust grains, and can thus be active elements in
the planetesimal formation process (Lesur et al., 2023)

In this paper, we focus on the COS. Several studies have investigated this
instability in terms of its linear growth phase and its nonlinear saturation,
which includes vortex formation (Lyra, 2014; Klahr & Hubbard, 2014; Latter,
2016; Teed & Latter, 2021). In addition, the concentration of dust in COS-
vortices was demonstrated in local shearing box simulations of gas including
Lagrangian dust particles (Lyra et al., 2018; Raettig et al., 2021). The
results of Raettig et al. (2021) suggest that vortices induced by the COS are
suitable sites for planetesimal formation, as their simulations revealed dust
densities exceeding the Roche density, depending on the disk radius under
consideration.

However, the aforementioned studies are based on a fully local shearing box
framework, where the radial buoyancy that powers the COS is treated as a
constant forcing. This means that the possible effects of the ensuing
instability on the disk structure are being ignored. The exception is the
simulations of Klahr & Hubbard (2014), which were conducted in a global
cylindrical setup, and more recently the simulations of Klahr et al. (2023),
which are in principle fully global. However, the former simulations merely
probed the linear growth of the instability, and, more importantly, they were
2D-axisymmetric, rather than 3D. The latter simulations were axisymmetric and
adopted small radial and vertical extents in favor of resolution.

The goal of this paper is to study the COS in a radially global disk setup, by
conducting both linear analyses and nonlinear hydrodynamic simulations (2D
axisymmetric and 3D). In contrast to local shearing box simulations, global
simulations self-consistently consider the disk’s background structure, and
possible changes thereof in response to the nonlinear evolution of
instabilities and their entailing angular momentum transport. This is expected
to be particularly relevant in 3D, where angular momentum can be transport
radially due to non-axisymmetric structures. Also the concentration of dust
grains should be affected by the global disk evolution. Furthermore, the most
unstable linear COS modes found in fully local calculations actually possess a
vanishing radial wavenumber (i.e. an infinite radial wavelength), which
somewhat contradicts the local approximation in the first place. Radially
global linear calculations presented below indeed confirm that the most
unstable modes have a radially global extent, but these modes are not captured
by a fully local model.

The paper is organised as follows. In Section 2 we introduce our hydrodynamic
model of a PPD. In Section 3 we summarize the theoretical background of the
COS, which will be useful in interpreting our results, and describe its
potential relevance to disk observations. In Section 4, we perform a radially
global linear analysis of the COS. In Section 5, we describe the setup and
results of our hydrodynamic simulations. In Section 6, we summarise and
discuss caveats of our work and prospects for future work.

##  2 Hydrodynamic Model

###  2.1 Basic equations

We consider a global hydrodynamic model of a PPD consisting of a compressible,
non-isothermal gas, governed by the set of fluid equations

$\left(\partial_{t}+\boldsymbol{v}\cdot\boldsymbol{\nabla}\right)\,\rho=-\rho\left(\boldsymbol{\nabla}\cdot\boldsymbol{v}\right),$
(1)
$\left(\partial_{t}+\boldsymbol{v}\cdot\boldsymbol{\nabla}\right)\,\boldsymbol{v}=-\frac{1}{\rho}\boldsymbol{\nabla}P+\frac{1}{\rho}\boldsymbol{\nabla}\cdot\boldsymbol{S}-\boldsymbol{\nabla}\Phi_{*},$
(2)
$\left(\partial_{t}+\boldsymbol{v}\cdot\boldsymbol{\nabla}\right)\,P=-\gamma
P\left(\boldsymbol{\nabla}\cdot\boldsymbol{v}\right)-\Lambda.$ (3)

In these equations $\rho$, $\boldsymbol{v}$ and $P$ are the gas volume mass
density, the three-dimensional velocity, and pressure, respectively, at time
$t$. We adopt a non-rotating frame with cylindrical coordinates
$(r,\varphi,z)$ and with origin on a central star of mass $M_{*}$ with
gravitational potential $\Phi_{*}$. In what follows, we neglect the vertical
variation of the potential and set $\Phi_{*}=-GM_{*}/r$, where $G$ is the
gravitational constant, which implies that we consider a region close to the
disk mid-plane such that $|z|\ll H$, where $H$ denotes the pressure scale
height defined below. Furthermore, the indirect gravitational term, self-
gravity, and magnetic fields are neglected. The remaining symbols in the above
equations are explained below.

We assume an ideal gas equation of state with pressure

$P=\frac{\mathcal{R}}{\mu}\rho T,$ (4)

where $\mathcal{R}$ is the gas constant and $\mu$ the mean molecular weight.
Following Klahr & Hubbard (2014), we assume gas cooling to occur in the
optically thin regime and define the cooling function

$\Lambda=\frac{\mathcal{R}}{\mu}\frac{\rho\delta T}{t_{c}},$ (5)

where $\delta T$ is the deviation of temperature from its equilibrium profile
specified below. We subsume viscous heating into this cooling function. In
what follows we work with the dimensionless constant cooling time

$\beta=\Omega_{\mathrm{K}}(r)t_{c}$ (6)

where

$\Omega_{\mathrm{K}}=\sqrt{\frac{GM_{*}}{r^{3}}}$ (7)

is the local Keplerian frequency. We fix the adiabatic index $\gamma=1.4$.

Finally,

$\boldsymbol{S}=\rho\nu\left[\boldsymbol{\nabla}\boldsymbol{v}+\left(\boldsymbol{\nabla}\boldsymbol{v}\right)^{\dagger}-\frac{2}{3}\boldsymbol{I}\boldsymbol{\nabla}\cdot\boldsymbol{v}\right]$
(8)

is the viscous stress tensor with the constant kinematic viscosity $\nu$. The
symbol $\dagger$ denotes the conjugate transpose and $\boldsymbol{I}$ stands
for the unit tensor. In our nonlinear simulations, viscous terms are only
included to ensure numerical stability, such that $\nu$ is chosen to be very
small. In particular, $\nu$ is much smaller than the typical turbulent
viscosities (defined in Section 5.1) measured in 3D simulations. Moreover,
below, we perform linear calculations, including a kinematic viscosity, to
describe the effect of external turbulence on the linear COS.

###  2.2 Disk equilibrium

We focus on the region close to the disk mid-plane $z=0$ and neglect vertical
gravity, so the equilibrium configuration does not depend on height $z$ away
from the mid-plane. In particular, vertical shear is absent, eliminating the
possible occurrence of the VSI. We assume an axisymmetric disk equilibrium
with radial power-law profiles for the gas volume mass density

$\rho(r)=\rho_{0}\left(\frac{r}{r_{0}}\right)^{-p}$ (9)

and temperature

$T(r)=T_{0}\left(\frac{r}{r_{0}}\right)^{-q}.$ (10)

with constant powerlaw indices $p$ and $q$, and reference density
$\rho_{0}=\rho(r_{0})$ and temperature $T_{0}=T(r_{0})$ at the reference
radius $r_{0}$.

Although the disk model is unstratified, we can define the characteristic disk
thickness or pressure scaleheight of a corresponding vertically stratified
model via

$H=c_{s}/\Omega_{\mathrm{K}},$ (11)

where

$c_{s}=\sqrt{\mathcal{R}T/\mu}$ (12)

is the sound speed.

The equilibrium velocity of the gas is given by

$\displaystyle\boldsymbol{v}(r)$
$\displaystyle=r\Omega(r)\boldsymbol{e}_{\varphi},$ (13)

with the orbital frequency

$\Omega(r)=\Omega_{\mathrm{K}}(r)\left[1-2\eta(r)\right]^{1/2},$ (14)

where we defined the dimensionless radial pressure gradient (Youdin & Goodman,
2005)

$\eta(r)\equiv-\frac{1}{2\rho(\Omega_{\mathrm{K}}r)^{2}}\frac{\partial
P}{\partial\ln r}=\frac{h^{2}}{2}\left(p+q\right)$ (15)

with the disk aspect ratio $h=H/r$. We also define $\Omega_{0}=\Omega(r_{0})$
and $H_{0}=H(r_{0})$. Furthermore, the squared radial buoyancy frequency is
given by

$N^{2}\equiv-\frac{1}{\gamma\rho}\frac{\partial P}{\partial r}\frac{\partial
S}{\partial
r}=-\frac{1}{\gamma}h^{2}\Omega^{2}\left(p+q\right)\left(q+\left[1-\gamma\right]p\right)$
(16)

with the dimensionless entropy

$\displaystyle S\equiv\ln{\frac{P}{\rho^{\gamma}}}.$ (17)

The second equality in Equations (15)-(16) assumes the radial power-law
profiles defined above.

##  3 Preliminaries on the COS

###  3.1 Theoretical aspects

The COS comprises an oscillatory destabilization of inertial waves of typical
frequency $\sim\Omega$, induced by a convectively unstable radial entropy
gradient, where epicyclic oscillations are amplified on account of a radial
buoyancy force which results from the cooling of gas occurring on a typical
time scale $t_{c}\sim 1/\Omega$ (see Latter 2016 for an illustrative
explanation and Teed & Latter 2021 for more details). Under typical conditions
where the radial pressure gradient $\partial_{r}P<0$, the COS requires a
negative radial entropy gradient, $\partial_{r}S<0$, since then the squared
radial buoyancy frequency (16), which is the central quantity for the linear
COS, is negative. The linear instability mechanism behind the COS was first
studied by Klahr & Hubbard (2014), Lyra (2014), and Latter (2016).

The nonlinear saturation of the COS was first explored in compressible 3D
shearing box simulations by Lyra (2014). These simulations revealed the
formation of large scale vortices and an entailing radial angular momentum
transport characterized by $\alpha$-viscous values on the order $\sim
10^{-3}.$ The nonlinear saturation of the axisymmetric COS was studied in
detail by Teed & Latter (2021) (citecite.teed2021TL21 henceforth) employing
local incompressible shearing box simulations. These authors found that two
dimensionless quantities largely determine the nonlinear outcome of the
instability. These are the Reynolds number

$Re=\frac{L^{2}\kappa}{\nu}\equiv\frac{1}{\alpha(r)},$ (18)

measuring the strength of viscous diffusion, and the (pseudo-)Richardson
number

$R=\frac{N^{2}}{\kappa^{2}}\equiv
h^{2}(r)\frac{\left(p+q\right)\left(q+\left[1-\gamma\right]p\right)}{\gamma},$
(19)

measuring the strength of the destabilizing radial entropy stratification
against the stabilizing effect of the radial angular momentum gradient. Here,
$L$ is a characteristic length scale of the flow and $\kappa$ denotes the
epicyclic frequency. Since our disk model is (as in citecite.teed2021TL21)
vertically unstratified, this length scale should be the vertical extent of
the simulation region. Our simulations below adopt a vertical domain
$L_{z}=0.5H$ (see §5), but for simplicity we set $L\sim H$, and we use
$\kappa\approx\Omega_{\mathrm{K}}$. Assuming this and considering the power-
law disks adopted above gives the second equalities in Eqs. (18)-(19). Note
that in contrast to citecite.teed2021TL21, our model is radially global, such
that $Re$ and $R$ are functions of disk radius. Specifically, $R$ ($Re$)
increases (decreases) with decreasing radius, implying a more vigorous COS at
smaller radii.

citecite.teed2021TL21 found that, depending on the values of $R$ and $Re$, the
COS saturates either into a weakly nonlinear state characterised by ordered
nonlinear waves, a state characterised by wave turbulence, or a state
supporting intermittent or persistent zonal flows and elevator flows (see
§5.3). They found that $R$ needs to exceed a critical value, roughly
delineated by the relation111The factor 0.1 instead of 0.05 (as in
citecite.teed2021TL21) arises since our scaling length is $H=2Lz$, rather than
$L_{z}$ .

$R_{\text{crit}}=0.1\left(Re/10^{4}\right)^{-1/2}$ (20)

in order for zonal flows to be formed in the saturated state. The point is
that the COS is only expected to play a role in dynamical process in PPDs if
zonal flows are formed, which for instance can concentrate dust grains.
Furthermore, in 3D simulations zonal flows are expected to wrap up into
vortices, which can induce radial angular momentum transport, and, which are
likewise able to concentrate dust (Raettig et al., 2021).

###  3.2 Occurrence in PPDs

The requirement of $N^{2}<0$ for the COS to operate in the disk mid-plane
typically translates to a rather steep temperature gradient and a flat density
profile (see Figure 1), which would be inconsistent with global disk profiles
inferred from early mm and sub-mm continuum observations of PPDs (e.g. Andrews
et al. 2009; Isella et al. 2009). However, flat density profiles could be
realized, for instance, near the edges of dead zones, planet-carved gaps, or
pressure bumps in general. As such, COS-active regions would be expected to
have typical radial widths $\Delta r\sim H$. On the other hand, more recent
ALMA imaging of HD 163296 (Mathews et al., 2013) and several disks in Lupus
(Tazzari et al., 2017) revealed very flat surface mass density profiles at
disk radii around tens of AU, which should in principle increase the
likelihood for the occurrence of the COS.

Whether the COS can operate also strongly depends on the gas’ cooling
timescale and is thus tied to the distribution of small dust grains tightly
coupled to the gas (Malygin et al., 2017; Barranco et al., 2018; Pfeil &
Klahr, 2019; Fukuhara et al., 2021). The above studies attempted to map out
the gas cooling time scale across the radial and vertical extent of a PPD.
However, such calculations require multiple simplifications, such that the
computed cooling times are subject to large uncertainties. Nevertheless, these
estimates suggest that one or more purely hydrodynamic instabilities are
likely active in planet-forming regions of the disk. In particular, the
results of Pfeil & Klahr (2019) suggest that the COS may occur in the
potential planet-forming region at $r=1-10\text{AU}$ if the structural
requirement $N^{2}<0$ is fulfilled.

Figure 1: Contour plot showing the minimal value of
$N^{2}/\Omega_{\mathrm{K}}(r_{0})^{2}$ within a radial region $r_{0}\pm
2H_{0}$ in a power law disk with density and temperature slopes $p$ and $q$,
respectively. The diagonal lines delineate the regions where
$\text{min}\left[N^{2}\right]$ takes negative or positive values, such that
the COS is in principle expected to occur in the parameter-space labeled
”COS”. Open (solid) circles represent values of $p$ and $q$ used in our
2D-simulations (2D as well as 3D simulations) below.

##  4 Radially global linear analysis

###  4.1 Linearised Equations

We now assume the disk equilibrium described in §2.2 and investigate the
stability of axisymmetric radially global, vertically local perturbations of
the form

$\delta\rho\to\delta\rho(r)\exp{\left(ik_{z}z-i\omega t\right)},$ (21)

and similar definitions for the perturbations $\delta v_{r}$, $\delta
v_{\varphi}$, $\delta v_{z}$ and $\delta P$. We define the complex mode
frequency $\omega=\omega_{R}+i\omega_{I}$, such that $\omega_{R}$ and
$\omega_{I}$ denote the real-valued mode’s frequency and growth rate,
respectively. Furthermore, $k_{z}$ is the real-valued vertical wavenumber and
we take $k_{z}>0$ without loss of generality. We also define

$\displaystyle W$ $\displaystyle=\frac{\delta P}{\rho},$ (22) $\displaystyle
Q$ $\displaystyle=c_{s}^{2}\frac{\delta\rho}{\rho},$ (23)

for convenience.

The linearised cylindrical equations then read

$\displaystyle i\omega\frac{Q}{c_{s}^{2}}=$
$\displaystyle\,\partial_{r}\ln\rho\delta v_{r}+ik_{z}\delta
v_{z}+\partial_{r}\delta v_{r}+\frac{\delta v_{r}}{r},$ (24) $\displaystyle
i\omega\delta v_{r}=$ $\displaystyle-2\Omega\delta
v_{\varphi}-\frac{\partial_{r}P}{\rho
c_{s}^{2}}Q+\partial_{r}W+\partial_{r}\ln\rho W,$ (25) $\displaystyle
i\omega\delta v_{\varphi}=$
$\displaystyle\left(2\Omega+r\partial_{r}\Omega\right)\delta v_{r},$ (26)
$\displaystyle i\omega\delta v_{z}=$ $\displaystyle\,ik_{z}W,$ (27)
$\displaystyle\begin{split}i\omega W=&\frac{\partial_{r}P}{\rho}\delta
v_{r}+\gamma c_{s}^{2}\left(ik_{z}\delta v_{z}+\partial_{r}\delta
v_{r}+\frac{\delta v_{r}}{r}\right)\\\
\quad&+\frac{1}{t_{c}}\left(W-Q\right).\end{split}$ (28)

Note that the driving force of the COS is the radial buoyancy force,
encapsulated in the term $\propto\partial_{r}P$ in the radial momentum
equation.

###  4.2 WKB dispersion relation

We can easily derive the WKB dispersion relation corresponding to Eqs.
(24)—(28) by ignoring all radial variations of the equilibrium disk structure,
as well as curvature terms. Furthermore, we introduce the complex radial WKB
wavenumber $k_{r}(r)$, assuming that $\partial_{r}\to ik_{r}(r)$. We then
obtain the fifth order dispersion relation

$-it_{c}\omega^{3}+\omega^{2}-c_{s}^{2}\left(k_{z}^{2}+\frac{\omega^{2}k_{r}^{2}}{\omega^{2}-\Omega^{2}}\right)\left(1-i\gamma\omega
t_{c}\right)=0.$ (29)

In the isothermal and adiabatic limits $t_{c}\to 0$ and $t_{c}\to\infty$ we
find the fourth order dispersion relation

$\omega^{2}-c_{s}^{2}\left(k_{z}^{2}+\frac{\omega^{2}k_{r}^{2}}{\omega^{2}-\Omega^{2}}\right)=0.$
(30)

and

$\omega^{2}-\gamma
c_{s}^{2}\left(k_{z}^{2}+\frac{\omega^{2}k_{r}^{2}}{\omega^{2}-\Omega^{2}}\right)=0,$
(31)

respectively. Note that (29) should not yield any growing solutions, which can
easily be shown in the isothermal and adiabatic limits, since in those cases
the dispersion relation yields a simple quadratic equation in $\omega^{2}$.
Furthermore, in the limit of high mode frequencies $|\omega|\gg\Omega$ we
isolate acoustic waves

$\omega^{2}-\gamma c_{s}^{2}k^{2}=0$ (32)

in the adiabatic limit, where $k=\sqrt{k_{r}^{2}+k_{z}^{2}}$. The isothermal
result can then be obtained by setting $\gamma$ to unity.

On the other hand, in the incompressible limit $c_{s}\to\infty$ we find

$\left(k_{z}^{2}+\frac{\omega^{2}k_{r}^{2}}{\omega^{2}-\Omega^{2}}\right)\left(1-i\gamma\omega
t_{c}\right)=0.$ (33)

The solution of the latter is the cooling mode

$\omega_{\text{COOL}}=-\frac{i}{\gamma t_{c}},$ (34)

decaying on cooling time scale $t_{c}$, as well as the pair of inertial waves

$\omega_{\text{IW}}=\pm\frac{k_{z}}{k}\Omega.$ (35)

Since the axisymmetric COS is captured within the incompressible Boussinesq
approximation (Latter 2016; Latter & Papaloizou 2017; citecite.teed2021TL21;
Lehmann & Lin 2023), we expect that short wavelength COS modes in our radially
global model essentially are inertial waves, described locally by Eq. (35).
This will indeed be confirmed below. The Boussinesq approximation can be
obtained from Eqs. (1)—(3) by assuming that the length scale of the phenomena
under investigation is much shorter than the radial and vertical disk
stratification lengths, that the flow is strongly subsonic
($|\delta\boldsymbol{v}|\ll c_{s0}$), and that the fractional perturbations of
density $\delta\rho/\rho\ll 1$ and pressure $\delta P/P\ll 1$, but with
$\delta\rho/\rho\gg\delta P/P$. The latter assumption is necessary to retain
the effect of the background entropy gradient in the model equations (see
Latter & Papaloizou 2017 for details).

###  4.3 Numerical methods

When supplemented with two boundary conditions accounting for the two first
order derivatives in $r$, such as

$\delta v_{z}(r_{\text{min}})=\delta v_{z}(r_{\text{max}})=0,$ (36)

Eqs. (24)—(28) constitute a linear eigenvalue problem of the form

$\boldsymbol{\mathcal{L}}\boldsymbol{b}=i\omega\boldsymbol{b}$ (37)

with the linear differential operator $\boldsymbol{\mathcal{L}}$, eigenvector
$\boldsymbol{b}=\left\\{Q,\delta v_{r},\delta v_{\varphi},\delta
v_{z},W\right\\}^{T}$ and complex eigenvalue $i\omega$. We use the spectral
code Dedalus (Burns et al., 2020) to solve the linear eigenvalue problem (37)
on a Chebyshev grid using a collocation point method (see, for instance, Lin
2021 for a similar treatment of a vertically stratified disk model).

Complementary to this, we verify the results obtained with the spectral code
using an alternative algorithm based on the shooting method. To this end, we
first recast the above equations as two first-order ordinary differential
equations in $r$:

$\begin{split}\partial_{r}\delta
v_{r}=&-\left(\frac{1}{r}+\frac{i\,\partial_{r}\ln\rho}{i+\gamma\omega
t_{c}}+\frac{\partial_{r}P\omega t_{c}}{P\left(i+\gamma\omega
t_{c}\right)}\right)\delta v_{r}\\\
\quad&+\left(-\frac{ik_{z}^{2}}{\omega}+\frac{\omega\rho\left(-1+i\omega
t_{c}\right)}{P\left(i+\gamma\omega t_{c}\right)}\right)W,\end{split}$ (38)
$\begin{split}\partial_{r}W=&i\bigg{(}\omega-\frac{2\Omega\left(2\Omega+r\partial_{r}\Omega\right)}{\omega}+\\\
\quad&\partial_{r}\ln P\frac{\left(\partial_{r}P-\gamma
P\partial_{r}\ln\rho\right)t_{c}}{\rho\left(i+\gamma\omega
t_{c}\right)}\bigg{)}\delta v_{r}\\\ \quad&+\left(\partial_{r}\ln
P\frac{\left(i+\omega t_{c}\right)}{\left(i+\gamma\omega
t_{c}\right)}-\partial_{r}\ln\rho\right)W,\end{split}$ (39)

which are then integrated using a fourth or fifth-order Runge-Kutta scheme
with pre-specified outer (we shoot from the outer to the inner domain
boundary) boundary values for $\delta v_{r}$ and $W$. The solution with the
correct value for $\omega$ is then obtained via a Newton-Raphson iteration
algorithm to match the corresponding inner two boundary values for one of the
fields (we use $\delta v_{z}$), obtained first with the spectral method. We
find that both methods generally are in excellent agreement. However, in
certain cases modes possess very short radial wavelengths ($k_{r}H_{0}\gtrsim
100$), such that the spectral method becomes prohibitively slow to obtain
fully resolved radial profiles of the eigenfunctions. Nevertheless, the method
can still find the correct eigenvalue $\omega$ corresponding to such
solutions, which we verified using example calculations with subsequently
increased radial resolution. For such cases, we present radial profiles
obtained with the shooting method, which can be computed efficiently even for
a large number of radial grid points $n_{r}\gtrsim 10^{4}$.

###  4.4 Results

####  4.4.1 Smooth disk region

Figure 2: Growth rates $\omega_{I}$ and frequencies $\omega_{R}$ as functions
of vertical wavenumber of linear, radially global COS modes in the fiducial
disk ($p=1.5$, $q=2$, $5\leq r/H_{0}\leq 15$).  Figure 3: Radial
eigenfunctions corresponding to the fastest growing COS modes in the fiducial
disk with $p=1.5$ and $q=2$ for increasing vertical wavenumbers
$k_{z}H_{0}=0.2,2,20$, as obtained with the spectral method. Notably, the
eigenfunctions become increasingly concentrated toward the inner domain
boundary with increasing $k_{z}$. The bottom panel shows the wavelet power
corresponding to the radial profile of $\operatorname{Re}\left[\delta
v_{z}\right]$. The red dashed curves in the wavelet power spectra correspond
to $k_{r}$ obtained from (35), with accordingly specified $k_{z}$ and using
$\omega_{\text{IW}}=\omega_{R}$. For the modes with $k_{z}H_{0}=2,20$ the
Lindblad resonance defined via (41) is located within the computational
domain, where the red dashed curves drop to zero.

Figure 2 shows frequencies $\omega_{R}$ and growth rates $\omega_{I}$ of
linear COS modes for vertical wavenumbers $0.2\leq k_{z}H_{0}\leq 1,000$,
obtained with the spectral code. Here we consider a radial domain $5\leq
r/H_{0}\leq 15$. As outlined in §1, typical regions that are unstable to the
COS are expected to be significantly smaller. However, in our simulations
below we will adopt the same radial domain, spanning $10H_{0}$. The reason is
that we wish to avoid spurious influences from the radial boundary conditions
on the long term nonlinear evolution in our simulations. Moreover, the linear
analysis of a smooth disk region does not depend qualitatively on the radial
domain size (cf. §4.4.2).

For each value of $k_{z}$ we generally find a band of modes within a range of
oscillation frequencies and growth rates, with the average frequencies and
growth rates increasing with increasing $k_{z}$. Moreover, larger growth rates
generally correspond to larger oscillation frequencies. Overall, oscillation
frequencies ranging from small values $\omega_{R}\sim 10^{-2}\Omega_{0}$ to
roughly the largest value of the Keplerian frequency within the considered
domain $\omega_{R}\approx 2.7\Omega_{0}$ (i.e. typical frequencies of local
inertial waves with $k_{z}\lesssim k$) can occur. Note that for smaller
$k_{z}$ the lowest occurring frequency is set by the radial grid resolution of
the spectral solver. Furthermore, modes with lower frequency generally extend
further out through the domain, whereas with increasing frequency (and hence
increasing growth rate) modes tend to be increasingly concentrated toward the
inner domain boundary. This is illustrated in Figure 3, which shows the radial
profiles of the various hydrodynamic quantities involved in the problem.

The behavior of the modes just described can be better understood if we
consider the WKB dispersion relation in the isothermal or adiabatic limit
(30)-(31). The latter can be rearranged to yield

$k_{r}^{2}=\frac{1}{\gamma
c_{s}^{2}}\left(\omega^{2}-\Omega^{2}\right)\left(1-\frac{\gamma
c_{s}^{2}k_{z}^{2}}{\omega^{2}}\right).$ (40)

This equation is analogous to Eq. (60) of Lubow & Pringle (1993) and Eq. (12)
of Svanberg et al. (2022), the latter applying to vertically stratified,
radially local disk models.

According to Eq. (40) neutral222with real-valued $k_{r}$ and $\omega$ waves
can only propagate in regions where the two conditions $\omega>(<)\Omega$
_and_ $\omega>(<)\sqrt{\gamma}c_{s}k_{z}$ are simultaneously fulfilled. The
first (second) equality applies to acoustic (inertial) waves. Since our linear
analysis above reveals no modes with $\omega>\Omega$, we conclude that COS
modes only exist radially inwards from their Lindblad resonances (LRs), whose
radii $r_{\text{LR}}$ are given by the condition

$\omega=\Omega(r_{\text{LR}}).$ (41)

This appears plausible, since in the Boussinesq approximation the linear COS
is known to excite inertial waves (Klahr & Hubbard, 2014; Lyra, 2014; Latter,
2016), and presumably there exist no linear (global) modes in our disk model
which connect both inertial and acoustic waves.

Nevertheless, modes that can appear in our simulations are restricted to
vertical wavenumbers $k_{z}H_{0}\gtrsim 10$, such that the modes in our
simulations can only fulfill the condition $\omega<\Omega$ and waves can only
propagate _inside_ of their Lindblad resonance radius. Note that generally,
our vertically unstratified analysis is only meaningful for modes with
$k_{z}H_{0}>1$.

The above implies that lower frequency modes have their LR at larger radii,
such that these modes possess smaller radial wavelengths within the
computational domain. Consequently, their growth rates are smaller, as their
coupling to the background entropy gradient becomes less efficient.

In addition, we compute a wavelet transform333We use the IDL ’wavelet’
procedure with a mother wavelet of type ’Paul’ (Torrence & Compo, 1998). of
the vertical velocity $\delta v_{z}$, which enables us to estimate the radial
wavenumber $k_{r}(r)$ of the eigenmodes. The bottom panel in Figure 3 shows
the wavelet power spectrum as a function of radius and radial wavenumber
$k_{r}H_{0}$. The mode with $k_{z}H_{0}=0.2$ shows a radial change of its
radial wavenumber that is characteristic of global waves, such as spiral
density waves in planetary rings (Shu, 1984). This radial change results from
the radially varying $\Omega_{\mathrm{K}}$. Note that the eigenfunctions for
the cases $k_{z}H_{0}=2$ and $k_{z}H_{0}=20$ are actually not wavelike, i.e.
they do not undergo at least one full oscillation in radial direction.
Therefore, the wavelet transforms for these cases cannot be expected to yield
a clear radial wavenumber. To improve the results we re-scaled the wavelet
power such that its minimal value is set to zero (by manual subtraction) for
these cases, as this procedure eliminates artificial residuals in the wavelet
power spectrum. The red dashed curve corresponds to the radial wavelength of
an inertial wave (35) with given vertical wavenumber $k_{z}$ and frequency
$\omega_{R}$.

Finally, we note that we checked the robustness of the computed modes by using
a number of different boundary conditions. Also we verified the existence of
some (including the fastest growing modes) with the shooting method.

####  4.4.2 Connection to local COS modes

Figure 4: Linear growth rates of radially global COS modes (solid curves) as
function of radial wavenumber $k_{r}$ for fixed $k_{z}H_{0}=10$ and for
different radial sizes of the computational domain $\Delta r$ are compared
with fully local growth rates (circles), computed via Eq. (42). The radial
wavenumbers $k_{r}$ are obtained from (29) and correspond to radial averages.
See the text for details. The insert figures illustrate contours of the
vertical velocity perturbation of the fastest growing modes in the cases
$\Delta r=5H_{0}$ and $\Delta r=H_{0}$.

We can attempt to compare the linear growth rates of the global COS modes
derived in the previous section with those of fully local modes. To this end,
we consider the expression for the growth rates of local modes with wavenumber
$\boldsymbol{k}=k_{r}\boldsymbol{e}_{r}+k_{z}\boldsymbol{e}_{z}$ (Klahr &
Hubbard, 2014; Lyra, 2014; Lehmann & Lin, 2023):

$\omega_{I,\text{loc}}=-\frac{\beta\mu^{2}N^{2}\Omega}{2\left(\Omega^{2}+\beta^{2}\mu^{2}\left[\Omega^{2}+N^{2}\right]\right)}$
(42)

where $\mu=k_{z}/\sqrt{k_{z}^{2}+k_{r}^{2}}$. This expression has been derived
using the incompressible Boussinesq approximation (Latter & Papaloizou, 2017),
which yields almost identical results to a local compressible calculation (see
Appendix A). We expect that the global COS modes derived in the previous
section increasingly resemble local COS modes as we decrease the radial domain
size, since then the variation of background quantities becomes smaller. This
is illustrated in Figure 4, where we present growth rates of global modes with
$k_{z}H_{0}=10$ for three different domain sizes $\Delta r=H_{0}$, $\Delta
r=5H_{0}$ and $\Delta r=10H_{0}$, as a function of radial wavenumber $k_{r}$,
represented by the solid lines. These are compared with local modes following
from (42), represented by open circles. The values of $k_{r}$ for each mode
are obtained using the WKB dispersion relation (29), where we assume values
for $c_{s}^{2},t_{c}$ and $\Omega$ corresponding to radial averages over the
entire domain. Similarly, radial averages of $N^{2}$ are used in (42). Note
that generally the WKB wavenumber $k_{r}$ is complex, and the values in Figure
4 are the corresponding real parts.

Thus, for domain sizes $\Delta r\lesssim H_{0}$, short wavelength modes with
$k_{r}H_{0}\gg 1$ resemble local COS modes, as indicated by the good match in
the corresponding growth rates. However, the fastest growing modes are still
global, as they do not extent throughout the entire domain unlike local modes,
but are concentrated toward the inner domain boundary. This is illustrated in
the insert figure, which shows the corresponding vertical velocity field for
$\Delta r=H_{0}$ and $\Delta r=5H_{0}$.

Figure 5: Maximal linear growth rates of radially global COS modes in a smooth
disk. Compared is the effect of a varying $\alpha$-viscosity (left panel),
cooling time (second panel), temperature slope (third panel) and density slope
(rightmost panel). If not indicated otherwise parameters correspond to the
fiducial disk ($q=2$, $p=1.5$, $\beta=1$, $\alpha=0$). Open circles in the
second panel from the left correspond to growth rates obtained from a fully
local model.

The reason is that the fastest growing modes’ LR lies within the computational
domain, implying that their frequencies $\omega_{R}\gtrsim\Omega_{0}$. Again,
we checked these modes’ robustness by using several different boundary
conditions. With increasing domain size the radial variation of the background
quantities (particularly $\Omega_{K}$) results in an increasing discrepancy
between local and global growth rates, as $k_{r}$ exhibits an increasing
variation throughout the domain. Obtaining local growth rates becomes less
meaningful for larger domain sizes.

####  4.4.3 Dependence on disk parameters

Figure 6: Illustration of the linear COS induced by a density structure in an
otherwise linearly stable disk. The Left panel shows radial profiles of
density $\rho$ (and the effective density slope $p_{\text{eff}}$) and squared
buoyancy frequency $N^{2}$ of the smooth disk, and disks containing either a
density bump or a density step with amplitude $A=0.4$. The right panel shows
linear growth rates of COS modes for the latter two cases for different
vertical wavenumbers $k_{z}$.

For completeness, Figure 5 shows growth rates of the fastest-growing COS modes
as a function of various disk parameters. The left panel compares different
values of an included $\alpha$-viscous stress444The indicated values of
$\alpha$ correspond to a disk radius $r=r_{0}$. (see Appendix B). Here, we
also indicate the region that is expected to be resolved in our hydrodynamic
simulations below (§5). The lower boundary is set by the vertical size of the
simulation domain ($\Delta z=0.5H_{0}$), whereas the upper boundary is set by
the grid resolution, assuming that 10 grid points per wavelength are required
to resolve a mode. Thus, based on the linear results, we expect the COS in our
simulations to be extinguished for $\alpha\gtrsim 10^{-4}$, which we will find
in reasonable agreement with the simulation results presented below. The
behavior of the curves in the remaining panels largely agrees with the fully
local calculations presented in previous works. That is, the growth rates in
the three rightmost panels are well explained by Eq. (42) using Eq. (16) for
the buoyancy frequency (its radial average). For reference, the open circles
in the second panel represent fully local growth rates computed using the
method explained in §4.4.2 for the fiducial case. It should be noted again
that modes $k_{z}H_{0}\lesssim 1$ are likely to be affected by the neglect of
vertical disk stratification.

####  4.4.4 Disk region containing a radial density structure

Figure 7: Hydrodynamic quantities (their real part) associated with the
fastest growing COS mode with $k_{z}H_{0}=20$ in a disk containing a pressure
bump (43) with $A=0.4$, as obtained with the shooting method. The mode is
excited in the region $0.7\lesssim r/r_{0}\lesssim 1$ where $N^{2}<0$ (see
Figure 6) but penetrates into the stable region at smaller radii. The bottom
panel shows the wavelet power spectrum of $\delta v_{z}$ as a function of
radial wavenumber $k_{r}$. The red dashed curve is the radial wavenumber
resulting from (35). The Lindblad resonance location where $\omega_{R}=\Omega$
is the location where the modes’ radial wavenumber drops to zero.

It is interesting to verify if the COS can operate in a more “realistic” disk
model. By this we mean a model of a disk region with temperature and density
slopes compatible with observational constraints, but containing a ”special”
structure that locally gives rise to the necessary temperature and density
gradients to yield a negative value of $N^{2}$ across a certain radial domain.
To this end we consider an initially convectively stable disk with $p=2$ and
$q=0.75$ (cf. Figure 1) on which we superimpose either a density bump

$\rho(r)\to\rho(r)\times\left[1+A\exp\left(-\frac{\left(r-r_{0}\right)^{2}}{2H_{0}^{2}}\right)\right],$
(43)

or a density step

$\rho(r)\to\rho(r)\times\left[1+\frac{A}{2}\tanh{\left(\frac{r-r_{0}}{H_{0}}\right)^{2}}\right],$
(44)

to render the region unstable to the COS. This is illustrated in Figure 6
(left panel) for an amplitude $A=0.4$. The upper panel displays the effective
density slope $p_{\text{eff}}=-\partial\log\rho/\partial\log r$, whereas the
lower panel shows $N^{2}$ for the three cases of a smooth disk, and a disk
containing a density bump and a density step, respectively. A linear analysis
adopting the latter two background disk structures yields growing modes for
$k_{z}H_{0}\gtrsim 5$. For increasing $k_{z}$ bands of growing modes are found
which become increasingly narrow and concentrated around the largest Keplerian
frequencies within the unstable region. The fastest growing modes have their
LR close to the outer boundary of the unstable region. However, note that even
for large $k_{z}H_{0}\gg 1$ the fastest growing modes possess growth rates
which are substantially smaller than in the smooth disk case discussed above.

It is interesting to note that in the present case most555We find some modes
with very small growth rates $\omega_{I}\lesssim 5\cdot 10^{-5}$ that traverse
the entire domain.  of the growing linear modes in this problem resemble to
some extent, resonantly forced (spiral) density waves (e.g., Shu 1984) since
they similarly appear to be excited within a narrow, unstable region and
propagate away from it with radially decreasing wavelength. As in the case of
the smooth disk, here we also find only modes radially inward of their
Lindblad resonances. In Figure 7 we show the fastest growing mode with
$k_{z}H_{0}=20$.

##  5 Hydrodynamic Simulations

###  5.1 Numerical setup

We use FARGO3D666http://fargo.in2p3.fr (Benítez-Llambay & Masset, 2016;
Benítez-Llambay et al., 2019) to perform hydrodynamic simulations in which we
integrate Eqs. (1)—(3) forward in time. As our disk is vertically
unstratified, we apply periodic boundary conditions to all quantities in the
azimuthal and vertical directions. Radial boundary conditions are such that
equilibrium values of gas density and azimuthal velocity are extrapolated into
the ghost zones, while gas radial and vertical velocities are set to zero at
the boundaries. Simulations are carried out in cylindrical coordinates as
defined in §2. Computational units are such that $r_{0}=M_{*}=G=1$. We adopt
$h_{0}=H_{0}/r_{0}=0.1$ in all runs. For our fiducial setup the numerical grid
covers $0.5\leq r/r_{0}\leq 1.5$ and $-0.25H_{0}\leq z\leq+0.25H_{0}$, and in
3D simulations $0\leq\varphi\leq 2\pi$. The grid resolution is $n_{r}\times
n_{z}\times n_{\phi}=2000\times 100\times 628$. Thus, the radial and vertical
resolutions are $200/H_{0}$, whereas the azimuthal resolution $\sim 10/H_{0}$
is significantly lower for reasons of computational costs777The typical size
of a single 3D simulation is $\approx 600$Gb., and since we expect structures
to form in our simulations to be predominantly axisymmetric or elongated in
azimuthal direction. 3D simulations are carried out on a GPU cluster, which is
necessary due to the high spatial resolution. All 3D simulations cover 1,000
reference orbits.

###  5.2 Diagnostics

Following previous studies (e.g. Lehmann & Lin, 2022), we describe the radial
and vertical turbulent angular momentum transport in 3D simulations via the
dimensionless quantities

$\alpha_{r}(t)=\frac{\langle\rho v_{r}\delta v_{\varphi}\rangle_{r\varphi
z}}{\langle P\rangle_{r\varphi z}}$ (45)

and

$\alpha_{z}(t)=\frac{\langle\rho v_{z}\delta v_{\varphi}\rangle_{r\varphi
z}}{\langle P\rangle_{r\varphi z}},$ (46)

respectively, where the brackets denote averaging over spatial dimensions as
indicated in the subscript, and where $\delta v_{\varphi}$ is the azimuthal
velocity deviation from its ground state value.

Furthermore, we define the deviation $\Delta X$ of a quantity $X(r,\varphi,z)$
from its vertical-azimuthal average across a cylindrical annulus at radius $r$
via

$\Delta X\equiv
X-\frac{\int\limits_{0}^{2\pi}\int\limits_{z_{\text{min}}}^{z_{\text{max}}}Xr\mathrm{d}\varphi\mathrm{d}z}{\int\limits_{0}^{2\pi}\int\limits_{z_{\text{min}}}^{z_{\text{max}}}r\mathrm{d}\varphi\mathrm{d}z}.$
(47)

If not stated otherwise, our radial diagnostic domain is defined by
$r_{\text{min}}/r_{0}=0.8$ and $r_{\text{max}}/r_{0}=1.2$.

###  5.3 Results of 2D axisymmetric simulations

In this section we present the results of 2D axisymmetric (radial-vertical)
simulations of the nonlinear saturation of the COS. In §5.3.1 we consider
simulations of a smooth power-law disk with $N^{2}<0$ for all radii. Next, in
§5.3.2 we consider a disk containing a pressure bump which locally gives rise
to the COS in an otherwise stable disk. The different values of $p$ and $q$
adopted throughout our simulations are shown in Figure 1 for illustration.

####  5.3.1 Nonlinear saturation of the COS in a smooth disk

Figure 8: a): The left panels show snapshots of the vertical velocity in the
course of the axisymmetric nonlinear saturation of the COS, where brighter
(darker) colors correspond to larger (smaller) values of $\delta v_{z}$. Note
that the color scaling of the bottom right frame at $1000$ orbits is the same
as in the left panel in the bottom frame b). The right frame shows the power
spectral density of the radial velocity in vertical wavelength space, and
illustrates the increased vertical elongation of structures, eventually
saturating into elevator flows and zonal flows with almost vanishing vertical
structure. The curves represent rms radial and vertical velocities as
indicated. The circles and the dashed line represent $\alpha_{r}(t)$ and its
time-average, respectively. See the text for details b): Saturated state of
the COS after 1,000 orbital periods, exhibiting persistent elevator flows
(left panel) and zonal flows (middle panel). The dashed curves in the middle
panel illustrate the geostrophic balance of the displayed zonal flows as
explained in the text. Figure 9: rms vertical velocities in 2D axisymmetric
hydrodynamic simulations. Compared are different cooling times $\beta$ (upper
left), initial temperature slopes $q$ (lower left), initial density slopes $p$
(lower right) and values of the dimensionless viscosity $\alpha$ defined via
Eqs. (18) and (8) at $r=r_{0}$. The dashed line near $100$ orbits in the upper
left panel represents the maximum growth rate predicted by the linear radially
global analysis above. It matches well with the largest linear growth rate
seen in the simulation for the same disk parameters. Note that the COS in
simulations with $\alpha\gtrsim 10^{-4}$ is fully damped for at least 10,000
orbits (not shown).

Figure 8 illustrates the nonlinear saturation of the COS in our fiducial
axisymmetric simulation. The upper left frame shows the evolution of initially
small scale, quasi-local modes into radially global modes. The latter
eventually develop persistent flow structures. The upper right frame shows the
evolution of the power spectral density of the radial velocity field
(contours), averaged over the radial diagnostic domain $0.8\leq r/r_{0}\leq
1.2$, as a function of vertical wavelength $\lambda_{z}$. Generally, we find
that the prevalent vertical length scale of flow structures increases with
time. Note that in the saturated state the prevalent vertical wavelength of
the radial velocity, which should represent secondary parasitic modes (Latter
2016; citecite.teed2021TL21) is $L_{z}/2$. This is consistent with our
numerical convergence analysis, presented in Appendix C, which shows that
smaller vertical domain sizes than our fiducial one result in a weakened
nonlinear saturation. The over-plotted solid curves are the rms vertical
(green) and radial (cyan) velocities. As explained in citecite.teed2021TL21,
the radial (vertical) rms velocities represent the power in the primary COS
(secondary parasitic) modes. The onset of parasitic modes eventually controls
the nonlinear saturation amplitude of the COS.

Furthermore, the circles represent the turbulent radial angular momentum flux
$\alpha_{r}$ computed via (45), such that open (solid) circles denote negative
(positive) values. The dashed line denotes its _negative_ time-averaged (over
the last 400 orbits) value $\langle\alpha_{r}\rangle_{t}\approx-1.3\cdot
10^{-6}$, and is roughly in agreement with corresponding values reported by
citecite.teed2021TL21 (their Figure 17) for $R\approx 0.05$. Note that since
our scaling length is $H$, rather than $L_{z}$ (as in citecite.teed2021TL21),
our value is to be multiplied by a factor of 4 for comparison.

The bottom frames in Figure 8 show snapshots of the vertical, azimuthal and
radial velocity field components, respectively, in the saturated state of the
simulation. The left and middle panels show the presence of persistent
elevator flows and zonal flows, respectively, which constitute radial patches
of alternating deviations of the corresponding velocity components (see
citecite.teed2021TL21 for details on these structures). Zonal flows are
accompanied by significant undulations in the radial pressure gradient, the
latter being required to establish geostrophic balance in these flow
structures. This is verified by the over-plotted curves of $2\Omega\delta
v_{\varphi}$ and $1/\rho\,\partial_{r}\delta P$, which are practically
identical, thereby demonstrating geostrophic balance between the Coriolis
force and the radial pressure gradient. Note that here, pressure variations
$\delta P$ are caused mainly due to variations in the density, as temperature
variations are found to be subdominant.

In Figure 9 we plot rms vertical velocities illustrating the dependence of the
nonlinear saturation of the COS on various physical simulation parameters.
These are the cooling time $\beta$ (upper left panel), the input viscosity
$\alpha$ defined via Eq. (18) (upper right panel), the initial radial density
slope $p$ (lower left panel) and the initial radial temperature slope $q$
(lower right panel). In all panels of Figure 9 the dashed curve represents our
fiducial simulation. Note the different simulation times in the left and right
panels . The fiducial values of $\beta$ and $p$ are chosen so as to maximize
linear growth rates and the saturation level of the COS reached within $1,000$
orbits. The adopted viscosity $\alpha=10^{-7}$ is roughly the largest possible
value which does not notably affect the evolution of the instability. All
simulations are initialized with low amplitude white noise $\sim
10^{-3}c_{s0}$ in all velocity components.

For $\alpha=10^{-7}$ we have $R_{\text{crit}}=0.0032$ by Eq. (20). The
simulations with varying density slope and temperate slope which develop
COS888For the simulations with $p=-2$ and $p=5$ the buoyancy frequency $N^{2}$
vanishes identically. Therefore, these simulations merely serve as a
consistency check. yield values in the range $0.037\lesssim R\lesssim 0.085$
and $0.008\lesssim R\lesssim 0.085$ within the diagnostic domain,
respectively, such that in all of these simulations persistent zonal flows are
expected to develop, at least after sufficient simulation time (see §3). This
is found to be largely the case, except for simulations with large cooling
times $\beta\gtrsim 100$. We suspect that in the latter simulations initially
another double diffusive instability develops, making use of viscous diffusion
(see Appendix D for details). However, after roughly 2,000 orbits we find the
emergence of elevator flows, indicating that the double diffusive instability
eventually gets supplanted by the COS. In the simulations presented here the
former instability is weak due to the small values of $\alpha$ adopted. We
hypothesize that after sufficient integration time the COS eventually
saturates to similar levels in all simulations shown in the upper left panel.

We also performed additional simulations using wave damping boundary
conditions (not shown). In these simulations we find a very similar nonlinear
saturation, albeit slightly delayed (by 100-200 orbits) compared to the
simulations shown in Figure 9. However, we find that the COS is almost
entirely suppressed (for at least 10,000 orbits) in simulations using damping
boundary conditions with $p=4.125$ and $p=-1.125$, and also in simulations
with cooling times $\beta\gtrsim 10$ (with final rms vertical velocity $\sim
10^{-6}$ which keeps decreasing). This is an interesting result, as it
suggests that the COS developing in our simulations is indeed a radially
global, rather than a local instability.

Figure 10: a): Nonlinear saturation of the COS in the vicinity of a pressure
bump. The upper left frame shows the evolution of the rms vertical velocity.
The dashed line indicates the linear growth predicted by our linear analysis.
The insert figure displays initial and final vertical velocities (at $z=0$).
The upper right frame shows contours of the vertical and azimuthal velocity
perturbations. The over-plotted curves are the initial radial buoyancy and
pressure profiles. b): Contours of the vorticity (59) defined below of a 3D
simulation for which we used the axisymmeteric saturation depicted in the top
panels as initial state. The dashed curves represent the vortencity
$\mathcal{L}$, defined in Eq. (48).

####  5.3.2 Nonlinear saturation of the COS in the vicinity of a pressure bump

In §4.4.4 we conducted a linear analysis of a disk region containing density
structures which give rise to the COS via a local modification of $N^{2}(r)$
in an otherwise stable disk, corresponding to the left most circle in Figure
1. Figure 10 illustrates the nonlinear evolution of the COS in the vicinity of
a pressure bump (the same as in Figure 6). In this simulation, we seeded the
linear eigenfunction with $k_{z}H_{0}=20$ (Figure 7) with a small initial
amplitude.

A similar outcome of the nonlinear saturation is obtained by seeding low
amplitude white noise as done before. However, seeding the fastest growing
mode significantly speeds up the evolution. Moreover, we can test its growth
rate against the linear analysis presented above. The vertical domain for this
simulation is set to be $\Delta z=\pi/5H_{0}$, such that the seeded mode fits
exactly two times in the vertical domain. The black curve in the upper left
frame displays the rms vertical velocity, with a linear growth phase of about
1,000 orbits, and the subsequent nonlinear saturation. The dashed green line
represents the linear growth rate as predicted by our linear analysis, and
which is indicated also in Figure 7. The insert figure shows the initial and
final radial profile of the vertical velocity (at $z=0$). Thus, even though
only the region $0.7\lesssim r/r_{0}\lesssim 1$ is unstable to the COS,
excited waves penetrate significantly into the stable region at smaller radii.
This is further illustrated in the upper right frames, which show contour
plots of the vertical and azimuthal velocity perturbations after 30,000
orbits. The over-plotted curves are the buoyancy and pressure profiles $N^{2}$
and $P$, respectively.

###  5.4 Results of 3D simulations

We now shift our focus to our 3D simulations. Given the significant
computational resources needed, we present a concise qualitative summary
instead of a comprehensive parameter exploration.

Figure 11: Nonlinear evolution of the COS in the fiducial 3D simulation. The
left panel shows snapshots of the vertically averaged vertical velocity. In
the right panel the contour plot represents the time evolution of the
azimuthal power spectral density of the radial velocity (planar cut at $z=0$).
The curve labeled ”$\langle P(n)\rangle_{t}$” is the time-averaged normalized
function (49), corresponding to the left plot-axis. The remaining curves are
the time evolution of the rms vertical velocity, rms radial velocity and
$\alpha_{r}$, corresponding to the right plot-axis.

In particular, we perform 3D simulations of smooth disks (§5.4.2) where we
consider three different values of the cooling time $\beta=0.1$, $\beta=1$ and
$\beta=10$, as well as four different values of the initial density slope
$p=-0.25$, $p=0.5$, $p=1.5$ and $p=2.5$. The latter values are also shown in
Figure 1. In the simulations adopting different $p$ we slightly adjust the
temperature slope $q$, such that the radial profile of $N^{2}$ is the same in
all simulations. Thus, essentially $q\approx 2$ in these simulations. In
addition, we consider one semi-3D simulation in the vicinity of a pressure
bump, which we discuss first in §5.4.1

####  5.4.1 Semi-3D simulation of the COS in the vicinity of a pressure bump

Before we describe our main results below, we briefly discuss one particular
3D simulation illustrated in the bottom panels of Figure 10. The latter shows
the contours of the vorticity of a 3D simulation that adopted the axisymmetric
nonlinear saturation shown in the top panels of the exact figure as the
initial state. In particular, the zonal flow formed around $r\approx 0.9$
becomes unstable to the Rossby wave instability (Lovelace et al., 1999) and
spawns several small vortices. This is consistent with the occurrence of a
maximum in the vortencity $\mathcal{L}$ (over-plotted dashed curve in all
panels) corresponding to the zonal flow, where

$\mathcal{L}=\frac{\Sigma}{2\omega_{z}}\left(\frac{\Pi}{\Sigma^{\gamma}}\right)^{2/\gamma}$
(48)

with the vertically integrate pressure $\Pi$. The vortices eventually merge
and strengthen, and the vortensity maximum disappears. The vortices then
excite (weak) spiral density waves (cf §5.4.4) and migrate radially inward.
The corresponding angular momentum flux is minimal, though, taking values
$\alpha_{r}\sim 10^{-5}$ (not shown). Note that this simulation’s non-
axisymmetric outcome is likely somewhat flawed due to the artificial
’switching on’ of the azimuthal dimension only after 30,000 orbits. It is
unclear if the zonal flow would have attained the same strength if the
simulation had been 3D from the outset.

####  5.4.2 Nonlinear saturation of the COS in a smooth disk

Figure 11 illustrates the nonlinear saturation of the COS in our fiducial 3D
simulation ($p=1.5$, $q=2$). Compared to Figure 8 we show here the evolution
of the vertically averaged vertical velocity (left panels) in the
$r\varphi$-plane. The vertical structure of the flow is similar to our 2D
simulations (5.3), but no notable _persistent_ zonal flows and elevator flows
form, as these are rapidly turned into vortices.

Figure 12:  Radial transport of angular momentum (top left), rms vertical
velocities (bottom left), radial mass accretion rate (top right), and rms
vertical transport of angular momentum (bottom right) in 3D simulations with
varying cooling time $\beta$. The quantity $M_{d}$ denotes the initial mass
contained within the diagnostic domain in each simulation.

The right panel shows the normalized azimuthal power spectrum of the mid-plane
radial velocity field

$P(n)\sim\left|\,\,\int\limits_{0}^{2\pi}\int\limits_{r_{\text{min}}}^{r_{\text{max}}}\mathrm{d}\varphi\mathrm{d}re^{in\varphi}\,\delta
v_{r}(z=0)\right|^{2},$ (49)

where $n$ is a positive integer denoting the azimuthal mode number. The curves
describe the rms vertical (red) and radial (green) velocities, as well as
$\alpha_{r}(t)$ (blue). Most power is located in modes with azimuthal mode
numbers $n=3-5$, owing to the excitation of spiral waves by vortices, which
are responsible for the radial angular momentum transport, measured through
$\alpha_{r}$ (see §5.4.4). It is noteworthy that the magnitude of the rms
radial and vertical velocities is reversed as compared with our 2D
simulations. This holds true for all simulations we conducted in this study.
We speculate that this difference is related to the presence of large-scale
vortices in 3D simulations, in contrast to persistent elevator flows in 2D
simulations.

In Figures 12 and 13 we show the time evolution of $\alpha_{r}$ (upper left
panels), the rms vertical velocities (lower left panels), the radial mass
accretion rate

$\dot{M}\equiv\int\limits_{0}^{2\pi}\int\limits_{z_{\text{min}}}^{z_{\text{max}}}r\delta
v_{r}\rho\mathrm{d}\varphi\mathrm{d}z$ (50)

(upper right panels), and $\text{rms}(\alpha_{z})$ (lower left panels), for
different cooling times $\beta$, and density slopes $p$, respectively. Our
measured $\alpha_{r}$-values are in the range $10^{-4}-10^{-3}$, roughly
consistent with values obtained from local compressible shearing box
simulations (Lyra, 2014).

####  5.4.3 Mass and angular momentum transport

In this section, we study the radial transport of mass and angular momentum
resulting from the COS in 3D simulations in more detail. Note that the
negative (time-averaged) angular momentum flux induced by the nonlinear COS in
our axisymmetric simulations is minimal (cf. §5.3.1), such that the disk’s
radial density structure is not altered to any significant extent (not shown).
This is consistent with the results of citecite.teed2021TL21.

Figure 13: Same as Figure 12 but now for simulations with varying initial
density slope $p$.

Despite its small magnitude though, the axisymmetric radial angular momentum
flux is the driving element for the formation of zonal flows, and hence
vortices (see §6 of citecite.teed2021TL21 for details). As mentioned above, in
our 3D simulations we generally find a significant outward radial transport of
angular momentum. Interestingly, in most cases this goes along with an
_outward_ radial transport of mass. For the case $p=1.5$ this is true at least
for cooling times $\beta=0.1$ and $\beta=1$. Indeed, using Eq. (50) we find
time-averaged accretion rates $\langle\dot{M}\rangle_{t}=-8.3\cdot 10^{-10}$,
$1.6\cdot 10^{-9}$, $2.3\cdot 10^{-7}$, $1.5\cdot 10^{-8}$, for simulations
adopting the initial slopes $p=-0.25,0.5,1.5,2.5$, in units
$M_{*}\,\Omega_{K0}$. We also find $\langle\dot{M}\rangle_{t}=1.6\cdot
10^{-8}$, $-8.9\cdot 10^{-8}$, for simulations adopting the cooling times
$\beta=0.1,10$ (and $p=1.5$). The finding of a positive accretion rate (50) is
contrary to the case of a laminar viscous accretion disk with a _constant_
kinematic viscosity. Nevertheless, our results are consistent with
axisymmetric viscous accretion disk theory (Lynden-Bell & Pringle, 1974) as we
show in the following.

Using (50) we can define the effective accretion velocity, as obtained from
our simulations, through

$v_{r,\text{sim}}\equiv\frac{\dot{M}}{2\pi r\Sigma},$ (51)

with the surface mass density

$\Sigma=\frac{1}{2\pi}\int\limits_{0}^{2\pi}\int\limits_{z_{\text{min}}}^{z_{\text{max}}}\rho\mathrm{d}\varphi\mathrm{d}z.$
(52)

For an axisymmetric laminar, viscous Keplerian accretion disk we have (Lynden-
Bell & Pringle, 1974)

$v_{r,\text{axi}}\equiv-\frac{3}{\sqrt{r}\Sigma}\partial_{r}\left(\overline{\nu}\Sigma\sqrt{r}\right),$
(53)

where in the current situation (i.e. for an $\alpha$-viscous prescription of a
turbulent, non-axisymmetric accretion disk) we define the axisymmetric
kinematic shear viscosity via

$\overline{\nu}(r,t)\equiv\alpha_{r}(r,t)\,\frac{c_{s0}^{2}}{\Omega_{\mathrm{K}}}$
(54)

with (cf §5.2)

$\alpha_{r}(r,t)=\frac{\langle\rho v_{r}\,\delta v_{\varphi}\rangle_{\varphi
z}}{\langle P\rangle_{\varphi z}}.$ (55)

In particular, $\alpha_{r}(r,t)$ and thus $\overline{\nu}$ may vary with
radius (and time), which explains our results, as shown below. Note that the
angular frequency $\Omega\approx\Omega_{\mathrm{K}}$ and the sound speed
$c_{s}^{2}$ given through (12) and (10) do not change significantly during our
simulations. Using (51) and (53) we find (for an axisymmetric laminar, viscous
Keplerian accretion disk)

$\dot{M}\sim\left[m(t)-\frac{1}{2}\right]\left(\frac{r}{r_{0}}\right)^{-m(t)}$
(56)

with

$m(t)\equiv a(t)+p(t)+3q/2-3,$ (57)

where we assumed

$\alpha_{r}(r,t)\sim\left(\frac{r}{r_{0}}\right)^{-a(t)}$ (58)

for simplicity. Note that $a(t)$ and $p(t)$, and hence also $m(t)$ vary with
time in our simulations. In contrast, as mentioned above, $q$ remains
essentially unchanged. In actuality, the density profile substantially
deviates from a simple power law at later stages of our simulations (cf.
Figures 14 and 16 below), but for simplicity we ignore this here. In the
current approximation a steady state density profile (such that the accretion
rate is radially constant) strictly requires $m(t)=0$. Furthermore, the pre-
factor $[m(t)-\frac{1}{2}]$ determines the sign and overall magnitude of the
accretion rate.

Figure 14 shows time- azimuth- and height-averaged radial profiles of
$\alpha_{r}$ (top panel) and $\rho$ (bottom panel) for simulations adopting
different initial density slopes $p$. The dashed curves are power law profiles
that have been fitted to the solid curves. The corresponding averaged power
law slopes are $\overline{a(t)}=0.74,1.40,3.68,0.69$ and
$\overline{p(t)}=-0.68,-0.23,-0.32,0.68$ for simulations adopting the initial
slopes $p=-0.25,0.5,1.5,2.5$. Using these values we find
$\overline{m(t)}=0.41,1.29,3.36,1.46$. Thus, once the COS saturates, none of
our simulated disks is predicted to remain in steady state, which is indeed
what we find in our simulations. Note that a larger absolute value of $m$
corresponds to an increasingly radially varying accretion rate, which should
result in an increasingly large modification of the radial density profile.
Indeed, the above $m$-values are qualitatively compatible with the amount of
radial mass-redistribution in our simulations, which is shown in Figure 16
(right frames). The left frames show simulations with different cooling times.
As seen in the panels for $N^{2}$, large portions of the simulation region are
stabilized with respect to the COS at later stages in some of the simulations.
Moreover, almost all simulations form at least one pressure bump, such that
$\partial_{r}\langle P\rangle_{\varphi z}=0$.

The aforementioned values of $\overline{m(t)}$ are also in qualitative
agreement with the values of the time-averaged accretion rate
$\langle\dot{M}\rangle_{t}$ given above. Indeed, using the values of $m(t)$
just derived, only the simulation with $p=-0.25$ is predicted to have a
negative accretion rate, whereas the other simulations are predicted to yield
a positive accretion rate. Moreover, the largest value is predicted for the
case $p=1.5$, followed by the case $p=2.5$, etc. , which is in agreement with
the values obtained from our simulations. Similarly, for the simulation with
$p=1.5$ and the longer cooling time $\beta=10$ we find
$\overline{a(t)}=-2.03$, $\overline{p(t)}=0.93$ and $\overline{m(t)}=-1.10$,
which is also in qualitative agreement with the measured accretion rate, which
is negative.

Finally, Figure 15 shows the time evolution of the radial accretion velocity,
obtained from Eqs. (51) and (53). Based on the degree of agreement between the
two quantities, we again conclude that the mass transport can _qualitatively_
be described by an axisymmetric laminar viscous accretion disk model.

Figure 14: The top and bottom panels show radial profiles (solid curves) of
the time-averaged $\alpha_{r}(r,t)$ as defined in (55) and mass density $\rho$
(averaged over time, azimuth and height), respectively, for simulations
adopting different initial density slopes $p$. The dashed lines are
corresponding linear fits. Note that the plots are double logarithmic. Figure
15: Testing the applicability of axisymmetric viscous accretion disk theory to
our 3D simulations of the COS. We compare the two different expressions for
the radial accretion velocity (51) and (53), labeled ”sim” and ”axi”,
respectively. The former is rigorously derived from our simulations, while the
latter assumes a laminar axisymmetric viscous disk with kinematic viscosity
(54).  Figure 16: Effect of the nonlinear COS on the radial disk structure.
From top to bottom we show initial (dashed curves) and final (after 1,000
orbits, solid curves) profiles of the radial pressure gradient (the vertically
and azimuthally averaged pressure), the density and the squared radial
buoyancy frequency. The left (right) panels compare different cooling times
$\beta$ (initial density slopes $p$).

####  5.4.4 Vortex migration and spiral density waves

Figure 17: Illustration of the excitation of spiral density waves by strong
vortices induced by the COS in the simulation with $p=2.5$. a): Vertically
averaged vertical component of vorticity defined in (59), as well as the
quantity $\Delta\mathcal{R}_{r\varphi}$ defined in (60). b): Vertical
structure of the upper left vortex in the left frames. The arrows indicate
velocities in the $rz$-plane and the dashed curve denotes the pressure
deviation $\Delta P$, showing the pressure maximum associated with the vortex.

Angular momentum transport in our 3D simulation is presumably realized through
the propagation and (shock-)damping of spiral density waves, excited by
sufficiently strong vortices (see Paardekooper et al. 2010 for an explanation
for the excitation mechanism). This is illustrated in Figure 17, where two
strong vortices are present in the computational domain, both of which excite
spiral waves. The left panel shows the vertically averaged vertical component
of the vorticity

$\boldsymbol{\omega}=\nabla\times\delta\boldsymbol{v},$ (59)

where $\delta\boldsymbol{v}$ is the velocity on top of the equilibrium
velocity (13). Furthermore, we consider the quantity (cf. §5.2)

$\Delta\mathcal{R}_{r\varphi}\equiv\Delta\left(\rho v_{r}\right)\Delta
v_{\varphi},$ (60)

representing the Reynolds stress induced by non-axisymmetric structures. This
quantity is displayed in the middle panel of Figure 17, showing that radial
angular momentum transport is entirely due to vortices. Moreover, the plot in
the right panel shows a vertical cut through the upper left vortex,
illustrating that strong COS-vortices can possess a significant vertical flow
structure.

The positive radial flux of angular momentum induced by the vortices leads to
radial inward migration of the latter (Paardekooper et al., 2010). This is
illustrated in Figure 18. Here, we follow Manger et al. (2020) and display the
quantity
$\text{min}[\langle\omega_{z}\rangle_{z}-\langle\omega_{z}\rangle_{z\varphi}]$,
where we additionally apply a smoothing filter to this quantity to suppress
the visibility of sharp features related to spiral waves, and to increase the
visibility of larger vortices. The dark diagonal streaks in the figure
represent vortices, which undergo rapid inward migration with typical speeds
of $0.01H$ per orbit. The faint horizontal bands visible in some panels are
spiral waves that survived our smoothing procedure. It appears that vortices
are more numerous in our simulations with larger $p\geq 1.5$. Also, in general
they more rapidly grow strong in the same simulations, which is reflected in
the earlier inward migration after their formation. This is consistent with
the values of $\alpha_{r}$ and $\dot{M}$ in figure 13. Furthermore, merging of
vortices can be seen in many cases. In addition, we find that the COS can, in
principle, generate large-scale, long-lived vortices that may be observable in
real PPDs based on their size. An example vortex with a radial width $\approx
3H$ is presented in Figure 19.

Finally, we note that the somewhat unexpected behavior of the simulation with
the longer cooling time $\beta=10$ (for which we expect the COS to develop
significantly later than in the case with $\beta=1$) could be related to the
fact that in the nonlinear state, sub-critical baroclinic vortex amplification
(Petersen et al., 2007a, b; Lesur & Papaloizou, 2010) should be active within
the vortices in our simulations, and the fact that the latter typically rotate
at a much slower rate than the local disk rotation rate, with a turnover
frequency $\sim\Omega_{\mathrm{K}}/\chi$, with the vortex aspect ratio $\chi$,
denoting the ratio of the vortex semi-major and semi-minor axes. Therefore,
the larger the vortex aspect ratio, the longer the rotation time, and as such
the longer the optimal cooling time for baroclinic amplification (Raettig et
al., 2013). Indeed, we find the occurrence of strong and large vortices in the
early stages of this simulation before 400 orbits (not shown), which produce
spiral density waves which undergo strong shocking. The gradual decay of
$\alpha_{r}$, $\dot{M}$ and $\alpha_{z}$ for the simulation with $\beta=10$
after 400 orbits is likely related to the strong modification of the radial
disk structure. Figure 16 shows that $N^{2}>0$ almost everywhere within the
diagnostic domain, such that the COS, and thus the production of new vortices,
should be strongly suppressed. Nevertheless, a detailed understanding of the
effect of the cooling time scale on COS vortices is beyond the scope of this
paper.

Figure 18: Illustration of the formation and subsequent migration of large
scale vortices in 3D simulations with different initial density slopes $p$.
The plots show contours of the quantity
$\text{min}[\langle\omega_{z}\rangle_{z}-\langle\omega_{z}\rangle_{z\varphi}]$
(we plot its square root for improved visibility). The dark small-scale
features, notably seen in the lower left panel, result from spiral wave shocks
which survived our applied smoothing algorithm.

####  5.4.5 Elliptic instability of vortices

Our simulation results suggest that vortices with a sufficiently small aspect
ratio $\chi$ are subjected to the elliptic instability (Lesur & Papaloizou,
2009), which should reduce their lifetimes.

Figure 19: Appearance of a large scale (radial width $\approx 3H$) long-lived
vortex with aspect ratio $\chi\approx 7$ in the 3D simulation adopting a short
cooling time $\beta=0.1$ (after $\approx 1,300$ orbits). This vortex exhibits
practically no vertical structure and vertical flow field (not shown), in
contrast to the vortex shown in Figure 17. The left panel shows the vertically
averaged $z$-component of vorticity (59). The arrows illustrate the retrograde
flow field associated with the vortex. The over-plotted dashed curve is the
scaled pressure deviation $\Delta P$ (at $z=0$), clearly showing a pressure
maximum associated with the vortex core. The right hand panel shows a polar
plot of the scaled density deviation $\Delta\rho$ . Spiral waves excited by
the vortex are clearly visible in the right panel.

To illustrate this, we compute the vertically averaged $z$-component of the
vorticity (59) and relate it to the vortex aspect ratio via (Lesur &
Papaloizou, 2010)

$\text{min}\left[\langle\omega_{z}\rangle_{z}\right]\equiv-C\frac{3}{2\chi}\left(\frac{\chi+1}{\chi-1}\right),$
(61)

which is valid for a Kida vortex (Kida, 1981) if the constant dimensionless
parameter $C=1$. Note that, strictly, this simplified correspondence assumes
that the minimum vorticity in our simulations occurs within the core of a
dominant (or representative) vortex at each time during the simulation. This
is a strongly idealized picture for two reasons. For one, several vortices are
present in the simulation domain during the nonlinear COS stage. Furthermore,
these vortices in general excite spiral waves. The latter form shocks in many
cases, which can dominate the minimum and maximum values of the vorticity. To
determine the proportionality factor $C$, we measure the aspect ratio of a
strong and resolved vortex in each of our simulations, using at least three
snapshots, and relate it to
$\text{min}\left[\langle\omega_{z}\rangle_{z}\right]$ via Eq. (61). To
suppress the influence of spiral waves, we apply the same smoothing procedure
as in Figure 18 above. We find $C=3.0,2.8,2.4,1.8$ for density slopes
$p=-0.25,0.5,1.5,2.5$, respectively.

For illustration Figure 20 shows the results for the cases $p=0.5$ and
$p=2.5$. In this figure, the solid curves represent $\chi$, whereas the dashed
curves show the quantity
$\text{max}\left[v_{z}\right]-\text{min}\left[v_{z}\right]$, which should be
an indicator for the onset of the elliptic instability, which is expected to
attack vortices with $\chi\lesssim 4$ (Lesur & Papaloizou, 2010). Indeed, in
many cases where the computed values of the _effective_ vortex aspect ratio
drop below this value in our simulations, we do find notable increases
(decreases) in the maximum (minimum vertical velocity). Note that the
correspondence mediated by Eq. (61) is in some cases compromised by the
presence of spiral wave shocks that survive the above-mentioned smoothing
procedure and which artificially reduce our computed $\chi$-values. The insert
figure in each panel shows an example of a vortex with $\chi<4$, which
contains fine-scale structures within its core, accompanied by strong vertical
velocity fluctuations, indicative of the elliptic instability.

##  6 Discussion

###  6.1 Summary

We studied the convective overstability near the mid-plane of radially global
protoplanetary discs. We first conducted a linear stability analysis.
Generally, we find that linear COS modes can only exist radially inward of
their corresponding Lindblad resonance. In particular, the fastest growing
modes are those whose Lindblad resonances lie within the computational domain,
such that their growth rate increases with decreasing resonance radius. Since
the modes are not resonantly excited, Lindblad resonances merely serve as
turning points in the present context.

The radially global, nonlinear saturation state of the COS in axisymmetric
simulations is similar to that in local incompressible shearing box
simulations (citecite.teed2021TL21). We find that the flow is dominated by
intermittent or persistent zonal flows and elevator flows, for sufficiently
large values of pseudo-Richardson number $R$. However, differences occur. In
§5.3.1 we found that for a viscosity $\alpha\gtrsim 10^{-5}$ the COS is
entirely suppressed in our simulations. For $Re=10^{5}$ we find
$R_{\text{crit}}\approx 0.1$ via Eq. (20), which is only marginally largr than
the largest values of $R$ (i.e. close to he inner radial domain boundary) for
our fiducial parameters. However, based on the results of
citecite.teed2021TL21 we would still expect the COS to generate some
hydrodynamic activity in form of a (disordered) nonlinear wave state in this
simulation. The reason for this discrepancy is unclear, but could be related
to the (inevitably) different boundary conditions applied in our simulations,
as compared to citecite.teed2021TL21. Nevertheless, the close resemblance of
the results for sufficiently large values of $R$ indicates that the nonlinear
saturation process of the axisymmetric COS in this parameter regime, which is
the most relevant one to planetesimal formation in PPDs, is well captured by a
radially local, incompressible model.

On the other hand, the radially global nature of the COS is not captured in
local shearing box models. This is expected to be relevant in more realistic,
non-smooth disc setups containing density features such as pressure bumps or
gaps. Here, we find, in agreement with our linear calculations, that the COS
may be locally excited but that its influence penetrates stable regions
radially inwards. Moreover, in simulations applying radial damping zones, we
find that the COS can be suppressed entirely (for at least 10,000 orbital
periods), which is not expected if it were a radially local instability.

The nonlinear evolution revealed by our 3D simulations appears to be rather
complex, and is - in contrast to axisymmetric simulations - no longer governed
only by the two quantities $Re$ and $R$ (§3). If the COS is sufficiently
strong, vortex formation ensues, which can result in significant radial mass
and angular momentum transport via the excitation of spiral density waves.
This transport modifies the radial disc structure, including forming pressure
bumps. Moreover, our results show that the (initial) radial disc structure
(described by the two slopes $p$ and $q$) affects the intensity of vortices
and spiral waves, and therefore also the magnitude of the ensuing radial mass
and angular momentum transport. This is another aspect that local shearing box
simulations cannot capture. Notably, the COS can induce _outward_ radial mass
transport, i.e., decretion, which is found to be the case in most of our 3D
simulations. In contrast, in simulations of VSI-active discs, one typically
finds an inward mass flux (Manger & Klahr, 2018; Lehmann & Lin, 2022). Note,
though, that such simulations are vertically stratified, unlike those
considered here.

Figure 20: Illustration of the elliptic instability in vortices resulting from
the COS. The solid curves show the effective aspect ratio $\chi$ of vortices
(61) occurring in 3D simulations with different initial density slopes $p$.
The dashed curves represent the maximum vertical velocity. The insert figures
show examples of vortices where the elliptic instability is suspected to be
active. In both cases, the displayed radial region is $0.7<r/r_{0}<1.1$ and
the azimuthal extent $\Delta\varphi=2$. Note that vortices around 300 orbits
in the simulation with $p=0.5$ are tiny, such that the elliptic instability
likely cannot be resolved, explaining the absence of notable increases of
$v_{z}$.

Our 3D simulations reveal the emergence of large-scale, long-lived vortices.
These generally migrate radially inward with rates $\sim 0.01H$ per orbit,
similar to vortices generated by the VSI (Manger et al., 2020; Lehmann & Lin,
2022). We also find that strong vortices with small aspect ratios are subject
to elliptic instability, which affects their lifetimes. Moreover, intense
vortices, such as the one in Figure 17, can contain significant vertical
structure in their flow field, whereas sufficiently weak vortices, such as the
one in Figure 19, lack any notable vertical structure. These differences are
expected to affect the concentration of dust by such vortices. We also find
that the COS can, in principle, spawn long-lived vortices (Figure 19) that
should be large enough to be observable by ALMA and could, therefore, be a
possible explanation for observed asymmetries in PPDs.

###  6.2 Implications for dust substructures

Small dust grains tightly coupled to the gas are expected to possess a
velocity

$\boldsymbol{v}_{\mathrm{d}}=\boldsymbol{v}+\tau_{\mathrm{s}}\frac{\nabla
P}{\rho_{\mathrm{tot}}},$ (62)

where $\tau_{\mathrm{s}}$ denotes the particle stopping time and
$\rho_{\mathrm{tot}}$ is total density of the dust-plus-gas mixture (Youdin &
Goodman, 2005). Most of our 3D simulations yield at least one radial location
where the radial pressure gradient vanishes (Figure 16). For sufficiently
small $\left|\boldsymbol{v}\right|$, dust grains respond to the local pressure
gradient and are thus expected to drift radially into these regions and, if
the pressure bump remains axisymmetric, produce dust rings. Dust rings may
also result from ‘traffic jams’ produced from differential dust drift owing to
a non-uniform radial pressure gradient, which is almost always the outcome of
the COS.

On the other hand, our 3D simulations also reveal that the nonlinear COS gives
rise to an averaged (turbulent) radial gas velocity (cf. Figure 15), which, as
explained above, in most cases results in mass decretion. Whether this
velocity is large enough to overcome pressure gradients and affect the dust
concentration in the pressure bumps described remains to be seen in future
work.

Finally, the COS readily produces vortices, which can trap dust radially and
azimuthally (Lyra & Lin, 2013; Raettig et al., 2015, 2021). The COS may,
therefore, explain dust rings and asymmetries observed in PPDs (Andrews et
al., 2018; Long et al., 2018; van der Marel et al., 2021). However, dust
trapping cannot be too significant; otherwise, planetesimal formation would
ensue, which is not observable.

###  6.3 Non-local particle stirring

In the more probable scenario where $N^{2}<0$ is only realized near particular
radii, our results indicate that the COS can still drive activity interior to
such regions due to its radially global character. This provides a background
turbulence to stir up dust grains. From the vertical Reynolds stress
$\text{rms}(\alpha_{z})$ (Figures 12-13), we can estimate the expected height
of a passive dust layer $H_{\mathrm{d}}$ to be (e.g. Lin 2021)

$\frac{H_{d}}{H}\approx\sqrt{\frac{\text{rms}(\alpha_{z})}{St}}\approx
10^{-3}St^{-1/2}.$ (63)

with the particle Stokes number $St=\tau_{\text{s}}\Omega_{\mathrm{K}}$, where
it is assumed that $St\gg\alpha_{z}$.

Considering a standard Minimum Mass Solar Nebula (Chiang & Youdin, 2010),
pebbles with an internal density of $1\,\text{g}\text{cm}^{-3}$ and a size of
$3\text{mm}$ would have $St\sim 10^{-3}$ at $1\text{AU}$, which are then
predicted to possess a scale height of $H_{\mathrm{d}}\approx 0.03H$. If the
vertically-integrated metallicity $Z=0.01$, then the resulting midplane dust-
to-gas ratio $\epsilon\sim 0.3$ would be insufficient for the streaming
instability to drive strong clumping and hence planetesimal formation should
not proceed (Johansen & Youdin, 2007). Furthermore, turbulent stirring in of
itself stabilizes the streaming instability (Umurhan et al., 2020; Chen & Lin,
2020). We can, therefore, expect a finite-sized pebble layer to be maintained
radially inward from where the COS is locally excited.

###  6.4 Caveats and outlook

In this study, we applied a simple optically thin gas cooling law with a
single parameter $\beta$. However, recent studies indicate that cooling times
favorable for the COS most likely imply that gas cooling is in the optically
thick regime close to the disk mid-plane. Future simulations should,
therefore, consider more realistic gas cooling prescriptions, incorporating
thermal diffusion. Furthermore, as our simulations are vertically
unstratified, they omit the possible presence of the VSI, which may occur for
the parameter regimes considered here, if vertical stratification and more
realistic cooling were accounted for. Fully global models are required to
study the possible interplay of the COS and the VSI in a realistic manner.

Our example simulation in which the COS results from the presence of a
pressure bump shows the formation of weak non-axisymmetries and vortices only
after ten thousands of orbits. It is, therefore, speculative whether the
required time scales are relevant to planetesimal formation in PPDs. It
further remains speculative if the resulting vortices can efficiently
concentrate dust. Also, more detailed observational constraints on the radial
density and temperature structure of PPDs are required to determine if the COS
is potentially relevant to planetesimal formation and whether or not it is
only expected to occur near certain density structures such as pressure bumps.

Due to limitations in computational resources, our results on the radial mass
transport in 3D simulations are not converged with respect to spatial
resolution (see Appendix C), and, perhaps more importantly, vertical domain
size (not shown). The latter implies that a quantitative study of the radial
mass transport due to the COS formally requires a vertically stratified disk
model. Such simulations will also need to adopt even higher grid resolutions
than in this work, which poses a great computational challenge.

We note that we do not fully understand the reasons behind the exceptionally
large radial mass flux occurring in our fiducial 3D simulation with $p=1.5$,
as compared to the other simulations. We find (not shown) that in the former
simulation the $\chi$-values computed via (61) seem to hover on the threshold
for marginal stability (i.e. $\chi\approx 4$), unlike the other simulations
(cf. Figure 20), where $\chi$ drops well below the threshold. Perhaps the
vortices in the simulation with $p=1.5$ are marginally unstable to the
elliptic instability, and might therefore survive longer as those in the other
simulations. Moreover, we find a relatively large pileup of mass in the
diagnostic domain in this simulation, compared to the other simulations (see
Figure 14). This explains to some extent the large mass flux. However, it
remains unclear how this behaviour is connected to the equilibrium density
slope $p$. We speculate that the radial vortencity profile $\rho/\omega_{z}$,
which yields a constant in the special case $p=1.5$ plays a role in the
evolution of vortices in our simulations. We therefore suggest that future
work should study the evolution of vortices in radially global disk models in
a more systematic manner.

Note also that our discussion on the elliptic instability in §5.4.5 is not
rigorous. As mentioned earlier, the correspondence mediated by Eq. (61) is in
some cases compromised by the presence of spiral wave shocks that survive our
smoothing procedure and which artificially reduce our computed $\chi$-values.
Moreover, the measurements of vortex aspect ratios are done by hand and are,
therefore, subject to significant potential inaccuracy. Furthermore, vortices
in our simulations are not necessarily Kida vortices, which are defined by a
constant patch of relative vorticity. Nevertheless, our results strongly
suggest that the elliptic instability occurs in the cores of COS-vortices of
sufficiently small aspect ratio.

In a follow-up study, we will investigate the impact of the nonlinear COS on
dust concentration in PPDs. In particular, we will clarify if COS vortices in
radially global PPD models can result in dust-to-gas density ratios large
enough to trigger vigorous streaming instability.

We thank an anonymous reviewer for useful comments. This work is supported by
the National Science and Technology Council (grants 112-2112-M-001-064-,
112-2124-M-002-003-, 113-2124-M-002-003-) and an Academia Sinica Career
Development Award (AS-CDA-110-M06). Simulations were performed on the _Kawas_
cluster at ASIAA and the _Taiwania-2_ cluster at the National Center for High-
performance Computing (NCHC). We thank NCHC for providing computational and
storage resources.

##  Appendix A COS in local compressible model

For completion we provide the linearised equations in a local compressible
model (cf. Lyra 2014). These are obtained from the radially global equations
(24)—(28) by replacing $\partial_{r}\to ik_{r}$ and using (9)-(10). The
resulting equations constitute a linear eigenvalue problem

$\displaystyle\boldsymbol{M}\boldsymbol{b}=\omega\boldsymbol{b},$ (A1)

with

$\boldsymbol{M}=\begin{bmatrix}-i\omega&cs^{2}\left(1+k_{r}-p\right)&0&-ik_{z}c_{s}^{2}&0\\\
\frac{2\eta}{c_{s}^{2}}&-i\omega&-2&0&ik_{r}-p\\\
0&2+r\partial_{r}\Omega&-i\omega&0&0\\\ 0&0&0&-i\omega&ik_{z}\\\
-\frac{1}{\beta}&-2\eta+c_{s}^{2}\gamma\left(1+ik_{r}\right)&0&ic_{s}^{2}\gamma
k_{z}&-i\omega+\frac{1}{\beta}\\\ \end{bmatrix},$

eigenvector $\boldsymbol{b}=\left\\{Q,\delta v_{x},\delta v_{y},\delta
v_{z},W\right\\}^{T}$ and eigenvalue $\omega$. The linear growth rates of
inertial waves resulting from this set of equations are largely similar to
those resulting from an incompressible Boussinesq model (Latter, 2016; Lehmann
& Lin, 2023) and will not be presented here for brevity.

##  Appendix B Linear problem including viscosity

The linearised radially global momentum equations including a _constant_
999this is chosen so as to facilitate comparison with our simulations.
kinematic shear viscosity $\nu$ are given by

$\begin{split}i\omega\delta v_{r}&=-2\Omega\delta
v_{\varphi}-\frac{\partial_{r}P}{\rho
c_{s}^{2}}Q+\partial_{r}W+\partial_{r}\ln\rho W\\\
\quad&\nu\big{[}\frac{4}{3r}\delta
v_{r}^{{}^{\prime}}+\frac{4}{3}\partial_{r}\delta
v_{r}^{{}^{\prime}}-(\frac{4}{3r^{2}}+k_{z}^{2})\delta
v_{r}+\frac{1}{3}ik_{z}\delta v_{z}^{{}^{\prime}}\\\
\quad&+\partial_{r}\ln\rho\left(\frac{4}{3}\delta
v_{r}^{{}^{\prime}}-\frac{2}{3r}\delta v_{r}-\frac{2}{3}ik_{z}\delta
v_{z}\right)\big{]}\end{split}$ (B1) $\begin{split}i\omega\delta
v_{\varphi}&=\left(2\Omega+r\partial_{r}\Omega\right)\delta v_{r}\\\
\quad&+\nu\big{[}\frac{1}{r}\delta
v_{\varphi}^{{}^{\prime}}+\partial_{r}\delta
v_{\varphi}^{{}^{\prime}}-(\frac{1}{r^{2}}+k_{z}^{2})\delta v_{\varphi}\\\
\quad&+\partial_{r}\ln\rho\left(\partial_{r}\delta
v_{\varphi}-\frac{1}{r}\delta v_{\varphi}\right)\\\
\quad&-\left(\partial_{r}\ln\rho\frac{\delta\rho}{\rho}-\frac{1}{\rho}\partial_{r}\delta\rho\right)r\partial_{r}\Omega\big{]},\end{split}$
(B2) $\begin{split}i\omega\delta
v_{z}&=ik_{z}W+\nu\big{[}-\frac{4}{3}k_{z}^{2}\delta v_{z}+\partial_{r}\delta
v_{z}^{{}^{\prime}}+\frac{1}{3}ik_{z}\delta v_{r}^{{}^{\prime}}\\\
\quad&+\frac{\delta v_{z}^{{}^{\prime}}}{r}+\frac{ik_{z}}{3r}\delta
v_{r}+\partial_{r}\ln\rho\left(\partial_{r}\delta v_{z}+ik_{z}\delta
v_{r}\right)\big{]},\end{split}$ (B3)

where a prime (′) indicates a radial derivative, and where we neglected the
equilibrium radial viscous inward drift:

$v_{r}=\frac{\nu}{2\Omega+r\partial_{r}\Omega}\left([3-p]\partial_{r}\Omega+r\partial_{r}^{2}\Omega\right),$
(B4)

which follows from the azimuthal component of Eq. (2) at equilibrium. In order
to solve the viscous problem (24), (28), (B1)—(B3) using the pseudo-spectral
method (§4.3) we now need to specify five additional boundary conditions. We
find that

$\displaystyle\delta v_{r}^{{}^{\prime}}(r_{\text{min}})$ $\displaystyle=0,$
(B5) $\displaystyle\delta v_{y}^{{}^{\prime}}(r_{\text{min}})$
$\displaystyle=\delta v_{y}^{{}^{\prime}}(r_{\text{max}})=0,$ (B6)
$\displaystyle\delta v_{z}^{{}^{\prime\prime}}(r_{\text{min}})$
$\displaystyle=0,$ (B7) $\displaystyle W(r_{\text{max}})$ $\displaystyle=0,$
(B8)

is a suitable set of boundary conditions.

On the other hand, to apply the shooting method we now consider the set of
first order ODEs:

$\begin{split}\partial_{r}U&=i\omega\delta v_{r}+2\Omega\delta
v_{\varphi}+\frac{\partial_{r}P}{\rho c_{s}^{2}}Q-\partial_{r}\ln\rho W\\\
\quad&-\left(\frac{1}{r}+\partial_{r}\ln\rho\right)\left(U-W\right)\\\
\quad&+\nu\bigg{[}\left(\frac{4}{3}\frac{1}{r^{2}}+\frac{2}{3}\frac{\partial_{r}\ln\rho}{r}+k_{z}^{2}\right)\delta
v_{r}-\frac{1}{3}ik_{z}\delta v_{z}^{{}^{\prime}}\\\
\quad&+\frac{2}{3}\partial_{r}\ln\rho ik_{z}\delta v_{z}\bigg{]},\end{split}$
(B9) $\begin{split}\partial_{r}\delta
v_{\varphi}^{{}^{\prime}}&=\frac{1}{\nu}\left[i\omega\delta
v_{\varphi}-\left(2\Omega+r\partial_{r}\Omega\right)\delta v_{r}\right]\\\
\quad&-\frac{Q}{c_{s}^{2}}\partial_{r}\Omega-\left(\frac{1}{r}+\partial_{r}\ln\rho\right)\delta
v_{\varphi}^{{}^{\prime}}\\\
\quad&+\left(\frac{1}{r^{2}}+k_{z}^{2}+\frac{\partial_{r}\ln\rho}{r}\right)\delta
v_{\varphi},\end{split}$ (B10) $\begin{split}\partial_{r}\delta
v_{z}^{{}^{\prime}}&=\frac{1}{\nu}\left[i\omega\delta
v_{z}-\frac{3}{4}ik_{z}W-\frac{1}{4}ik_{z}U\right]\\\
\quad&+\frac{4}{3}k_{z}^{2}\delta v_{z}-\frac{1}{r}\delta
v_{z}^{{}^{\prime}}-\frac{ik_{z}}{3r}\delta v_{r}\\\
\quad&-\partial_{r}\ln\rho\left(\delta v_{z}^{{}^{\prime}}+ik_{z}\delta
v_{r}\right),\end{split}$ (B11) $\displaystyle\partial_{r}\delta v_{r}$
$\displaystyle=\frac{3}{4\nu}\left(U-W\right),$ (B12)
$\displaystyle\partial_{r}\delta v_{\varphi}$
$\displaystyle=v_{\varphi}^{{}^{\prime}},$ (B13)
$\displaystyle\partial_{r}\delta v_{z}$ $\displaystyle=v_{z}^{{}^{\prime}},$
(B14)

where we defined the new variable $U$ via Eq. (B12), and where

$Q=\frac{A}{B}W+\frac{\gamma
c_{s}^{2}\partial_{r}\ln\rho-\frac{\partial_{r}P}{\rho}}{B}\delta v_{r},$
(B15)

following from (28), as well as

$\begin{split}W&=\left(1+\frac{4\nu i\omega
A}{3c_{s}^{2}B}\right)^{-1}\frac{4\nu}{3}\bigg{[}\bigg{(}-\frac{i\omega}{c_{s}^{2}}\frac{\gamma
c_{s}^{2}\partial_{r}\ln\rho-\frac{\partial_{r}P}{\rho}}{B}\\\
\quad&+(\partial_{r}\ln\rho+\frac{1}{r})\bigg{)}\delta
v_{r}-\left(ik_{z}\delta v_{z}+\frac{3}{4\nu}U\right)\bigg{]}\end{split}$
(B16)

following from (24). Furthermore, we defined the short notations

$\displaystyle A$ $\displaystyle\equiv i\omega-\frac{1}{t_{c}},$ (B17)
$\displaystyle B$ $\displaystyle\equiv i\omega\gamma-\frac{1}{t_{c}}.$ (B18)

##  Appendix C Numerical convergence of hydrodynamic simulations

###  C.1 2D axisymmetric simulations

Figure 21 shows results of 2D simulations with different values of numerical
parameters. These are the spatial resolution (upper panel), the radial domain
size (middle panel) and the vertical domain size (bottom panel). The plots in
the upper panel demonstrate that our 2D simulations, using a resolution of
$200/H$ are converged. The plots in the middle panel reveal that the
saturation level of the COS generally reduces with decreasing radial domain
size. Based on the results of our linear analysis in §4.4.1 this is expected,
and is due to the fact that the maximum value of $R$ defined in (19) is
reduced (the value at the inner domain boundary), as well as our finding that
for vertical wavenumbers $k_{z}H_{0}\gtrsim 10$ the fastest growing linear
modes (i.e. the modes that can ”access” the largest $R$-values) are strongly
concentrated towards the inner disk boundary. Thus, with decreasing radial
domain size these modes are gradually removed.

We also find a significant drop in the saturation level when decreasing the
vertical domain size as $\Delta z=0.5H_{0}\to 0.25H_{0}$, whereas the
corresponding drop is only modest for the reduction $\Delta z=H_{0}\to
0.5H_{0}$. These reductions imply that we increase the minimum vertical
wavenumbers $k_{z,\text{min}}=4\pi\to 8\pi$ and $k_{z,\text{min}}=2\pi\to
4\pi$, respectively. Considering Figure 2 (left panel) these findings are
plausible, since the largest growing modes in our simulations possess
$10\lesssim k_{z}H_{0}\lesssim 100$.

###  C.2 3D simulations

As pointed out in §6.1, the radial angular momentum flux in our 3D simulations
is not converged with respect to resolution in all spatial dimensions, which
is illustrated in Figure 22 (lower panel). This can likely be traced back to
the formation of spiral shocks (which generate very sharp gradients in the
hydrodynamic quantities), which are ultimately responsible for the deposition
of angular momentum in the disk. On the other hand, the plot in the upper
panel shows that rms radial (dashed curves) and vertical (solid curves)
velocities resulting from the COS are actually converged. Our chosen
resolution of 200/H in both radial and vertical directions proves
computationally demanding in 3D global simulations, making it currently
impractical to pursue an even higher resolution. Nevertheless, as already
noted in §6.1, a detailed study of the radial mass flux requires a vertically
stratified disk model.

Figure 23 illustrates the robustness of our (3D) simulation results with
respect to the radial boundary conditions. This is not a straight forward
issue in the present study, as the fastest growing COS modes are concentrated
towards the inner domain boundary. The plots show the time evolution of the
averaged density slope $\partial\ln\rho/\partial\ln r$ (upper panel) and the
$\alpha$-viscosity (lower panel). Compared are fiducial simulations with and
without wave damping radial boundary conditions. We find that the application
of a damping zone near the boundaries merely results in a delayed nonlinear
saturation of the COS, of equal magnitude as in the simulation with default
boundaries. The insert plots in the upper panel show the initial and final
radial profiles of the squared buoyancy frequency $N^{2}$ and pressure $P$, in
good agreement between the two simulations.

Figure 21: Testing convergence of 2D axisymmetric simulations with respect to
resolution (top panel), radial domain size (middle panel) and vertical domain
size (bottom panel). The simulations in the top panel were conducted using a
vertical domain size $\Delta z=H$. Note that all of these simulations adopted
wave damping boundaries. Figure 22: Testing convergence with respect to
spatial resolution in 3D simulations. We consider convergence simultaneously
in radial and vertical direction (labeled ”XZ”) as well as azimuthal direction
(labeled ”Y”). The upper panel shows rms vertical (solid curves) and radial
(dashed curves) velocities. The lower panel shows the time evolution of
$\alpha_{r}$. The curves in the lower panel have been smoothed for clarity.
The dashed curve in the lower panel describes the fiducial simulation. Figure
23: The plots illustrate the robustness of our simulation results with respect
to the radial boundary conditions.

##  Appendix D Double-diffusive instability in slowly cooled viscous disks

###  D.1 Occurrence in hydrodynamic simulations

In §5.3.1 we found that rms vertical velocities in 2D simulations do not show
the expected behavior of the COS with increasing cooling time. That is, for
large cooling times $\beta\gtrsim 100$ they exhibit an unexpected rise. We
suspect this to be caused by the presence of another instability. In the
following we briefly show that a ”double-diffusive” instability (not the COS)
can operate in the adiabatic limit and in the presence of an unstable entropy
gradient, as well as a considerable shear viscosity. It is double-diffusive in
that is relies on a large (viscous) diffusion of angular momentum and small
thermal diffusion.

We start by presenting in Figure 24 rms vertical velocities (right panel) of
adiabatic simulations ($\beta=10^{4}$) and constant viscosity taking values
$\alpha=10^{-5}-10^{-2}$. In these simulations we use $q=2$ and $p=0.5$, which
results in $N^{2}<0$, and a negligible equilibrium viscous accretion flow
according to Eq.(B4). We see that for $\alpha\gtrsim 10^{-4}$ the system
settles on appreciable rms velocity values. Furthermore, the left panels show
from top to bottom radial profiles of the vertically averaged radial velocity,
squared buoyancy frequency and entropy, respectively. As seen, the instability
gives rise to a radially varying radial outward velocity, such that (due to
the lack of cooling) entropy can be transported outwards, which results in a
flattening of the entropy profile. The buoyancy frequency develops strong,
stochastic deviations from the ground state profile. However, the averaged
value of $N^{2}$ (over height and radius) saturates to a finite negative value
(not shown), explaining the persistence of the instability in the nonlinear
regime. Note that the instability similarly develops for different values of
$p$. The difference is that the radial viscous accretion velocity (B4) adds up
to the radial velocity field induced by the instability. Depending on the sign
of the former (which depends on $p$) this somewhat speeds up or slows down the
flattening of the entropy profile.

###  D.2 Local incompressible description

The existence of such an instability was already predicted by Latter et al.
(2010), who studied a ”resistive double-diffusive” instability in magnetized
PPDs. The difference is that the version discussed here makes use of the
laminar Navier-Stokes viscosity instead of a magnetic field to diffuse angular
momentum of fluid parcels, in order to eliminate the stabilising effect of
rotation.

Let us consider the linearised equations for a gas within the local
incompressible Boussinesq model of Latter & Papaloizou (2017), as given in
Lehmann & Lin (2023) (their Eqs. (39)-(41)), for simplicity presented in the
limit101010It can be verified numerically that the linear instability persists
in this limit. $k_{x}\to 0$ and $\beta\to\infty$ (adiabatic limit), as well as
vanishing ground state velocity of the gas:

$\displaystyle\partial_{t}\delta\rho$
$\displaystyle=\frac{\rho_{0}}{H_{r}}\delta v_{x},$ (D1)
$\displaystyle\partial_{t}\delta v_{x}$ $\displaystyle=2\delta
v_{y}-H_{r}N^{2}\frac{1}{\rho_{0}}\delta\rho-k_{z}^{2}\alpha h_{0}^{2}\delta
v_{x},$ (D2) $\displaystyle\partial_{t}\delta v_{y}$
$\displaystyle=-\frac{1}{2}\delta v_{x}-k_{z}^{2}\alpha h_{0}^{2}\delta
v_{y},$ (D3)

where the radial entropy stratification is described by the quantity

$\frac{1}{H_{r}}\equiv\frac{1}{\gamma}\frac{\partial S}{\partial r}$ (D4)

and $N^{2}$ is given by (16). Furthermore, $\alpha$ is a dimensionless
viscosity parameter. The above equations are defined in a Cartesian frame
$(x,y,z)$ within small rectangular patch of the disk, centered on a fiducial
point $(r_{0},\phi_{0}-\Omega_{0}t,z=0)$, that rotates around the star at
approximately the local Keplerian frequency. The unit vectors
$\vec{\mathbf{e}}_{x},\,\vec{\mathbf{e}}_{y},\,\vec{\mathbf{e}}_{z}$ in the
box corresponding to the radial, azimuthal, and vertical directions in the
global disk.

In the _inviscid_ case $\alpha\to 0$, by taking the time derivative of Eq.
(D2) we find

$\partial_{t}^{2}\delta v_{x}=-\left(1+N^{2}\right)\delta v_{x}$ (D5)

which yields stable, buoyancy-adjusted epicycles. However, if we assume that
the viscosity is non-zero and consider an initial radial velocity perturbation
$\delta v_{x}>0$. Then, for sufficiently large $k_{z}$ the viscous term in the
azimuthal momentum equation prevents the azimuthal velocity $\delta v_{y}$
from becoming large enough to provide a restoring motion via the Coriolis term
in the radial momentum equation. At the same time, the buoyancy term is found
to be positive via (D1) and needs to exceed the viscous damping term in (D2)
in order for $\delta v_{x}$ to grow. The ensuing instability is thus expected
to be non-oscillatory, which can easily be confirmed as follows. The
dispersion relation for local modes $\propto\exp{\left(ik_{x}x+ik_{z}z-i\omega
t\right)}$ with $k_{x}=0$ following from above equations reads

$\omega^{3}+\omega^{2}\,2ik_{z}^{2}\alpha
h_{0}^{2}-\omega\left(1+N^{2}+k_{z}^{4}\alpha^{2}h_{0}^{4}\right)-ik_{z}^{2}\alpha
h_{0}^{2}N^{2}=0.$ (D6)

If we assume purely growing modes with $\omega\equiv i\omega_{I}$ we find that
all coefficients of the resulting dispersion relation are positive, such that
the condition

$k_{z}^{2}\alpha h_{0}^{2}N^{2}<0$ (D7)

is a sufficient condition for instability, expressing the necessity of a non-
vanishing viscosity _and_ an unstable entropy gradient with $N^{2}<0$.

Nevertheless, even though the linear instability can be understood within the
local Boussinesq framework, its nonlinear saturation involves considerable
changes of the background disk structure, as shown above. Furthermore, the
relevance of this instability to PPDs is questionable, since the considered
values of $\alpha$ would require a significant level of turbulence generated
by some other source in the first place.

Figure 24: Illustration of the nonlinear saturation of the double diffusive
instability powered by viscous momentum diffusion and a negative radial
entropy gradient. The right panel shows rms vertical velocities in simulations
with increasing kinematic shear viscosity $\alpha$, defined via (18). The left
panels show the evolution of the vertically averaged radial velocity, squared
buoyancy frequency and entropy, in the simulation with $\alpha=10^{-3}$.

## References

  * Andrews et al. (2009) Andrews, S. M., Wilner, D. J., Hughes, A. M., Qi, C., & Dullemond, C. P. 2009, ApJ, 700, 1502, doi: 10.1088/0004-637X/700/2/1502
  * Andrews et al. (2018) Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41, doi: 10.3847/2041-8213/aaf741
  * Armitage (2011) Armitage, P. J. 2011, ARA&A, 49, 195, doi: 10.1146/annurev-astro-081710-102521
  * Bai & Stone (2010) Bai, X.-N., & Stone, J. M. 2010, ApJ, 722, 1437, doi: 10.1088/0004-637X/722/2/1437
  * Balbus & Hawley (1991) Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214, doi: 10.1086/170270
  * Barker & Latter (2015) Barker, A. J., & Latter, H. N. 2015, MNRAS, 450, 21, doi: 10.1093/mnras/stv640
  * Barranco et al. (2018) Barranco, J. A., Pei, S., & Marcus, P. S. 2018, ApJ, 869, 127, doi: 10.3847/1538-4357/aaec80
  * Benítez-Llambay et al. (2019) Benítez-Llambay, P., Krapp, L., & Pessah, M. E. 2019, ApJS, 241, 25, doi: 10.3847/1538-4365/ab0a0e
  * Benítez-Llambay & Masset (2016) Benítez-Llambay, P., & Masset, F. S. 2016, The Astrophysical Journal Supplement Series, 223, 11.  http://stacks.iop.org/0067-0049/223/i=1/a=11
  * Blum (2018) Blum, J. 2018, Space Sci. Rev., 214, 52, doi: 10.1007/s11214-018-0486-5
  * Burns et al. (2020) Burns, K. J., Vasil, G. M., Oishi, J. S., Lecoanet, D., & Brown, B. P. 2020, Phys. Rev. Res., 2, 023068, doi: 10.1103/PhysRevResearch.2.023068
  * Carrera et al. (2021) Carrera, D., Simon, J. B., Li, R., Kretke, K. A., & Klahr, H. 2021, AJ, 161, 96, doi: 10.3847/1538-3881/abd4d9
  * Chen & Lin (2020) Chen, K., & Lin, M.-K. 2020, ApJ, 891, 132, doi: 10.3847/1538-4357/ab76ca
  * Chiang & Youdin (2010) Chiang, E., & Youdin, A. N. 2010, Annual Review of Earth and Planetary Sciences, 38, 493, doi: 10.1146/annurev-earth-040809-152513
  * Drazkowska et al. (2022) Drazkowska, J., Bitsch, B., Lambrechts, M., et al. 2022, arXiv e-prints, arXiv:2203.09759.  https://arxiv.org/abs/2203.09759
  * Fukuhara et al. (2021) Fukuhara, Y., Okuzumi, S., & Ono, T. 2021, ApJ, 914, 132, doi: 10.3847/1538-4357/abfe5c
  * Gammie (1996) Gammie, C. F. 1996, ApJ, 457, 355, doi: 10.1086/176735
  * Isella et al. (2009) Isella, A., Carpenter, J. M., & Sargent, A. I. 2009, ApJ, 701, 260, doi: 10.1088/0004-637X/701/1/260
  * Johansen & Youdin (2007) Johansen, A., & Youdin, A. 2007, ApJ, 662, 627, doi: 10.1086/516730
  * Kida (1981) Kida, S. 1981, Journal of the Physical Society of Japan, 50, 3517, doi: 10.1143/JPSJ.50.3517
  * Klahr et al. (2023) Klahr, H., Baehr, H., & Melon Fuksman, J. D. 2023, arXiv e-prints, arXiv:2305.08165, doi: 10.48550/arXiv.2305.08165
  * Klahr & Hubbard (2014) Klahr, H., & Hubbard, A. 2014, ApJ, 788, 21, doi: 10.1088/0004-637X/788/1/21
  * Latter (2016) Latter, H. N. 2016, MNRAS, 455, 2608, doi: 10.1093/mnras/stv2449
  * Latter et al. (2010) Latter, H. N., Bonart, J. F., & Balbus, S. A. 2010, MNRAS, 405, 1831, doi: 10.1111/j.1365-2966.2010.16556.x
  * Latter & Papaloizou (2017) Latter, H. N., & Papaloizou, J. 2017, MNRAS, 472, 1432, doi: 10.1093/mnras/stx2038
  * Lehmann & Lin (2022) Lehmann, M., & Lin, M. K. 2022, A&A, 658, A156, doi: 10.1051/0004-6361/202142378
  * Lehmann & Lin (2023) Lehmann, M., & Lin, M.-K. 2023, MNRAS, 522, 5892, doi: 10.1093/mnras/stad1349
  * Lesur & Papaloizou (2009) Lesur, G., & Papaloizou, J. C. B. 2009, A&A, 498, 1, doi: 10.1051/0004-6361/200811577
  * Lesur & Papaloizou (2010) Lesur, G., & Papaloizou, J. C. B. 2010, A&A, 513, A60, doi: 10.1051/0004-6361/200913594
  * Lesur et al. (2023) Lesur, G., Ercolano, B., Flock, M., et al. 2023, in Protostars and Planets VII, ed. S.-i. Inutsuka, T. Muto, K. Tomida, & M. Tamura, 465, doi: 10.48550/arXiv.2203.09821
  * Li & Youdin (2021) Li, R., & Youdin, A. N. 2021, ApJ, 919, 107, doi: 10.3847/1538-4357/ac0e9f
  * Lin (2021) Lin, M.-K. 2021, ApJ, 907, 64, doi: 10.3847/1538-4357/abcd9b
  * Lin & Youdin (2015) Lin, M.-K., & Youdin, A. N. 2015, ApJ, 811, 17, doi: 10.1088/0004-637X/811/1/17
  * Long et al. (2018) Long, F., Pinilla, P., Herczeg, G. J., et al. 2018, ApJ, 869, 17, doi: 10.3847/1538-4357/aae8e1
  * Lovelace et al. (1999) Lovelace, R. V. E., Li, H., Colgate, S. A., & Nelson, A. F. 1999, ApJ, 513, 805, doi: 10.1086/306900
  * Lubow & Pringle (1993) Lubow, S. H., & Pringle, J. E. 1993, ApJ, 409, 360, doi: 10.1086/172669
  * Lynden-Bell & Pringle (1974) Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603, doi: 10.1093/mnras/168.3.603
  * Lyra (2014) Lyra, W. 2014, ApJ, 789, 77, doi: 10.1088/0004-637X/789/1/77
  * Lyra & Klahr (2011) Lyra, W., & Klahr, H. 2011, A&A, 527, A138, doi: 10.1051/0004-6361/201015568
  * Lyra & Lin (2013) Lyra, W., & Lin, M.-K. 2013, ApJ, 775, 17, doi: 10.1088/0004-637X/775/1/17
  * Lyra et al. (2018) Lyra, W., Raettig, N., & Klahr, H. 2018, Research Notes of the American Astronomical Society, 2, 195, doi: 10.3847/2515-5172/aaeac9
  * Malygin et al. (2017) Malygin, M. G., Klahr, H., Semenov, D., Henning, T., & Dullemond, C. P. 2017, A&A, 605, A30, doi: 10.1051/0004-6361/201629933
  * Manger & Klahr (2018) Manger, N., & Klahr, H. 2018, MNRAS, 480, 2125, doi: 10.1093/mnras/sty1909
  * Manger et al. (2020) Manger, N., Klahr, H., Kley, W., & Flock, M. 2020, MNRAS, 499, 1841, doi: 10.1093/mnras/staa2943
  * Mathews et al. (2013) Mathews, G. S., Klaassen, P. D., Juhász, A., et al. 2013, A&A, 557, A132, doi: 10.1051/0004-6361/201321600
  * Nelson et al. (2013) Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610, doi: 10.1093/mnras/stt1475
  * Paardekooper et al. (2010) Paardekooper, S.-J., Lesur, G., & Papaloizou, J. C. B. 2010, ApJ, 725, 146, doi: 10.1088/0004-637X/725/1/146
  * Petersen et al. (2007a) Petersen, M. R., Julien, K., & Stewart, G. R. 2007a, ApJ, 658, 1236, doi: 10.1086/511513
  * Petersen et al. (2007b) Petersen, M. R., Stewart, G. R., & Julien, K. 2007b, ApJ, 658, 1252, doi: 10.1086/511523
  * Pfeil & Klahr (2019) Pfeil, T., & Klahr, H. 2019, ApJ, 871, 150, doi: 10.3847/1538-4357/aaf962
  * Raettig et al. (2015) Raettig, N., Klahr, H., & Lyra, W. 2015, ApJ, 804, 35, doi: 10.1088/0004-637X/804/1/35
  * Raettig et al. (2013) Raettig, N., Lyra, W., & Klahr, H. 2013, ApJ, 765, 115, doi: 10.1088/0004-637X/765/2/115
  * Raettig et al. (2021) Raettig, N., Lyra, W., & Klahr, H. 2021, ApJ, 913, 92, doi: 10.3847/1538-4357/abf739
  * Safronov (1972) Safronov, V. S. 1972, Evolution of the protoplanetary cloud and formation of the earth and planets.
  * Shu (1984) Shu, F. H. 1984, in Planetary Rings, 513–561
  * Simon et al. (2016) Simon, J. B., Armitage, P. J., Li, R., & Youdin, A. N. 2016, ApJ, 822, 55, doi: 10.3847/0004-637X/822/1/55
  * Svanberg et al. (2022) Svanberg, E., Cui, C., & Latter, H. N. 2022, MNRAS, 514, 4581, doi: 10.1093/mnras/stac1598
  * Tazzari et al. (2017) Tazzari, M., Testi, L., Natta, A., et al. 2017, A&A, 606, A88, doi: 10.1051/0004-6361/201730890
  * Teed & Latter (2021) Teed, R. J., & Latter, H. N. 2021, MNRAS, doi: 10.1093/mnras/stab2311
  * Torrence & Compo (1998) Torrence, C., & Compo, G. P. 1998, Bulletin of the American Meteorological Society, 79, 61
  * Turner & Drake (2009) Turner, N. J., & Drake, J. F. 2009, ApJ, 703, 2152, doi: 10.1088/0004-637X/703/2/2152
  * Turner et al. (2014) Turner, N. J., Fromang, S., Gammie, C., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 411, doi: 10.2458/azu_uapress_9780816531240-ch018
  * Umurhan et al. (2020) Umurhan, O. M., Estrada, P. R., & Cuzzi, J. N. 2020, ApJ, 895, 4, doi: 10.3847/1538-4357/ab899d
  * Urpin (2003) Urpin, V. 2003, A&A, 404, 397, doi: 10.1051/0004-6361:20030513
  * Urpin & Brandenburg (1998) Urpin, V., & Brandenburg, A. 1998, MNRAS, 294, 399, doi: 10.1046/j.1365-8711.1998.01118.x
  * van der Marel et al. (2021) van der Marel, N., Birnstiel, T., Garufi, A., et al. 2021, AJ, 161, 33, doi: 10.3847/1538-3881/abc3ba
  * Youdin & Johansen (2007) Youdin, A., & Johansen, A. 2007, ApJ, 662, 613, doi: 10.1086/516729
  * Youdin & Goodman (2005) Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459, doi: 10.1086/426895