# “Halfway to Rayleigh” and other Insights to the Rossby Wave Instability

Eonho Chang Graduate Interdisciplinary Program in Applied Mathematics,
University of Arizona, Tucson, AZ 85721, USA  Department of Astronomy and
Steward Observatory, University of Arizona, Tucson, AZ 85721, USA  Andrew N.
Youdin Department of Astronomy and Steward Observatory, University of Arizona,
Tucson, AZ 85721, USA  Lunar and Planetary Laboratory, University of Arizona,
Tucson, AZ 85721, USA

###### Abstract

The Rossby wave instability (RWI) is the fundamental non-axisymmetric radial
shear instability in disks. The RWI can facilitate disk accretion, set the
shape of planetary gaps and produce large vortices. It arises from density
and/or temperature features, such as radial gaps, bumps or steps. A general,
sufficient condition to trigger the RWI is lacking, which we address by
studying the linear RWI in a suite of simplified models, including
incompressible and compressible shearing sheets and global, cylindrical disks.
We focus on enthalpy amplitude and width as the fundamental properties of disk
features with various shapes. We find analytic results for the RWI boundary
and growth rates across a wide parameter space, in some cases with exact
derivations and in others as a description of numerical results. Features
wider than a scale-height generally become unstable about halfway to Rayleigh
instability, i.e. when the squared epicyclic frequency is about half the
Keplerian value, reinforcing our previous finding. RWI growth rates
approximately scale as enthalpy amplitude to the 1/3 power, with a weak
dependence on width, across much of parameter space. Global disk curvature
affects wide planetary gaps, making the outer gap edge more susceptible to the
RWI. Our simplified models are barotropic and height-integrated, but the main
results should carry over to more complex and realistic scenarios.

Astrophysical fluid dynamics(101) — Planet formation(1241) — Protoplanetary
disks(1300) — Hydrodynamics(1963)

††journal: ApJ

##  1 Introduction

The Rossby wave instability arises when radial disk structures, such as bumps
or gaps, induce strong pressure gradients and non-Keplerian radial shear
(Lovelace et al., 1999; Li et al., 2000). The RWI can generate large vortices
(Li et al., 2001), for instance at the edges of planetary gaps (de Val-Borro
et al., 2007), which affects planet migration (Lin & Papaloizou, 2010). The
RWI also helps transport matter falling onto accretion disks (Kuznetsova et
al., 2022).

Dust is trapped in both RWI-produced vortices and the rings that trigger the
RWI, in agreement with the disk structures observed by ALMA (Pinilla & Youdin,
2017). The RWI thus constrains observable rings and vortices (Chang et al.,
2023), for instance by regulating planet-carved gaps (Cimerman & Rafikov,
2023). Dust trapped in such rings and vortices can trigger planet formation
(Chiang & Youdin, 2010; Dra̧żkowska & Dullemond, 2018; Hu et al., 2018; Lyra
et al., 2024).

These significant consequences arise from simple considerations. The RWI does
not require vertical motions, baroclinicity or cooling, in contrast to the
vertical shear instability (Nelson et al., 2013; Lin & Youdin, 2015) and other
thermal disk instabilities (Lesur et al., 2023; Klahr et al., 2023). The RWI
can be triggered by zonal flows arising from these hydrodynamic (Manger et
al., 2020), or magnetohydrodynamic (MHD; Johansen et al., 2009),
instabilities. RWI analyses that include 3D motions (Meheut et al., 2012; Lin,
2013), cooling (Huang & Yu, 2022), dust feedback (Liu & Bai, 2023) and non-
ideal MHD (Cui et al., 2024) are crucial for a complete understanding, and
generally find modest corrections to idealized cases.

Even for simple cases, a general criterion for the onset of the RWI has been
elusive. Ono et al. (2016) found that the RWI was triggered partway between
the Lovelace and Rayleigh criteria, for a variety of barotropic disk features.
The Lovelace criterion, equivalent to a vortensity extrema in isentropic
disks, is necessary but insufficient for the RWI (Lovelace et al., 1999). The
Rayleigh criterion gives axisymmetric instability for disks with radially
decreasing angular momentum somewhere, i.e. negative squared epicyclic
frequency, $\kappa^{2}$.

Chang et al. (2023) found that disk bumps (barotropic and baroclinic)
triggered RWI when $\kappa^{2}$ was locally reduced to $\sim 60\%$ of the
Keplerian value. We colloquially refer to this criterion as “halfway to
Rayleigh” instability.

This work aims to develop a more fundamental understanding of the RWI boundary
and growth rates, including the “halfway to Rayleigh” criterion. We develop
scaling relations using the strength and width of disk features. We start with
simplified shearing sheet models and test against global disk models. This
approach is motivated by previous shearing sheet models studying
incompressible (Lithwick, 2007) and compressible (Vanon & Ogilvie, 2016) shear
instability, linear Rossby modes (Umurhan et al., 2016) and non-linear RWI
with cooling (Fung & Ono, 2021).

We present our method for studying the RWI with shearing sheet models in §2.
Sections 3 and 4 present our results for the incompressible and compressible
sheets, respectively. We compare to global disks in §5. A suggesting starting
point is the summary of our main results in §6.

##  2 Shearing Sheet RWI Models

###  2.1 The Compressible Shearing Sheet

The shearing sheet models a disk patch centered at radius $R_{\mathrm{c}}$,
rotating at the local Keplerian frequency, $\varOmega$, with cartesian $x,y,z$
coordinates oriented radially, azimuthally and vertically. Vertical averaging
gives the equations of motion (Goldreich & Tremaine, 1978; Johnson & Gammie,
2005)

$\displaystyle\frac{D\varSigma}{Dt}$
$\displaystyle=\varSigma\nabla\cdot\boldsymbol{v}$ (1a)
$\displaystyle\left(\frac{D}{Dt}+2\varOmega\hat{z}\times\right)\boldsymbol{v}$
$\displaystyle=3\varOmega^{2}x\hat{x}-\frac{1}{\varSigma}\nabla P$ (1b)
$\displaystyle D(P/\varSigma^{\gamma})/Dt$ $\displaystyle=0$ (1c)

for fluid velocity $\boldsymbol{v}$, surface density $\varSigma$, and (height-
averaged) pressure $P$, with $D/Dt=\partial/\partial
t+\boldsymbol{v}\cdot\nabla$. An ideal gas with adiabatic index $\gamma$,
adiabatic motions and no self-gravity or viscosity are assumed.

Combinging Eqs. (1a, 1b),

$\displaystyle\frac{Dq}{Dt}$
$\displaystyle=\frac{\nabla\varSigma\times\nabla{P}}{\varSigma^{3}}\cdot\hat{z}\,,$
(2)

shows that vortensity,
$q\equiv(2\varOmega+\hat{z}\cdot\nabla\times\boldsymbol{v})/\varSigma$, is
conserved in the absence of baroclinic effects.

We consider an axisymmetric equlibrium with linear perturbations (using $0,1$
subscripts, repectively) as
$\varSigma=\varSigma_{0}(x)+\varSigma_{1},\,P=P_{0}(x)+P_{1},\,\boldsymbol{v}=v_{0}(x)\hat{y}+u_{1}\hat{x}+v_{1}\hat{y}$.
Perturbed quantities have a Fourier dependence
$\propto\exp[\imath(k_{y}y-\omega{t})]$, and $x$-dependent amplitudes.

The equilibrium orbital motion is

$\displaystyle v_{0}$
$\displaystyle=-\frac{3}{2}\varOmega{x}+\Delta{v}_{0}=-\frac{3}{2}\varOmega{x}+\frac{1}{2\varOmega}\frac{d\varPi_{0}}{dx}$
(3)

where $\varPi_{0}$, the equilibrium enthalpy, $\varPi=\int{d}P/\varSigma$
gives the non-Keplerian motion, $\Delta{v}_{0}$. The equilibrium vortensity,
$q_{0}=\kappa^{2}/(2\varOmega\varSigma_{0})$, depends on the squared epicyclic
frequency:

$\displaystyle\kappa^{2}$
$\displaystyle=2\varOmega\left(2\varOmega+\frac{dv_{0}}{dx}\right)=\varOmega^{2}+\frac{d^{2}\varPi_{0}}{dx^{2}}\,.$
(4)

The linear equations of motion for the Fourier amplitudes (given the same
symbols as perturbed quantities for simplicity) are

$\displaystyle-\imath\Delta\omega\varSigma_{1}$
$\displaystyle=-\frac{d}{dx}\left(\varSigma_{0}u_{1}\right)-\imath
k_{y}\varSigma_{0}v_{1}$ (5a)
$\displaystyle-\imath\Delta\omega{u_{1}}-2\varOmega v_{1}$
$\displaystyle=-\frac{1}{\varSigma_{0}}\frac{dP_{1}}{dx}+\frac{\varSigma_{1}}{\varSigma_{0}^{2}}\frac{dP_{0}}{dx}$
(5b)
$\displaystyle-\imath\Delta\omega{v_{1}}+\frac{\kappa^{2}}{2\varOmega}u_{1}$
$\displaystyle=-\imath{k_{y}}\frac{P_{1}}{\varSigma_{0}}$ (5c)
$\displaystyle-\imath\Delta\omega\left(\frac{P_{1}}{P_{0}}-\gamma\frac{\varSigma_{1}}{\varSigma_{0}}\right)$
$\displaystyle=-u_{1}\frac{d\ln(P_{0}/\varSigma_{0}^{\gamma})}{dx}$ (5d)

with Doppler shifted frequency, $\Delta\omega\equiv\omega-v_{0}(x)k_{y}$. We
define a squared sound speed $c_{0}^{2}\equiv\gamma{P}_{0}/\varSigma_{0}$,
scale-height $H_{0}\equiv c_{0}/\varOmega$, (inverse) entropy lengthscale

$\displaystyle
L_{S}^{-1}\equiv\frac{1}{\gamma}\frac{d\ln(P_{0}/\varSigma_{0}^{\gamma})}{dx}$
(6)

and radial buoyancy frequency

$\displaystyle N^{2}$
$\displaystyle\equiv-\frac{1}{\gamma\varSigma_{0}}\frac{dP_{0}}{dx}\frac{d\ln(P_{0}/\varSigma_{0}^{\gamma})}{dx}=-\frac{c^{2}}{\gamma{L}_{S}}\frac{d\ln(P_{0})}{dx}\,.$
(7)

Manipulations yield an ODE for $\Psi\equiv{P}_{1}/\varSigma_{0}$,

$\displaystyle\Psi^{\prime\prime}+B(x)\Psi^{\prime}$ $\displaystyle=C(x)\Psi,$
(8)

the shearing sheet version of Eq. (15) in Li et al. (2000) with primes for
$x$-derivatives, $B\equiv{d}\ln\mathcal{F}/dx$ and

$\displaystyle\mathcal{F}$
$\displaystyle\equiv\frac{\varSigma_{0}\varOmega^{2}}{\kappa^{2}+N^{2}-\Delta\omega^{2}},$
(9a) $\displaystyle C$
$\displaystyle\equiv{k_{y}^{2}}+\frac{\varSigma_{0}}{\mathcal{F}H_{0}^{2}}+\frac{2\varOmega
k_{y}B}{\Delta\omega}+C_{2},$ (9b) $\displaystyle C_{2}$
$\displaystyle\equiv\frac{1-L_{S}^{\prime}}{L_{S}^{2}}+\frac{B}{L_{S}}+\frac{4\varOmega{k}_{y}}{\Delta\omega{L}_{S}}-\frac{k_{y}^{2}N^{2}}{\Delta\omega^{2}}.$
(9c)

This work considers isentropic equilibria with $C_{2}=1/L_{S}=N^{2}=0$.

The corotation resonance at $\Delta\omega=0$ defines a corotation radius,
$x_{c}$, where $\Re[\Delta\omega(x_{c})]=0$. At the Lindblad resonances, where
$\Delta\omega^{2}=\kappa^{2}+N^{2}$ and $1/\mathcal{F}=0$, $B$ is singular.

The Schrödinger form of Eq. (8) uses $\Xi=\sqrt{\mathcal{F}}\Psi$ to obtain
(Ono et al., 2016):

$\displaystyle\Xi^{\prime\prime}$ $\displaystyle=D(x)\Xi$ (10a) $\displaystyle
D$ $\displaystyle=\frac{B^{\prime}}{2}+\frac{B^{2}}{4}+C\,.$ (10b)

We solve Eq. (8), since $\Xi$ is singular at Linblad resonances, but $D$ is a
useful effective potential.

Figure 1: Radial profiles of shearing sheet equilibria for (_from left to
right:_) bumps, gaps, drops and jumps of width $W$ and amplitude
$\mathcal{J}=\Delta\varPi/(\varOmega W)^{2}=0.5,1,2$ (_blue, orange and green
curves_). _Top row:_ Enthalpy. _Middle row_ Epicyclic frequency squared,
$\kappa^{2}$. Minima (_dots_) and maxima (_pluses_) of vorticity (and
equivalently $\kappa^{2}$) are marked, in all rows. _Bottom row:_ The
effective potential $D_{\mathrm{inc}}$ for marginally stable RWI (offset by
$k_{y}^{2}$, see text).

###  2.2 The Incompressible Shearing Sheet

For the incompressible shearing sheet (Latter & Papaloizou, 2017) we take the
limit $\gamma\rightarrow\infty$, so that Eqs. (1a, 1c) give
$\nabla\cdot\boldsymbol{v}=0$. We replace $\nabla{P}/\varSigma=\nabla\varPi$
in equation (1b). The equilibrium is set by the choice of $\varPi_{0}(x)$,
from which $v_{0}(x)$ and $\kappa^{2}$ follow Eqs. (3, 4). The perturbed flow
obeys a stream function, $\psi$, as
$u_{1}=-\imath{k_{y}}\psi,\,v_{1}=\psi^{\prime}$.

The vorticity $\zeta=(\nabla\times\boldsymbol{v})\cdot\hat{z}$, with
equilibrium $\zeta_{0}=v_{0}^{\prime}(x)$ and perturbation
$\zeta_{1}=\psi^{\prime\prime}-k_{y}^{2}\psi$ is conserved $D\zeta/Dt=0$. Thus

$\displaystyle-\imath\Delta\omega\zeta_{1}$
$\displaystyle=-\zeta_{0}^{\prime}u_{1}$ (11)

which gives

$\displaystyle\psi^{\prime\prime}=\left(k_{y}^{2}+\frac{v_{0}^{\prime\prime}}{v_{0}-\omega/k_{y}}\right)\psi\equiv{D}_{\mathrm{inc}}(x)\psi\,,$
(12)

the famous Rayleigh equation for non-rotating incompressible shear flows.
Coriolis forces set $v_{0}$, but rotation is otherwise absent (see Lithwick,
2007).

There is a vast literature on this equation (Drazin & Reid, 2004). Relevant
results include Rayleigh’s theorem that a vorticity extrema,
$\zeta_{0}^{\prime}(x)=0$, is required for instability. Fjørtoft’s theorem
further states that this inflection point must be a maximum in
$|\zeta_{0}(x)|$. Since disks have $\zeta_{0}(x)<0$, instability requires a
(signed) vorticitiy minimum.

Fjørtoft’s theorem agrees with the interpretation of $D_{\mathrm{inc}}$ as a
potential, since for corotation at a vorticity minimum
$\Re(D_{\mathrm{inc}})<0$ near corotation for long wavelengths,
$k_{y}\rightarrow 0$. We further see that long wavelengths are the most
unstable. When applied to compressible, barotropic disks, a vortensity minimum
is required for instability.

Comparing to the compressible case, we might expect
$D\rightarrow{D}_{\mathrm{inc}}$ in some incompressible limit. Despite a
shared $k_{y}^{2}$ term, we find that for $k_{y}W\ll 1$, the incompressible
limit has $\mathcal{F}\propto 1/\kappa^{2}$ and
$2\varOmega{B}\rightarrow-4v_{0}^{\prime\prime}$. Thus the compressbile
corotation term is 4 times larger. This surprising result is possible since
$D$ has additional relevant terms and is a potential for a different fluid
quantity. Despite this difference, our compressible results have a well-
behaved incompressible limit (§4).

###  2.3 Disk Features

To understand the universal features of RWI, we consider various disk
structures, including bumps, gaps and step. Our compressible and
incompressible models share a common equilibrium enthalpy $\varPi_{0}(x)$ and
thus $v_{0}(x)$. Our parameterization

$\displaystyle\varPi_{0}(x)={\varPi}_{\mathrm{b}}+\Delta\varPi{S}(x/W)$ (13)

has two constants, the reference value ${\varPi}_{\mathrm{b}}$ (which only
affects compressible models) and amplitude $\Delta\varPi>0$. This work
considers the shapes:

$\displaystyle S(X)=\begin{cases}G(X)&\text{bump}\\\ 1-G(X)&\text{gap}\\\
\dfrac{1-\tanh(X)}{2}&\text{drop}\\\
\dfrac{1+\tanh(X)}{2}&\text{jump}\end{cases}$ (14)

for $G(X)=\exp(-X^{2}/2)$ and scaled distance $X\equiv{x}/W$. All shapes vary
from $0$ to $1$ over a radial width $\sim{W}$. Since
$\varPi_{\mathrm{b}}=\min(\varPi_{0}(x))>0$, all $\varPi_{0}(x)>0$.

We describe some properties of our shape functions next, then apply them to
compressible models in §2.3.2.

####  2.3.1 Shape Functions

Figure 1 plots (in the top row) our enthalpy features. The scaled amplitude

$\displaystyle\mathcal{J}\equiv\Delta\varPi/(\varOmega{W})^{2}$ (15)

measures a feature’s vorticity amplitude.

The middle row of Figure 1 plots $\kappa^{2}$. The location of vorticity (and
$\kappa^{2}$) minima is $x_{\mathrm{m}}=0$ for bumps, $\pm\sqrt{3}W$ for
gaps—which have a pair of vorticity minima—and $W\ln(2\pm\sqrt{3})/2\simeq\pm
1.32W$, for jumps and drops, respectively. The location of vortensity minima,
relevant for compressible flows, will be slightly shifted.

Inner and outer gap edges are symmetrically equivalent in the shearing sheet.
So are the drop and jump cases. Henceforth the “step” case refers to both.

Rayleigh instability occurs for $\min(\kappa^{2})<0$ and requires vertical
motions, absent from our model. The Rayleigh instability is still highly
relevant and occurs for $\mathcal{J}>\mathcal{J}_{\kappa}\gtrsim 1$.
Specifically, from Equations (4, 13)

$\displaystyle\frac{\kappa^{2}}{\varOmega^{2}}$
$\displaystyle=1+\mathcal{J}\frac{d^{2}S}{dX^{2}}\,.$ (16)

and $\mathcal{J}_{\kappa}\equiv 1/\max(-d^{2}S/dX^{2})=1$ for bumps, 2.241 for
gaps, and 2.598 for steps. Thus
$\min(\kappa^{2})=(1-\mathcal{J}/\mathcal{J}_{\kappa})\varOmega^{2}$ in the
shearing sheet. For global models, $\kappa^{2}$ also depends on
$W/R_{\mathrm{c}}$ (§5.2). This dependence vanishes in the shearing sheet
limit, $W/R_{\mathrm{c}}\ll 1$.

The bottom row of Figure 1 plots the incompressible effective potential as
$D_{\mathrm{inc}}-k_{y}^{2}$ (Eq. 12). The corotation radius is at a vorticity
minimum, $x_{\mathrm{m}}$, with phase speed
$c_{\omega}=\omega/k_{y}=v_{0}(x_{\mathrm{m}})$. This choice removes the
corotation singularity and gives (consistent with Fjørtoft’s theorem) a
negative potential well for trapped modes. The compressible potential $D$
behaves similarly, but waves also propagate exterior to Lindblad resonances,
where $D<0$ (see Fig. 2).

Figure 2: _Top left:_ Along the RWI boundary for $k_{y}H=0.1$ in a
compressible shearing sheet bump, the numbers
\raisebox{-0.9pt}{1}⃝-\raisebox{-0.9pt}{3}⃝ mark the modes investigated.
_Other panels:_ The effective potential $D$ (_orange curves, with negative
regions shaded_) and pressure perturbation $\Psi$ (_blue curves for magnitude
[thick], real and imaginary parts [thin solid and dashed]_) of the numbered
modes, with the bump width (_green dotted lines_) and nominal (Keplerian)
location of Lindblad resonances (_gray dotted lines_) marked.
\raisebox{-0.9pt}{1}⃝: Distant Lindblad resonances, with a trapped mode in the
Rossby zone. \raisebox{-0.9pt}{2}⃝: The trapped Rossby mode couples to density
waves exterior to nearby Lindblad resonances. \raisebox{-0.9pt}{3}⃝: No
Lindblad resonances since $\kappa^{2}<0$, and a “leaky” potential (negative
everywhere). This Rayleigh unstable region is not our focus.

####  2.3.2 Compressible Shearing Sheet Features

The compressible shearing sheet model requires not just
$\varPi_{0}^{\prime}(x)$ but also $\varSigma_{0}$ and $P_{0}$. We consider
polytropic models with
$P_{0}/P_{\mathrm{b}}=(\varSigma_{0}/\varSigma_{\mathrm{b}})^{\Gamma}$, with
reference values $\varSigma_{\mathrm{b}},\,P_{\mathrm{b}}$. The structure
index $\Gamma$ could differ from the adiabatic index $\gamma$ (but doesn’t
here, see below).

The polytropic enthalpy

$\displaystyle\varPi_{0}$
$\displaystyle=\int\frac{dP_{0}}{\varSigma_{0}}=\frac{\Gamma}{\Gamma-1}\frac{P_{0}}{\varSigma_{0}}$
(17)

matches Eq. (13) for

$\displaystyle\varSigma_{0}$
$\displaystyle={\varSigma}_{\mathrm{b}}\left[1+\frac{\Delta\varPi}{{\varPi}_{\mathrm{b}}}S(x/W)\right]^{\frac{1}{\Gamma-1}}\,,$
(18)

and
${\varPi}_{\mathrm{b}}=\Gamma{P}_{\mathrm{b}}/[(\Gamma-1)\varSigma_{\mathrm{b}}]$.

This compressible polytropic model requires 3 additional parameters, besides
$k_{y}W$ and $\mathcal{J}$: $\gamma,\,\Gamma$ and $H\equiv{c}/\varOmega$ with

$\displaystyle{c}^{2}\equiv\gamma{P}_{\mathrm{b}}/{\varSigma}_{\mathrm{b}}=\gamma(\Gamma-1)\varPi_{\mathrm{b}}/\Gamma\,.$
(19)

We drop $\mathrm{b}$ subscripts from reference $H$ and $c$ values for
convenience. We don’t need $\varSigma_{\mathrm{b}}$ or $P_{\mathrm{b}}$
independently, as Eq. (8) only depends on logarithmic derivatives of
$\varSigma_{0}$ and $P_{0}$.

To reduce parameter space, we fix $\Gamma=\gamma=4/3$ for an adiabatic sheet
with $N^{2}=1/L_{S}=0$. A diatomic gas with $\gamma_{\mathrm{3D}}=7/5$
corresponds to our height integrated
$\gamma=(3\gamma_{\mathrm{3D}}-1)/(\gamma_{\mathrm{3D}}+1)=4/3$ (Goldreich et
al., 1986; Li et al., 2000). Thus $H$ is the only additional free parameter
our compressible models.

The limits $\Gamma\rightarrow 1,\infty$ describe constant temperature and
$\varSigma_{0}$ features, respectively (Chang et al., 2023). For completeness,
the $\Gamma\rightarrow 1$ limit of Eq. (18) is

$\displaystyle\varSigma_{0}$
$\displaystyle=\gamma\frac{P_{0}}{c^{2}}={\varSigma}_{\mathrm{b}}\exp\left[\frac{\gamma\Delta\varPi}{c^{2}}S\left(\frac{x}{W}\right)\right]\,.$
(20)

with $\varPi_{\mathrm{b}}(\Gamma-1)\rightarrow{c}^{2}/\gamma$ remaining
finite.

###  2.4 Boundary Conditions and Solution Methods

Solving our second order ODEs requires a pair of boundary conditions, applied
at large distances $|x|\gg{W},\,1/k_{y}$, and (for the compressible case)
$|x|\gg{H}$.

For the incompressible case (Eq. 12), $D_{\mathrm{inc}}\rightarrow{k}_{y}^{2}$
at large $|x|$. Physical solutions decay exponentially, with boundary
conditions,

$\displaystyle\psi^{\prime}=\pm{k}_{y}\psi\,,$ (21)

at large $\mp|x|$.

For the compressible case, boundary conditions exterior to the Lindblad
resonances should match onto outgoing density waves. We seek WKB solutions of
the form $\Psi\sim{A}(x)\exp(\imath\int^{x}k_{x}(\chi)d\chi)$. Compared to
previous works (Li et al., 2000; Ono et al., 2016; Chang et al., 2023) who
used just $k_{x}$, we find $A(x)$ to lowest order, which improves some
numerical results.

First we confirm that outgoing waves have $k_{x}(x)>0$. The large $|x|$,
Keplerian limit gives $\Delta\omega\rightarrow
3\varOmega{k}_{y}x/2,\,\mathcal{F}/\varSigma_{0}\rightarrow-\varOmega^{2}/\Delta\omega^{2},\,B\rightarrow-3\varOmega{k}_{y}/\Delta\omega\rightarrow-2/x$
and
$C\simeq{D}\rightarrow-(\Delta\omega/c_{0})^{2}\rightarrow[3k_{y}x/(2H_{0})]^{2}$.
To lowest WKB order, Eq. (8) gives $k_{x}=\pm\sqrt{-C}$, i.e.
$\Delta\omega^{2}=(k_{x}c_{0})^{2}$. The group velocity

$\displaystyle\frac{\partial\omega}{\partial{k}_{x}}$
$\displaystyle=\frac{\partial\Delta\omega}{\partial{k}_{x}}\approx\frac{k_{x}c^{2}}{\Delta\omega}\approx\frac{2k_{x}H_{0}}{3k_{y}x}c$
(22)

confirms that $k_{x}>0$ for outgoing waves (as $k_{y}>0$ by convention).

For more accuracy, we adopt the physical optics solution to Eq. (10a),

$\displaystyle\Xi$
$\displaystyle\sim\frac{c_{\Xi}}{\sqrt{k_{x,D}}}\exp\left(\imath{\int}^{x}k_{x,D}(\chi)d\chi\right)$
(23)

with $k_{x,D}=\sqrt{-D}$ (the desired positive root) and $c_{\Xi}$ an
arbitrary (complex) constant.

Taking the derivative gives the boundary condition

$\displaystyle\Xi^{\prime}$
$\displaystyle=\left(\imath\sqrt{-D}-\frac{1}{4D}\frac{dD}{dx}\right)\Xi\,.$
(24)

The desired boundary condition for $\Psi=\Xi/\sqrt{\mathcal{F}}$ follows as

$\displaystyle\Psi^{\prime}$
$\displaystyle=\left(\imath\sqrt{-D}-\frac{B}{2}-\frac{1}{4D}\frac{dD}{dx}\right)\Psi\,.$
(25)

At large $|x|$, $d\ln(D)/dx/4\rightarrow 1/(2x)$, so that $|\Xi|\propto
1/|x|^{1/2}$ and $|\Psi|\propto|x|^{1/2}$, in agreement with our numerical
solutions.

Our numerical solutions use the shooting method. At the inner boundary,
$x_{i}$, we pick an arbitrary $\psi(x_{i})$ or $\Psi(x_{i})$ and set the
derivative with the boundary condition, Eq. (21) or Eq. (25). We integrate
with the Dormand-Prince method (“DOP853”) implemented in
scipy.integrate.solve_ivp. The integrated solution deviates from the outer
boundary condition. Using Muller’s method, we minimize the residual error and
find the complex eigenvalue $\omega\equiv\omega_{r}+\imath{s}$. The shooting
method requires good initial guesses. We use known solutions to gradually
explore parameter space.

For global models, we apply the same method but solve Eq. (15) in Li et al.
(2000) instead of Eq. (8).

We have validated our numerical result several ways, including adjusting the
outer boundary positions, finding the incompressible limit of compressible
results and using different methods for the RWI stability boundary (below).
Similar to Li et al. (2000), we derive an energy equation from Eq. (5), which
after azimuthal averaging (denoted by brackets) is

$\displaystyle\frac{\partial}{\partial{t}}\left[\frac{\varSigma_{0}}{2}\left(\langle|\boldsymbol{v}_{1}|^{2}\rangle+\frac{\langle\Psi^{2}\rangle}{c^{2}}\right)\right]$
$\displaystyle=$ (26)
$\displaystyle-\frac{dv_{0}}{dx}\varSigma_{0}\langle{u}_{1}v_{1}\rangle-\frac{d}{dx}\langle{P}_{1}u_{1}\rangle$

We verified that growth rates and eigenfunctions found by our numerical method
satisfy this relation, over a range of parameters $\mathcal{J},\,W/H,\,$ and
$k_{y}H$.

####  2.4.1 Locating the Stability Boundary

We find marginally stable modes using a simplified method (Ono et al., 2016).
With $s=0$, we fix the corotation radius, $x_{\mathrm{c}}$, to vorticity (or
vortensity) minima for incompressible (or compressible) models, which sets
$\omega=k_{y}v_{0}(x_{\mathrm{c}})$. With this choice, the ODE has real
coefficients, and no corotation singularity, as shown in Figure 1 for the
incompressible case.

One physical parameter, usually $k_{y}W$ or $W/H$, varies as the shooting
parameter (and eigenvalue). With other parameters held fixed, this method
finds marginally stable solutions. For the incompressible model, this method
uses $\psi(x)$ purely real. For the compressible models, $\Psi(x)$ has a
complex boundary condition (Eqs. 24, 25). Thus the eigenvalue ($W/H$) can
acquire an imaginary part, which is unphysical. Usually this imaginary part is
negligibly small ($\lesssim 10^{-3}$ of the real part), which validates the
method. Fig. 2 shows example solutions obtained with this method. Growth rates
away from the stability boundary are also mapped, using the usual method.

For more extreme parameters—near Rayleigh instability and for $k_{y}W\simeq 1$
(placing Lindblad resonances in the Rossby zone)—this method can fail, as it
does (for different reasons) in baroclinic disks (Chang et al., 2023). In
these cases, we simply measure where $s$ drops to small values.

It is numerically difficult to find growth rates with $s/\Omega\lesssim
10^{-3}$. Since both methods agree on the stability boundary location (when
this simplified method works), the stability boundary is relatively sharp.

Figure 3: Incompressible RWI growth rate $s$ for (left to right) bump, gap and
step features against scaled enthalpy amplitude
$\mathcal{J}\equiv\Delta\varPi/(W^{2}\varOmega^{2})$ and the ratio of feature
width to azimuthal wavelength, $k_{y}W$. The RWI boundary (dotted yellow),
$\mathcal{J}=k_{y}W$ reference line (dashed yellow), and Rayleigh unstable
regions (gray shaded) are shown.

##  3 Incompressible Results

Figure 3 maps RWI growth rates for various shapes in the incompressible
shearing sheet. For a given shape, the incompressible RWI is completely
described by the parameters for amplitude,
$\mathcal{J}=\Delta\varPi/(\varOmega W)^{2}$, and width (scaled to
wavenumber), $k_{y}W$. We describe the incompressible stability boundary,
growth rates, and eigenfunctions below.

###  3.1 Incompressible stability boundary

The dotted yellow curves in Figure 3 show the stability boundary, found as
described in section 2.4.1. RWI occurs for larger $\mathcal{J}$ or smaller
$k_{y}W$ than this boundary, and no modes (stable, unstable or damped) exist
on the other side.

The stability boundary is best understood as smoothly connected
$\mathcal{J}\ll 1$ and ${J}\gg 1$ limits. For $\mathcal{J}\ll 1$, the
stability boundary follows $\mathcal{J}\simeq{f}_{\mathrm{MS}}k_{y}W$, or

$\displaystyle\Delta\varPi$
$\displaystyle=f_{\mathrm{MS}}\varOmega^{2}k_{y}W^{3}$ (27)

with $f_{\mathrm{MS}}\simeq 1.20,\,2.39,\,2.65$ for the bumps, gaps and steps,
respectively.

For $\mathcal{J}\gg 1$, the stability boundary is simply
$k_{y}W=g_{\mathrm{MS}}$, with $g_{\mathrm{MS}}=\sqrt{2},1.05,2.0$ for the
bump, gap and step cases, respectively. While large $\mathcal{J}$ values are
Rayleigh unstable, this limiting behavior explains why the stability curve
steepens for $\mathcal{J}\gtrsim 1$.

These limiting behaviors can be understood several ways, described below.

Figure 4: Similar to Figure 3, except incompressible growth rates are scaled
by a characteristic rate $(\varOmega k_{y}^{2}\Delta\varPi)^{1/3}$. Away from
the stability boundary, these scaled rates vary only moderately.

####  3.1.1 Intuitive Explanations

The $\mathcal{J}\gg 1$ instability condition, $k_{y}W<g_{\mathrm{MS}}$ follows
the idea that counter-propagating Rossby waves (CRWs) drive shear instability
(Heifetz et al., 1999). For simplicity, we consider bumps and examine the
approximate condition for CRWs at $x\simeq\pm{W}$ to maintain stationary
phase, with phase speed $c_{\omega}=\omega/k_{y}=0$, as illustrated in §3.3.

With a $\psi(x)\propto\exp(\imath{k}_{x}x)$ WKB approximation, Eq. (12) gives
$c_{\omega}=v_{0}+v_{0}^{\prime\prime}/(k_{x}^{2}+k_{y}^{2})$. At $x=\pm{W}$,
$v_{0}/(\varOmega W)\sim\mp(1+\mathcal{J})$ roughly accounts for Keplerian and
non-Keplerian flow, and $v_{0}^{\prime\prime}\sim\pm\mathcal{J}\varOmega/W$.
Taking $k_{x}W\simeq 1$ matches the local wave packet to feature size, giving

$\displaystyle\frac{v_{0}^{\prime\prime}}{k_{x}^{2}+k_{y}^{2}}$
$\displaystyle\sim\pm\frac{\mathcal{J}}{1+(k_{y}W)^{2}}(\varOmega{W})\,.$ (28)

Thus $c_{\omega}=0$ requires $\mathcal{J}\sim(1+\mathcal{J})(1+(k_{y}W)^{2})$.

For $\mathcal{J}\gg 1$ this rough analysis requires $1\sim 1+(k_{y}W)^{2}$ or
$k_{y}W\lesssim 1$ for phase matching and instability, as desired. For
$\mathcal{J}\ll 1$, this analysis fails.

Instead, for the $\mathcal{J}\ll 1$ boundary, another WKB analysis applies.
Since $k_{y}W\ll 1$, waves have a shallow decay at large $|x|/W$, as
$\psi\propto\exp(-k_{y}|x|)$. To match onto this decay, the slope across the
Rossby zone must change sign, but only change magnitude by a small amount,
$\Delta\Phi\equiv{W}\psi^{\prime}|_{-W}^{W}/\psi\sim-k_{y}W$.

Across corotation, the slope change from WKB oscillations,
$\psi\propto\exp(\imath\sqrt{-D_{\mathrm{inc}}}x)$, is

$\displaystyle\Delta\Phi$
$\displaystyle=W\int_{-W}^{W}\psi^{\prime\prime}dx/\psi\sim{D}_{\mathrm{inc}}W^{2}\sim-\mathcal{J},$
(29)

Where the depth of the potential near corotation
$D_{\mathrm{inc}}\simeq-\mathcal{J}/W^{2}$ (Fig. 1). A trapped mode thus
requires $\mathcal{J}\sim{k}_{y}W$, in agreement with the stability boundary.
The small change in wave phase
$\sqrt{-D_{\mathrm{inc}}}W\sim\sqrt{\mathcal{J}}\ll 1$ explains the failure of
standard WKB theory for $\mathcal{J}\ll 1$, as noted above.

For a more physical explanation of the $\mathcal{J}\ll 1$ stability boundary,
we briefly summarize the analysis of shearing waves by Lithwick (2007).
Shearing waves interact with axisymmetric disk features of width $W$ and
vorticity amplitude $\Delta\zeta_{0}\simeq\Delta\Pi/(\varOmega{W}^{2})$.

A leading wave with initial radial wavenumber $k_{x}(t=0)\simeq-1/W$ and fixed
$k_{y}>0$ swings through a radial orientation, $k_{x}(t_{\mathrm{sw}})=0$, in
time $t_{\mathrm{sw}}=(2k_{x}(0)/(3\varOmega{k}_{y})\simeq
1/(\varOmega{k}_{y}W)$, since for $\mathcal{J}\ll 1$ the flow is nearly
Keplerian (Goldreich & Lynden-Bell, 1965). While swinging, the wave couples to
the disk feature and spawns a new leading wave. The amplitude of successive
waves increases if $\Delta\zeta_{0}t_{\mathrm{sw}}\gtrsim 1$ or
$\Delta\varPi\gtrsim\varOmega^{2}k_{y}W^{3}$, reproducing the $\mathcal{J}\ll
1$ instability criterion.

####  3.1.2 More quantitative explanations

The above arguments can be made more rigorous. For $\mathcal{J}\ll 1$,
Lithwick (2007) couples the physical argument ($\Delta\zeta
t_{\mathrm{sw}}\gtrsim 1$) to a stability boundary given by the integral

$\displaystyle k_{y}\varOmega$
$\displaystyle=\frac{1}{3}\int_{-\infty}^{\infty}\frac{d\zeta_{0}/dx}{x-x_{\mathrm{c}}}\,dx$
(30)

with the vorticity minimum at $x_{\mathrm{c}}$. This result reproduces Eq.
(27), and, for our shapes, precisely gives $f_{\mathrm{MS}}=6/I_{\mathrm{MS}}$
with

$\displaystyle
I_{\mathrm{MS}}=\int_{-\infty}^{\infty}\frac{S^{\prime\prime\prime}(X)}{X-X_{\mathrm{c}}}\,dX$
(31)

where $X_{\mathrm{c}}=x_{\mathrm{c}}/W$. Integrating $I_{\mathrm{MS}}$
reproduces our numerical results. For bumps, $f_{\mathrm{MS}}=3/\sqrt{2\pi}$,
and gaps, $f_{\mathrm{MS}}=6/\sqrt{2\pi}$. The numerically integrated
$I_{\mathrm{MS}}$ for steps is also consistent.

For $\mathcal{J}\gg 1$, the stability boundary $k_{y}W=\sqrt{2}$ for the bump
case can be derived exactly. The $\mathcal{J}\rightarrow\infty$ limit gives a
parabolic potential well
$D_{\mathrm{inc}}W^{2}\rightarrow(k_{y}W)^{2}+3-(x/W)^{2}$. This potential has
quantized bound states of “energy” $E=3-(k_{y}W)^{2}=2n+1$ for $n=0,1,$… (e.g.
Manasreh, 2012, Appendix E). For $(k_{y}W)^{2}>0$ only the $n=0$ bound state
exists, which demonstates the lack of RWI modes with higher radial order. This
bound state has $k_{y}W=\sqrt{2}$, as claimed.

For all shapes, a necessary condition for RWI follows from the requirement
that $D_{\mathrm{inc}}<0$. This necessary condition is only close to the
stability boundary for $\mathcal{J}\gtrsim 1$. For $\mathcal{J}\gg 1$ this
necessary condition is $k_{y}W<\sqrt{3}$ for gaps and bumps and $k_{y}W<2$ for
steps. These simple necessary conditions are close to, but less strict than,
the sufficient conditions for instability.

Figure 5: A map of perturbed vorticity, $\zeta_{1}$, for the RWI in the
incompressible shearing sheet, versus $x,y$ position, with arrows for the the
perturbed velocity $\boldsymbol{v}_{1}$. Rows from _top to bottom_ have growth
rates $s/\varOmega\approx 10^{-3},10^{-2},$ and $10^{-1}$, respectively.
Columns from _left to right_ consider enthalpy bumps, gaps and steps,
respectively. Solid vertical lines mark the feature width, $x=\pm{W}$. Dotted
vertical lines locate the minima of the equilibrium vorticity, $\zeta_{0}$.
The enthalpy amplitude $\mathcal{J}$ was chosen to produce the desired growth,
with wavenumber $k_{y}=0.1/W$ in all cases.

###  3.2 Incompressible growth rates

Figure 3 shows that RWI growth rates increase with $\mathcal{J}$ and with
$k_{y}W$, except near marginal stability. While it is easier to trigger RWI
for long wavelengths (small $k_{y}W$), growth rates well into the unstable
region are faster for smaller wavelengths (large $k_{y}W$). The smooth
variations in growth rates, especially away from the stability boundary,
suggest an analytic scaling.

Figure 4 plots growth rates scaled by the characteristic rate

$\displaystyle s_{\mathrm{inc}}$
$\displaystyle\equiv(\varOmega{k}_{y}^{2}\Delta\varPi)^{1/3}\,.$ (32)

This approximation is reasonably good, aside from the rapid decay near the
stability boundary.

In the $k_{y}W\ll 1$ limit, Lithwick (2007) derives the growth rate
$s/\varOmega=(3/2)\alpha{k}_{y}W$, where $\alpha$ follows from the complex
integral constraint on $\alpha$ and $\beta$:

$\displaystyle k_{y}\varOmega$
$\displaystyle=\frac{1}{3W}\int_{-\infty}^{\infty}\frac{d\zeta_{0}/dX}{(X-\beta)^{2}+\alpha^{2}}(X-\beta+\imath\alpha)\,dX.$
(33)

This result reduces to Eq. (30) for marginal stability, where
$\beta\rightarrow{X}_{\mathrm{c}}$. Analysis is simplest for bumps, where
symmetry about the vorticity minimum gives $\beta=0$, and the imaginary part
of the integral vanishes.

Our parameterization, with $S_{\mathrm{B}}$ for the bump shape, gives:

$\displaystyle\frac{k_{y}W}{\mathcal{J}}$
$\displaystyle=\frac{1}{f_{\mathrm{US}}(\alpha)}\equiv\frac{1}{6}\int_{-\infty}^{\infty}\frac{XS_{\mathrm{B}}^{\prime\prime\prime}(X)}{X^{2}+\alpha^{2}}\,dX,$
(34)

The function $f_{\mathrm{US}}(\alpha)$ for unstable modes, gives
$f_{\mathrm{US}}(0)=f_{\mathrm{MS}}$ at marginal stability. Eq. (34) only
gives simple expressions in limiting cases.

For $\alpha\ll 1$, $\pi\alpha/2\rightarrow
1/f_{\mathrm{MS}}-k_{y}W/\mathcal{J}$ gives the rise in growth rates near the
stability boundary as

$\displaystyle\frac{s}{\varOmega}$
$\displaystyle\approx\frac{3}{\pi}k_{y}W\left(\frac{1}{f_{\mathrm{MS}}}-\frac{k_{y}W}{\mathcal{J}}\right)$
(35)

The $\alpha\rightarrow\infty$ limit gives
$f_{\mathrm{US}}\rightarrow\alpha^{4}/\sqrt{2\pi}$ and

$\displaystyle\frac{s}{\varOmega}\rightarrow\frac{3}{2}(2\pi)^{1/8}\left(\frac{\Delta\varPi{k}_{y}^{3}W}{\varOmega^{2}}\right)^{1/4}$
(36)

Unfortunately this limit does not directly apply. We are mainly interested in
$k_{y}W/\mathcal{J}\gtrsim 0.01$, corresponding to $\alpha\lesssim 3.4$, i.e.
at most order unity.

Our approximate Eq. (32) corresponds to $f_{\mathrm{US}}\sim\alpha^{3}$, a
good approximation for order unity $\alpha$.

###  3.3 Incompressible Eigenfunctions

To visualize RWI modes, Figure 5 maps perturbed vorticity $\zeta_{1}$ and
velocity vectors $\boldsymbol{v}_{1}$ for various growth rates and feature
types. These incompressible eigenfunctions are similar to the standard global,
compressible RWI (Ono et al., 2016). Corotation is near the vorticity minima
marked with the dotted line.

The RWI mechanism is clearest for larger growth rates (bottom row). The pair
of CRWs across corotation (analyzed in §3.1.1) is evident. This wave pair is
shifted in azimuthal phase, but radially symmetric for the bumps, and
asymmetric for gaps and steps, consistent with their asymmetric potentials
(Fig. 1). The aximuthal phase shift causes flow though the vorticity minima to
primarily enter regions of negative perturbed vorticity. This explanation of
the growth mechanism is well known for general shear flows (Heifetz et al.,
1999) and the RWI (Ono et al., 2016).

At lower growth rates (middle row), the phase shift decreases. The feeding of
negative vorticity from the background into perturbations is less direct,
entering narrow fingers near corotation. Even closer to marginal stability
(top rows) the phase shift is nearly gone, and feeding via narrow glitches
near corotation is harder to see.

For marginal stability, $\zeta_{1}$ is non-zero and smooth through
corototation. However all growing modes have $\zeta_{1}=0$ at vorticity
extrema (as shown by Eq. 11). This fact explains the necessity of small
glitches near marginal stability, and the width of prominent CRWs, which fit
between a vorticity maximum and minimum (cf. Fig. 1).

##  4 Compressible shearing sheet results

Figure 6: Marginal stability curves for the RWI of bumps in the compressible
shearing sheet, for different values of $k_{y}H$. The incompressible,
$k_{y}H\rightarrow\infty$, limit (_dotted yellow curve_) and
$\min(\kappa^{2})=0.6$ reference (_pink line with x’s_) are shown. Axes and
Rayleigh unstable region as in Fig. 3.

We now analyze the compressible shearing sheet model of sections 2.1, 2.3.2.
Compressible effects are captured by the value of $k_{y}H$ (see below). Our
incompressible results roughly correspond to the $k_{y}H\rightarrow\infty$
limit.

For an effective Mach number, we use the Keplerian shear across a lengthscale
of $1/k_{y}$ to define

$\displaystyle\mathcal{M}_{\mathrm{eff}}\equiv\frac{\varOmega/k_{y}}{c}=\frac{1}{k_{y}H}\,.$
(37)

RWI modes with $\mathcal{M}_{\mathrm{eff}}\lesssim 1$ behave incompressibly,
which is expected of subsonic flows. In global protoplanetary disks, the RWI
is moderately compressible for $m=1$ modes, and more incompressible for higher
$m$ (§6, point 3).

###  4.1 Compressible Stability Boundary

Figure 6 shows the effect of compressibility, measured by $k_{y}H$, on the RWI
boundary. The bump feature is chosen, and is representative, with quantitative
shape effects are noted below.

For $k_{y}H=1$ the stability boundary overlaps the incompressible limit
($k_{y}H\rightarrow\infty$). As $k_{y}H$ decreases, compressibility effects
increase.

For sufficiently small $k_{y}H\lesssim 0.01$, the stability boundary breaks
into three distinct regions, approximately:

$\displaystyle\frac{\Delta\varPi}{\varOmega^{2}W^{2}}$
$\displaystyle\approx\begin{cases}f_{\mathrm{MS}}W\left(k_{y}+\frac{1}{4H}\right)&\text{if}\,\,W\lesssim{H}\\\
0.4j_{\mathrm{MS}}&\text{if}\,\,H\lesssim{W}\lesssim{k}_{y}^{-1}\\\
\infty&\text{if}\,\,k_{y}W\gtrsim{g}_{\mathrm{MS}}.\\\ \end{cases}$ (38)

The shape-dependent factors $f_{\mathrm{MS}},\,g_{\mathrm{MS}}$ (§3.1) and
$j_{\mathrm{MS}}$ (below) are order unity. For incompressible parameters,
$k_{y}H\gtrsim 0.3$, this stability boundary reverts to the incompressible
case, with no intermediate width region. For marginal compressibility,
$k_{y}H\simeq 0.1$, these regions are not as distinct, with overlapping
transitions.

For small widths, $W<H$, the compressible ($k_{y}H\ll 1$) stability boundary
follows $\Delta\varPi\approx{f}_{\mathrm{MS}}\varOmega^{2}W^{3}/(4H)$,
independent of $k_{y}$. Compared to the incompressible $k_{y}W\ll 1$ boundary,
$\Delta\varPi\propto{W}^{3}$ is identical, but compressible enthalpy features
must be $\simeq 1/(4k_{y}H)$ larger for instability. This stabilizing effect
generally arises from the fact that some of the energy is used to compress the
flow (Blumen, 1970).

Wide features and/or short wavelength modes,
$k_{y}W\gtrsim{g}_{\mathrm{MS}}\sim 1$, are RWI-stable, like the
incompressible case. However, compressible modes are more unstable between
$0.3\lesssim{k_{y}}W\lesssim{g}_{\mathrm{MS}}$ (see Fig. 6). This effect
arises because Lindblad resonances, absent from the incompressible limit,
approach the Rossby zone, as described below. Ultimately, the widest features
require a global treatment (§5).

For intermediate widths, with $W\gtrsim{H}$ but $W\lesssim{k}_{y}^{-1}$, the
stability boundary is approximately given by
$\mathcal{J}=\Delta\varPi/(W^{2}\varOmega^{2})=0.4j_{\mathrm{MS}}$ with
$j_{\mathrm{MS}}\simeq 1,2.8,3.0$ for bumps, gaps and steps, respectively. The
value of $\min(\kappa^{2})=1-\mathcal{J}/\mathcal{J}_{\kappa}\simeq
0.6,0.51,0.54$ for bumps, gaps and steps (respectively) is more similar,
emphasizing that $\kappa^{2}$, and being “halfway to Rayleigh” is more
fundamental.

Figure 7: For marginal RWI, the minimum value of $\kappa^{2}/\varOmega^{2}$
caused by bumps, gaps and steps (_blue, orange and green curves_ ,
respectively) vs. feature width. _Solid_ curves show compressible shearing
sheet models for $k_{y}H=0.03$. Global models of bumps, outer gap edges and
jumps (_dashed curves_) and of inner gap edges and drops (_dot-dashed curves_)
have $m=1,\,h=0.03$ (matching $mh=k_{y}H$). Global models break the symmetry
between inner and outer gap edges and between jumps and drops.

Figure 7 plots $\min(\kappa^{2})/\varOmega^{2}$ for marginal stability,
moderate compressibility, $k_{y}H=0.03$, and different shapes. (The global
models in this figure are discussed in §5.3.3.) For $W/H\gtrsim 2$, the
stability boundary is “halfway to Rayleigh” with $\min(\kappa^{2})\simeq
0.5{-}0.6\varOmega^{2}$. For stronger stronger compressibility (smaller
$k_{y}H$), $\min(\kappa^{2})$ values would be more strictly constant (Figure
6).

In Figure 7, the $k_{y}W\sim 1$ stability boundary is off scale at $W/H\sim
30$. For $W/H\lesssim 1$, the stability boundary approaches
$\min(\kappa^{2})/\varOmega^{2}=1-(f_{\mathrm{MS}}/\mathcal{J}_{\kappa})(W/H)/4$,
following Eq. (38).

While shape effects are minor in the shearing sheet, bumps most readily
trigger RWI, at larger $\min(\kappa^{2})$ values (Fig. 7) and smaller enthalpy
amplitudes (smaller $f_{\mathrm{MS}}$ and $j_{\mathrm{MS}}$ values, see Fig.
3).

We next examine origin of the three limits in Eq. (38).

####  4.1.1 Small width compressible boundary

In §3.1.1, incompressible instability for $k_{y}W\ll 1$ is given as the
Lithwick (2007) wave shearing time criteria
$\Delta\zeta_{0}t_{\mathrm{sh}}\gtrsim 1$. The corresponding $W\ll{H}$
compressible instability criterion is that the sound crossing time
$t_{\mathrm{sc}}\equiv{W}/(\varOmega H)\lesssim\Delta\zeta_{0}/\varOmega^{2}$.
Rotation appears explicitly in the compressible (but not the incompressible)
instability condition, consistent with the discussion after Eq. (12).

We can also adapt the WKB derivation of the $\mathcal{J},\,k_{y}W\ll 1$
incompressible stability boundary in §3.1.1 to the compressible $k_{y}H\ll
1,\,W\lesssim{H}$ case. The main difference in this case is that, from Eq.
(9b), the decay outside the Rossby zone follows $\Psi\propto\exp(-|x|/H)$.
Thus the slope change across corotation (now for $\Psi$) becomes
$\Delta\Phi\sim-W/H$. In this limit, the potential depths are similar (to
order unity), so the induced $\Delta\Phi\sim-\mathcal{J}$. Matching these two
gives $\mathcal{J}\sim{W}/H$, the desired compressible boundary.

####  4.1.2 Intermediate width stability boundary

The “halfway to Rayleigh” instability criterion is given above as the scaled
enthalpy condition, $\mathcal{J}\gtrsim 0.4j_{\mathrm{MS}}\sim 1$. In absolute
terms, relative to $\varPi_{\rm b}\sim c^{2}$, this condition becomes
$\Delta\varPi/c^{2}\gtrsim(W/H)^{2}$. Thus widths larger than $H$ require
increasingly strong enthalpy features, a relevant point for the astrophysical
origin of these features.

To explain this stability boundary, a negative potential at corotation
$D(x_{\mathrm{c}})<0$ gives a useful necessary condition for instability
(similar to §3.1.2). While simple to state, there are many terms in $D$ to
evaluate. These terms are stabilizing (or destabilizing) if they make a
positive (or negative) contribution to $D(x_{\rm c})$.

We focus on the bump case with $x_{\mathrm{c}}=0$ for simplicity, and take the
$k_{y}W\ll 1$ and $\Delta\varPi\gg{c}^{2}$ (equivalent to $W\gg{H}$ as noted
above) limits. These limits avoid the transitions to neighboring stability
regimes.

The main stabilizing term is $\varSigma_{0}/(\mathcal{F}H_{0}^{2})$, the usual
source of the corotation barrier in disks. With
$\varSigma_{0}/\mathcal{F}\propto\kappa^{2}$ this term is reduced near
Rayleigh instability, and also, via $H_{0}$, by disk heating. In our limits

$\displaystyle\left.\frac{\varSigma_{0}}{\mathcal{F}H_{0}^{2}}\right|_{x=0}$
$\displaystyle\rightarrow\frac{3}{W^{2}}\left(\frac{1}{\mathcal{J}}-1\right)$
(39)

The main destabilizing term is the corotation term
$C_{\mathrm{cor}}=2\varOmega{k}_{y}B/\Delta\omega$, though $B^{\prime}/2$ also
contributes. At corotation $B\propto{d}\ln(q_{0})/dx$ diverges approaching
Rayleigh instability.

Thus a simple explanation of the “halfway to Rayleigh” result is that the
stabilizing corotation barrier vanishes, and the destabilizing corotation
resonance diverges, for $\kappa^{2}\rightarrow 0$. Thus instability occurs
somewhat before this point.

Our limits give

$\displaystyle C_{\mathrm{cor}}(0)+\frac{B^{\prime}(0)}{2}$
$\displaystyle\rightarrow-\frac{3}{W^{2}}\left(\frac{3}{2(1-\mathcal{J})}+\frac{1}{3+\mathcal{J}}\right)$
(40)

with the advertised $\mathcal{J}\rightarrow 1$ divergence.

Combining Eqs. (39, 40), the necessary criterion $D(0)<0$ becomes
$\mathcal{J}>0.29$. This condition is close to, but naturally below, the
sufficient condition for bumps, $\mathcal{J}>0.4$.

One insight from this analysis is that equation of state effects should have a
modest effect on this stability boundary, via $H_{0}$ and
$\varSigma_{0}^{\prime}$. However $\kappa^{2}$ is the dominant effect. We
defer a more detailed study of thermodynamic, including baroclinic, effects.

####  4.1.3 Large width stability boundary

The $k_{y}W\gtrsim 1$ condition for stability matches the incompressible case
which was physically justified in §3.1.1. We do generalize that argument to
include compressibility, for reasons explored below.

The enhanced instability of compressible models for
$0.3\lesssim{k_{y}}W\lesssim 1$ is due to the proximity of Lindblad
resonances, as noted above. While limited to a small region of parameter
space, this result does go against the usual trend of compressibility
hindering instability.

We expect nearby Lindblad resonances to enhance RWI because the outer wave
propagation zones approach the corotation amplifier (Narayan et al., 1987;
Tsang & Lai, 2008). A similar effect is the reduction of the forbidden zone
width, i.e. Toomre $Q$ barriers, in self-gravitating disks (Mark, 1976;
Goldreich & Tremaine, 1978).

We defer a detailed study of this effect, but note the basic properties of
Lindblad resonances in our models. Their location, where
$\Delta\omega^{2}=\kappa^{2}$, is at $|x|=\pm 2/(3k_{y})$ in limit of pure
Keplerian flow in the shearing sheet (for corotation at $x=0$). This location
clearly approaches $|x|\lesssim{W}$ for $k_{y}W\gtrsim 1$. Non-Keplerian flow
affects the exact location of Lindblad resonances in the Rossby zone.

To understand why Lindblad resonances only affect compressible modes, note
that density waves only propagate where $D<0$. For Keplerian flow, this
propagation region follows from the first two terms in Eq. (9b) as
$|x|>(2/3)\sqrt{1/k_{y}^{2}+H^{2}}$, i.e. always with $|x|>2H/3$ (Artymowicz,
1993). This effective Lindblad resonance location is far from the Rossby zone
for incompressible modes with $k_{y}H\gtrsim 1$.

The simplified analyses offered in other regimes is complicated by the
presence of Lindblad singularities (where $\mathcal{F}\rightarrow\infty$) in
the Rossby zone (see Fig. 2). Lindblad singularities can be removed from the
ODE (Goldreich et al., 1986). However, they are replaced by “sonic”
singularities at $|x|\simeq 2H/3$, which also lie in the Rossby region in this
$W>H$ regime. We thus defer further analytic exploration of this regime.

Figure 8: RWI growth rate, $s$, for bumps in the compressible shearing sheet,
mapped versus scaled bump amplitude, $\Delta\varPi$, and bump width relative
to disk scale-height, $W/H$. Different values of $k_{y}H$ are shown from left
to right. Compressible (_dashed orange_) and incompressible (_yellow dotted_)
marginal stability curves are compared.

###  4.2 Compressible growth rates

Figure 8 plots the growth rates for three different levels of compressibility,
$k_{y}H=0.03,0.1,0.3$, from stronger to weaker, which also corresponds to a
range of wavelengths from long to short. The bump case is shown, but other
shapes are similar. The width is plotted as $W/H$, compared to $k_{y}W$ in
Figure 6 for a different perspective. The characteristic value of $k_{y}W=0.3$
(where compressible effects transition from stabilizing to destabilizing, as
described above) lies at $W/H=10,3,1$ (respectively) in these plots.

For fixed values of $W/H\lesssim 1$ the stability boundary is similar for the
longer wavelength (more compressible) cases $k_{y}H=0.03,0.1$ but higher for
$k_{y}H=0.3$ (pushed by the incompressible limit). However in the unstable
region, larger $k_{y}H$ modes grow faster.

For larger fixed values of $W/H\gtrsim 1$ models with smaller $k_{y}H$ have an
extended region of instability, out to $W/H\sim 1/(k_{y}H)$, but only have
faster growth (than larger $k_{y}H$ modes) close to the stability boundary.

More quantitatively, we showed that incompressible growth rates are
approximately $s_{\mathrm{inc}}\sim(\varOmega{k}_{y}^{2}\Delta\varPi)^{1/3}$.
For $k_{y}H\lesssim 0.1$, the compressible growth rates in Figure 8 are better
characterized by a reduced rate:

$\displaystyle s_{\mathrm{comp}}$
$\displaystyle\sim(\varOmega{k}_{y}^{3}H\Delta\varPi)^{1/3}\,,$ (41)

for regions not too close to the stability boundary. As with the
incompressible case, the growth rate has a weak dependence on feature width,
$W$, and increases with $k_{y}$ in the unstable region.

For quick estimates of RWI growth, Eq. (38) can be used to determine if
parameters are unstable, while the growth rate can be estimated as
$\min(s_{\rm inc},s_{\rm comp}$) for parameters a factor $\gtrsim 2$ from the
stability boundary.

##  5 Comparison to Global Models

To understand the validity and limitations of our shearing sheet models, we
compare to global, cylindrical disk models. We demonstrate that shearing box
models are valid for narrow features with sufficiently small
$W/R_{\mathrm{c}}$, and investigate the role of global disk parameters
$m=k_{y}R_{\mathrm{c}}$ and $h=H/R_{\mathrm{c}}$.

We describe our global enthalpy features, and compare to other
parameterizations in §5.1. We analyze the Rayleigh stability boundary in
global models in §5.2. We compare the RWI in shearing sheet and global models
in §5.3, and finally address RWI in dust traps in §5.4.

###  5.1 Global disk models

To best compare to our shearing sheet models, our global models use an
enthalpy feature

$\displaystyle\varPi_{0}(R)$
$\displaystyle=\left(\frac{R}{R_{\mathrm{c}}}\right)^{q}\left[\varPi_{\mathrm{b}}+{\Delta\varPi}S(\Delta{R}/W)\right]\,.$
(42)

with $\Delta{R}=R-R_{\mathrm{c}}$, powerlaw $q$, and using the same shape
functions $S$ as Eq. (13).

For the same polytropes, $P_{0}\propto\varSigma_{0}^{\Gamma}$, Eq. (18)
generalizes to

$\displaystyle\varSigma_{0}(R)$
$\displaystyle=\varSigma_{\mathrm{b}}\left(\frac{R}{R_{\mathrm{c}}}\right)^{n}\left[1+\frac{\Delta\varPi}{\varPi_{\mathrm{b}}}S(\Delta
R/W)\right]^{\tfrac{1}{\Gamma-1}}$ (43)

with $n=q/(\Gamma-1)$. We again set $\Gamma=\gamma=4/3$ and set $n=q=0$ for
simplicity in this work.

Works not using our enthalpy formulation can still be compared to our results.
Surface density features (Ono et al., 2016; Chang et al., 2023),

$\displaystyle\varSigma_{0,\varSigma}(R)$
$\displaystyle=\varSigma_{\mathrm{b}}\left(\frac{R}{R_{\mathrm{c}}}\right)^{n}\left[1+\frac{\Delta\varSigma}{\varSigma_{\mathrm{b}}}S_{\varSigma}(\Delta
R,W)\right]\,.$ (44)

are equivalent to our enthalpy formulation, with
$\varSigma_{0}\simeq\varSigma_{0,\varSigma}$, for

$\displaystyle\frac{\Delta\varSigma}{\varSigma_{\mathrm{b}}}S_{\varSigma}$
$\displaystyle\simeq\begin{cases}\dfrac{1}{\Gamma-1}\dfrac{\Delta\varPi}{\varPi_{\mathrm{b}}}S&\text{if}\,\,\dfrac{\Delta\varPi}{\varPi_{\mathrm{b}}}\ll
1,\\\
\left(\dfrac{\Delta\varPi}{\varPi_{\mathrm{b}}}S\right)^{\tfrac{1}{\Gamma-1}}&\text{if}\,\,\dfrac{\Delta\varPi}{\varPi_{\mathrm{b}}}\gg
1.\end{cases}$ (45)

For small amplitudes the shapes are the same and the amplitudes are similar.
Even the isothermal $\Gamma=1$ case follows, from
$(\Gamma-1)\varPi_{\mathrm{b}}\rightarrow{c}^{2}/\gamma$ (Eq. 19). For larger
amplitudes, neither amplitudes nor shapes are the same, so comparisons require
more care.

Instead of the above analytic comparison, for general disk features, the
corresponding enthalpy profile $\varPi_{0}(R)$ can be computed, and the
amplitude and width measured. For barotropic models, the profile is simply
$\varPi_{0}=\int{d}P_{0}/\varSigma_{0}$. For non-barotropic disks,

$\displaystyle\varPi_{0}(R)=\varPi_{0}(R_{\mathrm{r}})+\int_{R_{\mathrm{r}}}^{R}\frac{dP_{0}/dR^{\prime}}{\varSigma_{0}}dR^{\prime}\,,$
(46)

for arbitrary values of the reference disk location, $R_{\mathrm{r}}$, and
$\varPi_{0}(R_{\mathrm{r}})$. The reference enthalpy is irrelevant, as a
reference sound speed $c$ (which is not arbitrary) can be used instead. Thus,
with some effort, the enthalpy properties of any disk feature can be measured,
and applied to the results of this work.

###  5.2 Rayleigh Instability in Global Models

The results for our global models are shown in Fig. 9 for bumps and Fig. 10
for other shapes. In global models, the shearing sheet symmetry of inner vs.
outer gap edges and jumps vs. drops is broken. We show these additional cases.

We first discuss the location of the Rayleigh stability boundary. In shearing
sheet models, the Rayleigh stability boundary is at fixed $\mathcal{J}$ (Eq.
16, Figs. 3, 6, 8). For global models, the critical $\mathcal{J}$ value
changes for larger widths, $W/R_{\mathrm{c}}$.

We describe Rayleigh instability as being “enhanced” (or “reduced”), when the
critical $\mathcal{J}$ value drops (or increases) for wider features. Of the
shapes considered, most show enhanced Rayleigh instability, to quite different
degrees. Only the inner gap edge shows obviously reduced Rayleigh instability.
The difference with the outer gap edge (which has the strongest enhancement)
is striking. Curiously, the drop shape differs from the inner gap edge, which
also can be considered a drop.

We wish to understand these effects, since the Rayleigh stability boundary is
crucial for our RWI analysis. We start with the orbital frequency

$\displaystyle\varOmega_{0}^{2}(R)$
$\displaystyle=\varOmega_{\mathrm{K}}^{2}(R)+\frac{1}{R}\frac{d\varPi_{0}}{dR}$
(47)

where the Keplerian $\varOmega_{\mathrm{K}}=\varOmega_{\mathrm{Kc}}R_{\rm
c}^{3/2}/R^{3/2}$. Using Eq. (42) with $q=0$ and
$\kappa^{2}=R^{-3}d(\varOmega_{0}^{2}R^{4})/dR$ gives

$\displaystyle\frac{\kappa^{2}}{\varOmega_{\mathrm{Kc}}^{2}}$
$\displaystyle=\left(\frac{R_{\mathrm{c}}}{R}\right)^{3}+\mathcal{J}\left[S^{\prime\prime}(X)+\frac{3W}{R}S^{\prime}(X)\right]$
(48)

for $X=\Delta R/W$, which reduces to the local limit, Eq. (16), for $W,\Delta
R\ll{R}_{\mathrm{c}}$.

There are two main global effects in equation (48). The first “Keplerian”
effect is the $(R_{\mathrm{c}}/R)^{3}$ term, which enhances Rayleigh
instability if $R>R_{\mathrm{c}}$ at $\min(\kappa^{2})$. For example, the
outer edge of gaps and jumps have vorticity (and $\kappa^{2}$) minima at
$R>R_{\mathrm{c}}$. Conversely, this effect reduces Rayleigh stability for
inner gap edges and drops.

The second “non-Keplerian” effect is given by the term $3(W/R)S^{\prime}(X)$.
Sub-Keplerian speeds ($S^{\prime}(X)<0$) contribute to lower vorticity and
enhanced Rayleigh stability. This term is positive (negative) for jumps
(drops) and outer (inner) gap edges. Thus this second effect counteracts the
first for these shapes (but not bumps, as discussed last).

Which effect dominates depends on shape details, especially how far the
vorticity (and $\kappa^{2}$) minimum is from $R_{\mathrm{c}}$. For a
quantitative criterion, we Taylor expand equation (48) about $W=0$, and
evaluate at $R=R_{\mathrm{m}}=R_{\mathrm{c}}+X_{\mathrm{m}}$, the location of
the local $\kappa^{2}$ minimum. This expansion shows that $\min(\kappa^{2})$
is lower, and Rayleigh instability enhanced if

$\displaystyle f_{W}$ $\displaystyle\equiv
X_{\mathrm{m}}-\mathcal{J}S^{\prime}(X_{\mathrm{m}})=X_{\mathrm{m}}+\frac{S^{\prime}(X_{\mathrm{m}})}{S^{\prime\prime}(X_{\mathrm{m}})}>0\,.$
(49)

The final expression uses
$\mathcal{J}=\mathcal{J}_{\kappa}=-1/S^{\prime\prime}(X_{\mathrm{m}})$, the
small $W$ Rayleigh stability boundary. Since $f_{W}>0$ for outer gap edges
(the “Keplerian” effect dominates) and for drops (the “non-Keplerian” effect
dominates) the enhanced Rayleigh instability of these shapes is explained.

Shapes with $f_{W}<0$ (inner gap edges and jumps) are expected to show reduced
Rayleigh instability. This expectation holds for inner gap edges, but jumps
are more complicated. Rayleigh stability is indeed reduced as $W\rightarrow
0$, but this small effect is not visible in Figure 10. At larger
$W/R_{\mathrm{c}}$, the first order Taylor expansion is insufficient for
jumps. As $X_{\mathrm{m}}$ increases for larger $W/R_{\mathrm{c}}$, the
“Keplerian” effect dominates, explaining the enhanced Rayleigh instability
seen for jumps. The stronger competition between the two effects explains why
the jump case shows a weaker enhancement, starting at larger $W$, compared to
other shapes.

The bump case is special with $f_{W}=0$, which is marginal by our simple
criterion. However since bumps have $S^{\prime}(X)<0$ for $X>0$, both global
effects combine constructively for the bump case, unlike the other cases. Thus
bumps have enhanced Rayleigh stability at larger $W$, with $\min(\kappa^{2})$
shifting to $X>0$.

From this analysis, we come to a better understanding of the Rayleigh
stability boundary for wide disk features.

For even more extreme parameters than Rayleigh instability, outward pressure
gradients can exceed stellar gravity, giving $\varOmega_{0}^{2}<0$ (Eq. 47).
The $\min(\varOmega_{0}^{2})=0$ boundary appears in Figures 9 and Fig. 10.
This boundary is shown to emphasize how extreme this region of parameter space
is.

Figure 9: RWI marginal stability curves for bumps in a $h=0.1$ global disk
with azimuthal mode number $m=1,2,3$ (_solid blue, orange, green curves_
respectively). _Dashed curves_ show shearing sheet results for comparison. The
greyed out region is Rayleigh unstable. The _dotted yellow_ curve gives the
minimum amplitude of dust-trapping bumps.

###  5.3 Global RWI vs. Shearing Sheet

We now compare the results of our shearing sheet models to the equivalent
compressible global models, described in §5.1. Global models introduce the
lengthscale, $R_{\mathrm{c}}$ and thus one additional dimensionless parameter.
Our shearing sheet parameters $\mathcal{J},k_{y}W$ and $k_{y}H$ we add the
mode number $m=k_{y}R_{\mathrm{c}}$. Removing the local wavenumber $k_{y}$, an
equivalent set—$\mathcal{J},m,h,W/R_{\mathrm{c}}$—uses the aspect ratio
$h=H/R_{\mathrm{c}}$.

####  5.3.1 Effect of Mode Number $m$

Figure 9 plots the marginal stability curves for $m=1,2$ and 3 RWI modes in a
$h=0.1$ disk versus $W/R_{\mathrm{c}}$ for a bump shape (other shapes are
addressed next). The equivalent compressible shearing sheet models have
$k_{y}H=mh=0.1,0.2,0.3$, and their stability curves are shown for comparison.

Shearing sheet and global results agree very well for
$W/R_{\mathrm{c}}\lesssim 0.1$, as expected. This agreement is excellent even
for the most global $m=1$ modes. For $W/R_{\mathrm{c}}\gtrsim 0.1$ global and
shearing sheet results differ, moreso for lower $m$. At larger widths, global
bumps are more susceptible to the RWI than shearing sheet bumps.

In Fig. 9, the Rayleigh stability boundary deviates from constant
$\mathcal{J}$, as described in §5.2. The $m=1$ mode curves to avoid the
Rayleigh stability boundary, for $W/R_{\mathrm{c}}\gtrsim 0.1$, another
manifestation of the “halfway to Rayleigh” result.

The $m=2,3$ modes don’t share this behavior, crossing the Rayleigh boundary.
Since these modes are only weakly compressible, with $mh=k_{y}H=0.2,0.3$, the
“halfway to Rayleigh” behavior is not expected, and is also not seen in the
comparison shearing sheet models.

The shearing sheet stability boundary at $k_{y}W\simeq 0.7$ is at
$W/R_{\mathrm{c}}\simeq 0.7/m\simeq 0.7,.35,.23$ for the modes in Fig. 9. The
global models are more unstable, i.e. to larger widths than this boundary.
This destabilizing effect diminishes for smaller $W/R_{\mathrm{c}}$
boundaries, as expected.

Figure 10: RWI growth rates of $m=1$ modes in our $h=0.1$ global disk for
different shapes. Curves for marginal stability (_solid_), compressible
shearing sheet marginal stability (_dashed_), Rayleigh marginal stability
(_black dot-dashed_) and marginal dust traps (_dotted yellow_) are shown.

####  5.3.2 Effect of Feature Shape

Figure 10 plots marginal stability curves and growth rates for $m=1,h=0.1$ RWI
modes for a range of shapes (but not bumps, shown in Fig. 9). Global marginal
stability is compared to shearing sheet results for $k_{y}H=mh=0.1$. The
agreement is again good for $W/R_{\mathrm{c}}\lesssim 0.1$, i.e. $W/H\lesssim
1$ on this plot.

The inner and outer gap edges show the largest differences between global and
shearing sheet models, starting for $W/H\gtrsim 0.5$. The large distance
between the gap center and vorticity minima ($\simeq 1.7W$, see Fig. 1) is a
natural explanation.

For most shapes, global models are more susceptible to the RWI than shearing
sheet models of equivalent parameters. Inner gap edges are the only exception
(of the shapes considered). Inner gap edges are also the only shape to show
reduced Rayleigh instability at larger widths (§5.2). The “halfway to
Rayleigh” behavior of the RWI boundary thus also applies to global models as
they deviate from the shearing sheet approximation.

For narrow widths, $W/H\lesssim 1$, the growth rates in Fig. 10 are consistent
with shearing sheet (cf. Fig. 8 for $k_{y}H=0.1$), similar for all shapes, and
given approximately by Eq. (41). For wider features, we do not offer a global
correction to this analytic approximation, as the effects seem shape
dependent. The basic behavior is that growth rates steadily increase away from
the RWI boundary.

####  5.3.3 Halfway to Rayleigh, Globally

We refer back to Fig. 7 for the the minimum value of
$\kappa^{2}/\varOmega_{0}^{2}$ along the RWI boundary for global models with
$m=1,\,h=0.03$ (a slightly thinner disk than above). The results are generally
consistent with the equivalent $k_{y}H=mh=0.03$ shearing sheet models. Note
that global models compare to the orbital frequency $\varOmega_{0}(R)$, not
the fixed $\varOmega$ of the shearing sheet.

For $W\gtrsim 2H$, all models have the RWI boundary occurring “halfway to
Rayleigh” with $\min(\kappa^{2}/\varOmega_{0}^{2})\sim 0.5{-}0.6$. The inner
and outer gap edges again show the largest global corrections at larger
widths. Inner gap edges are again the most special case and the most resistant
to RWI, especially for wide gaps. Inner gap edges require the lowest values of
$\min(\kappa^{2}/\varOmega_{0}^{2})$ and the largest enthalpy amplitudes (Fig.
10) to trigger RWI.

For thinner disks, $h\lesssim 0.01$, global corrections are less significant
(for fixed $W/H$) and $\min(\kappa^{2}/\varOmega_{0}^{2})$ more constant along
the RWI boundary. This effect is due to stronger compressibility, with
$mh=k_{y}H<0.01$ (Fig. 6). Disks with moderate thickness,
$0.03\lesssim{h}\lesssim 0.3$, are more realistic but more complicated, due to
intermediate compressibility and stronger global curvature effects.

###  5.4 RWI in Dust Traps

Chang et al. (2023) examined which dust trapping rings became unstable to –
and would thus be modified by – the RWI. The condition for dust trapping is a
maximum in the midplane pressure,

$\displaystyle P_{\rm mid}$ $\displaystyle=\varOmega_{\rm
K}\frac{P_{0}}{c_{0}}=\frac{\varOmega_{\rm
K}}{\sqrt{\gamma}}\sqrt{P_{0}\varSigma_{0}}$ (50)

assuming a vertically isothermal structure.

Figures 9 and 10 show the minimum amplitude needed for dust trapping. Note
that no dust traps exist for inner gap edges or drops since they reinforce
$dP_{\rm mid}/dR<0$ instead of reversing it. Dust traps that are stable to RWI
lie in the parameter space above the yellow dotted dust trapping curves and
below the solid RWI boundaries.

As in Chang et al. (2023), which only considered bumps, stable (to RWI) dust
traps exist above a minimum width and for a range of intermediate amplitudes.
This parameter space is larger for outer gap edges and jumps, vs. bumps, for
reasons that can be explained by an analysis of $P_{\rm mid}$ similar to that
of $\kappa^{2}$ in §5.2.

We defer a more detailed study of dust trap stability that further extends
Chang et al. (2023). We mainly note that such an analysis is facilitated by
the insights to the RWI boundary established in this work.

##  6 Conclusions

We examine the linear Rossby Wave instability (RWI) with a suite of simplified
models, to gain a basic understanding of the conditions for instability, and
unstable growth rates. The disk features that trigger the RWI are best
characterized by their enthalpy amplitude, $\Delta\varPi$ and width $W$. When
different combinations of temperature and density produce the same enthalpy
profile, the equilibrium velocity and vorticity profiles are the same (Eqs. 3,
4). We apply enthalpy features with various shapes (§2.3) to a suite of
models, in the incompressible shearing sheet (§3), the compressible shearing
sheet (§4) and global models (§5). Our main insights, explored in detail in
the text, follow.

  1. 1.

The RWI in the incompressible shearing sheet (ISS) is simply characterized by
two dimensionless parameters: the scaled enthalpy amplitude
$\mathcal{J}=\Delta\varPi/(\varOmega W)^{2}$ and $k_{y}W$ (wavenumber times
width). The ISS RWI can be understood analytically, including the stability
boundary (Eqs. 27, 30) and growth rate (Eqs. 32, 34). The ISS RWI has a
similar mechanism and eigenfunctions to the full disk RWI and to generic shear
instabilities (Fig. 5).

  2. 2.

The RWI in the compressible shearing sheet (CSS) requires the additional
parameter, $k_{y}H$ (wavenumber times scale-height), an inverse Mach number.
Modes with $k_{y}H>1$ behave incompressibly. Smaller $k_{y}H$ values show
stronger compressibility effects (Fig. 6).

  3. 3.

The RWI is moderately compressible in typical protoplanetary disks with aspect
ratios $0.03\lesssim{h}\lesssim 0.3$ (Chiang & Youdin, 2010). Specifically,
$m=1$ azimuthal modes have $k_{y}H=mh\rightarrow{h}$, and are compressible,
while $m\gtrsim 1/h$ modes are incompressible.

  4. 4.

The RWI is usually most readily triggered by the longest wavelength, $m=1$
modes (Ono et al., 2016; Chang et al., 2023; §5.3.1). However, in very thin
disks, modes with different $m$ but $mh=k_{y}H\lesssim 0.01$ will have nearly
the same RWI boundary, due to strong compressibility (Eq. 38).

  5. 5.

Only disk features with widths $W\lesssim 1/k_{y}$ can trigger the RWI (Figs.
3, 6). In global models with a feature at radius $R_{\mathrm{c}}$, this limit,
$W/R_{\mathrm{c}}\lesssim 1/m$, is relevant (i.e. smaller than the disk) for
$m>1$ (Fig. 9 and §5.3.1). This limit is roughly derived in §3.1.1.

  6. 6.

The RWI boundary often lies “halfway to Rayleigh instability” in that
$\min(\kappa^{2})$ drops to $\sim 0.5{-}0.6\,\varOmega^{2}$. This behavior
occurs for widths $H\lesssim{W}\lesssim 1/k_{y}=R_{\mathrm{c}}/m$ (Fig. 7), a
range that expands for thinner disks (and thus stronger compressibility, Fig.
6). This boundary is roughly derived in §4.1.2.

  7. 7.

For narrow features, with $W\lesssim H$ the RWI boundary follows
$\Delta\varPi\propto{W}^{3}$ (Eq. 38). This scaling agrees with the low
amplitude behavior in Ono et al. (2016), see Eq. 45. We explain the relevant
factors that turn this previously known proportionality into an equality.

  8. 8.

The stability boundary for the RWI of localized disk feature (with $W\lesssim
0.2R_{\mathrm{c}}$) can be approximated by equation (38). The enthalpy
amplitude and width must be calculated (see §5.1). For wider disk features,
the “halfway to Rayleigh” criterion is a good approximation for $m=1$ modes
(Figs. 9, 10).

  9. 9.

Shape effects are generally minor when comparing the same enthalpy amplitude
and width. However bumps are the most susceptible to RWI. Wide gaps show the
largest global corrections, compared to shearing sheet models. The inner edges
of wide gaps are the least susceptible to RWI (Fig. 7, 10).

This final point implies that wide, symmetric planetary gaps the outer gap
edge should generally support more vigorous RWI and vortex formation. Vortices
at the outer edges of gaps could be more prominent in simulations and
observations for other reasons as well, including: larger area, longer orbital
and viscous timescales, numerical resolution and more dust trapping (Fu et
al., 2014; Lobo Gomes et al., 2015; Regály et al., 2017; Hammer et al., 2017).
Alternately, our results imply that over longer times, the RWI should make
wide planetary gaps more asymmetric, with closer and steeper inner edges. The
radial powerlaw of the background disk also affects gap asymmetries (Cimerman
& Rafikov, 2023).

Our simplified models neglect many physical effects, notably baroclinicity,
cooling, 3D motions and self-gravity. Previous works have studied the RWI with
these effects and shown their importance (§1). Studying these effects with
methods similar to this work—e.g. parameter studies of models of increasing
complexity—could yield further insights.

## Acknowledgements

We thank Leonardo Krapp, Juan Garrido-Deutelmoser, Gordon Ogilvie, Roman
Rafikov, Kaitlin Kratter, Wlad Lyra, Orkan Umurhan and other participants in
the sptheory (University of Arizona Star and Planet formation Theory) and
pfits+ (Planet Formation in the Southwest) group meetings for inspiring
discussions and advice. We acknowledge support from the NASA Theoretical and
Computational Astrophysical Networks (TCAN) via grant #80NSSC21K0497.

## References

  * Artymowicz (1993) Artymowicz, P. 1993, ApJ, 419, 166
  * Blumen (1970) Blumen, W. 1970, Journal of Fluid Mechanics, 40, 769, doi: 10.1017/S0022112070000435
  * Chang et al. (2023) Chang, E., Youdin, A. N., & Krapp, L. 2023, ApJ, 946, L1
  * Chiang & Youdin (2010) Chiang, E., & Youdin, A. N. 2010, Annual Review of Earth and Planetary Sciences, 38, 493
  * Cimerman & Rafikov (2023) Cimerman, N. P., & Rafikov, R. R. 2023, MNRAS, 519, 208
  * Cui et al. (2024) Cui, C., Tripathi, A., Yu, C., Lin, M.-K., & Youdin, A. 2024, arXiv e-prints, arXiv:2407.02103
  * de Val-Borro et al. (2007) de Val-Borro, M., Artymowicz, P., D’Angelo, G., & Peplinski, A. 2007, A&A, 471, 1043
  * Drazin & Reid (2004) Drazin, P. G., & Reid, W. H. 2004, Hydrodynamic Stability
  * Dra̧żkowska & Dullemond (2018) Dra̧żkowska, J., & Dullemond, C. P. 2018, A&A, 614, A62
  * Fu et al. (2014) Fu, W., Li, H., Lubow, S., & Li, S. 2014, ApJ, 788, L41
  * Fung & Ono (2021) Fung, J., & Ono, T. 2021, ApJ, 922, 13
  * Goldreich et al. (1986) Goldreich, P., Goodman, J., & Narayan, R. 1986, MNRAS, 221, 339
  * Goldreich & Lynden-Bell (1965) Goldreich, P., & Lynden-Bell, D. 1965, MNRAS, 130, 125
  * Goldreich & Tremaine (1978) Goldreich, P., & Tremaine, S. 1978, ApJ, 222, 850
  * Hammer et al. (2017) Hammer, M., Kratter, K. M., & Lin, M.-K. 2017, MNRAS, 466, 3533
  * Heifetz et al. (1999) Heifetz, E., Bishop, C. H., & Alpert, P. 1999, Quarterly Journal of the Royal Meteorological Society, 125, 2835, doi: 10.1002/qj.49712556004
  * Hu et al. (2018) Hu, X., Tan, J. C., Zhu, Z., et al. 2018, ApJ, 857, 20
  * Huang & Yu (2022) Huang, S., & Yu, C. 2022, MNRAS, 514, 1733
  * Johansen et al. (2009) Johansen, A., Youdin, A., & Klahr, H. 2009, ApJ, 697, 1269
  * Johnson & Gammie (2005) Johnson, B. M., & Gammie, C. F. 2005, ApJ, 626, 978
  * Klahr et al. (2023) Klahr, H., Baehr, H., & Melon Fuksman, J. D. 2023, arXiv e-prints, arXiv:2305.08165
  * Kuznetsova et al. (2022) Kuznetsova, A., Bae, J., Hartmann, L., & Mac Low, M.-M. 2022, ApJ, 928, 92
  * Latter & Papaloizou (2017) Latter, H. N., & Papaloizou, J. 2017, MNRAS, 472, 1432
  * Lesur et al. (2023) Lesur, G., Flock, M., Ercolano, B., et al. 2023, in Astronomical Society of the Pacific Conference Series, Vol. 534, Protostars and Planets VII, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 465
  * Li et al. (2001) Li, H., Colgate, S. A., Wendroff, B., & Liska, R. 2001, ApJ, 551, 874
  * Li et al. (2000) Li, H., Finn, J. M., Lovelace, R. V. E., & Colgate, S. A. 2000, ApJ, 533, 1023
  * Lin (2013) Lin, M.-K. 2013, ApJ, 765, 84
  * Lin & Papaloizou (2010) Lin, M.-K., & Papaloizou, J. C. B. 2010, MNRAS, 405, 1473
  * Lin & Youdin (2015) Lin, M.-K., & Youdin, A. N. 2015, ApJ, 811, 17
  * Lithwick (2007) Lithwick, Y. 2007
  * Liu & Bai (2023) Liu, H., & Bai, X.-N. 2023, MNRAS, 526, 80
  * Lobo Gomes et al. (2015) Lobo Gomes, A., Klahr, H., Uribe, A. L., Pinilla, P., & Surville, C. 2015, ApJ, 810, 94
  * Lovelace et al. (1999) Lovelace, R. V. E., Li, H., Colgate, S. A., & Nelson, A. F. 1999, ApJ, 513, 805
  * Lyra et al. (2024) Lyra, W., Yang, C.-C., Simon, J. B., Umurhan, O. M., & Youdin, A. N. 2024, arXiv e-prints, arXiv:2406.17934
  * Manasreh (2012) Manasreh, M. O. 2012, Introduction to nanomaterials and devices, 1st edn. (Hoboken, N.J: Wiley)
  * Manger et al. (2020) Manger, N., Klahr, H., Kley, W., & Flock, M. 2020, MNRAS, 499, 1841
  * Mark (1976) Mark, J. W. K. 1976, ApJ, 205, 363
  * Meheut et al. (2012) Meheut, H., Yu, C., & Lai, D. 2012, MNRAS, 422, 2399
  * Narayan et al. (1987) Narayan, R., Goldreich, P., & Goodman, J. 1987, MNRAS, 228, 1
  * Nelson et al. (2013) Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610
  * Ono et al. (2016) Ono, T., Muto, T., Takeuchi, T., & Nomura, H. 2016, ApJ, 823, 84
  * Pinilla & Youdin (2017) Pinilla, P., & Youdin, A. 2017, in Astrophysics and Space Science Library, Vol. 445, Astrophysics and Space Science Library, ed. M. Pessah & O. Gressel, 91
  * Regály et al. (2017) Regály, Z., Juhász, A., & Nehéz, D. 2017, ApJ, 851, 89
  * Tsang & Lai (2008) Tsang, D., & Lai, D. 2008, MNRAS, 387, 446
  * Umurhan et al. (2016) Umurhan, O. M., Shariff, K., & Cuzzi, J. N. 2016, ApJ, 830, 95
  * Vanon & Ogilvie (2016) Vanon, R., & Ogilvie, G. I. 2016, MNRAS, 463, 3725