# Emergence of vortices at the edges of planet-driven gaps in protoplanetary
discs

Nicolas P. Cimerman1, Roman R. Rafikov1,2
1Department of Applied Mathematics and Theoretical Physics, University of
Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK
2Institute for Advanced Study, Einstein Drive, Princeton, NJ 08540, USA
E-mail<EMAIL_ADDRESS>(NPC)

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Young planets embedded in protoplanetary discs (PPDs) excite spiral density
waves, which propagate, shock and deposit angular momentum in the disc. This
results in gap opening around the planetary orbit, even for low (sub-thermal)
mass planets, provided that the effective viscosity in the disc is low. The
edges of these planet-induced gaps are known to be prone to emergence of
observable vortices via the Rossby Wave Instability (RWI). We study timescales
for the development of vortices driven by low mass planets in inviscid discs.
We employ a recently developed semi-analytical theory of vortensity production
by the planet-driven shock to predict vortensity evolution near the planet,
from which we derive the radial profile of the planet-induced gap as a
function of time (this procedure can have multiple other uses, e.g. to study
dust trapping, suppression of pebble accretion, etc.). We then analyze the
linear stability of the gap edges against the RWI, obtaining the timescales
for the first appearance of unstable modes and (later) fully developed
vortices at gap edges. We present useful formulae for these timescales as
functions of planetary and disc parameters and provide their physical
justification. We also thoroughly test our semi-analytical framework against
high resolution 2D hydrodynamic simulations, confirming the accuracy of our
theoretical predictions. We discuss ways in which our semi-analytical
framework can be extended to incorporate additional physics, e.g. planetary
accretion, migration, and non-zero disc viscosity. Our results can be used to
interpret observations of PPDs and to predict emergence of vortices in
simulations.

###### keywords:

hydrodynamics – instabilities – shock waves – accretion discs – planets and
satellites: formation – methods: numerical

††pubyear: 2020††pagerange: Emergence of vortices at the edges of planet-
driven gaps in protoplanetary discs–E

##  1 Introduction

Recent observations of protoplanetary discs (hereafter PPDs) have shown a
plethora of substructures in the millimetre continuum emission that probes the
spacial distribution of large dust grains (e.g Andrews, 2020). While many
substructures are axisymmetric, ring- or gap-like, there are some that show
intriguing non-axisymmetric lobes or arcs (van der Marel et al., 2016; Kraus
et al., 2017; Dong et al., 2018; Pérez et al., 2018), which have been
interpreted as dust traps inside vortices (e.g. van der Marel et al., 2013;
Baruteau et al., 2019). Supporting this interpretation, velocity measurements
using CO emission lines have revealed kinematic structures in HD 142527, which
have been tentatively attributed to vortices (Boehler et al., 2021). Providing
further information, vortices can also reveal themselves in near-infrared
scattered light images of discs (Zhu et al., 2015; Marr & Dong, 2022). Beyond
these promising observable features, vortices have also been suggested to
provide efficient particle traps that may gravitationally collapse and enhance
planetesimal formation rates (Meheut et al., 2012; Zhu et al., 2014; Zhu &
Stone, 2014). They can also affect orbital migration of planets via their
mutual gravitational coupling (e.g Li et al., 2009; Yu et al., 2010; Lin &
Papaloizou, 2010; McNally et al., 2019). These observations and theoretical
ideas provide ample motivation to study the origin of vortices in PPDs,
especially since they can act as signposts of ongoing planet formation.

Indeed, young planets that are still embedded in their natal PPDs are known to
gravitationally excite spiral density waves that carry energy and angular
momentum across the disc. Wave angular momentum can be transferred to the
background disc material via the linear or non-linear wave damping, driving
disc evolution and resulting in gap formation around the planetary orbit
(Goodman & Rafikov, 2001; Rafikov, 2002a, 2016; Miranda & Rafikov, 2020a, b).
The edges of deep planet-driven gaps can become susceptible to the Rossby Wave
Instability (RWI, Lovelace et al., 1999; Li et al., 2000), a linear shear
instability associated with growing non-axisymmetric perturbations, that
eventually breaks the gap edge into vortices (de Val-Borro et al., 2007; Li et
al., 2009; Yu et al., 2010; Lin & Papaloizou, 2010). Thus, vortices can
potentially reveal the planets that induce them in the first place.

The details of the vortex production process depend on the disc and planet
properties, as well as the wave damping mechanism. The efficiency of linear
wave damping mechanisms depends on disc thermodynamics and effective viscosity
(Goodman & Rafikov, 2001; Miranda & Rafikov, 2020a, b), while non-linear
damping — wave steepening due to its non-linearity (a finite amplitude effect)
resulting in the formation of a shock — is robust and unavoidable. Partly for
that reason, in this work we will focus on the non-linear wave damping.

This type of damping is very efficient even for weakly non-linear waves
launched by relatively low-mass planets (Goodman & Rafikov, 2001; Rafikov,
2002a), with masses $M_{\mathrm{p}}$ below the so-called thermal mass:

$\displaystyle
M_{\mathrm{th}}=\frac{c_{\mathrm{s}}^{3}}{\Omega_{\mathrm{p}}G}=\left(\frac{H_{\mathrm{p}}}{R_{\mathrm{p}}}\right)^{3}M_{\star}=h_{\mathrm{p}}^{3}\,M_{\star},$
(1)

where $c_{\mathrm{s}}$ is the sound speed, $\Omega_{\mathrm{p}}$ is the
orbital angular frequency at the planetary distance $R_{\mathrm{p}}$,
$H_{\mathrm{p}}$ and $h_{\mathrm{p}}=H_{\mathrm{p}}/R_{\mathrm{p}}$ are the
disc scale height and aspect ratio at $R_{\mathrm{p}}$, and $M_{\star}$ is the
stellar mass. It has been shown that even such sub-thermal planets can be
responsible for producing prominent features in protoplanetary discs, both
gaps/rings (Dong et al., 2017; Bae et al., 2017; Zhang et al., 2018; Miranda &
Rafikov, 2019b, 2020a, 2020b) and vortices (Hallam & Paardekooper, 2020;
Hammer et al., 2021). For this reason, as well as to enable (semi-)analytical
progress in understanding vortex formation, in this work we will focus on
planets with $M_{\mathrm{p}}\lesssim M_{\mathrm{th}}$.

Non-linear damping of planet-driven density waves proceeds via their
inevitable evolution into shocks. Planetary shock fronts are known (Kevlahan,
1997; Lin & Papaloizou, 2010; Dong et al., 2011b) to produce a jump of
vortensity (or potential vorticity)
$\zeta=\Sigma^{-1}\mathbf{\nabla}\times\mathbf{u}$ (here $\Sigma$ is the disc
surface density and $\mathbf{u}$ is the fluid velocity), a quantity that is
otherwise conserved along streamlines in two-dimensional (2D), inviscid and
barotropic flows. As fluid elements periodically cross the planet-induced
shock, disc vortensity near the planet steadily evolves, a process closely
related to gap opening (Muto et al., 2010). Eventually, this vortensity
evolution triggers the RWI (Lovelace et al., 1999) and vortices emerge in the
disc.

The timescale on which this pathway for vortex formation operates is an
important diagnostic, which can help us constrain the properties of putative
planets triggering vortices. Indeed, higher mass planets launch stronger
density waves, which drive faster vortensity evolution and cause vortices to
form earlier. Thus, if a PPD is observed to have a gap suggestive of a planet
opening it and features a vortex at the gap edge, this would imply that either
(i) planetary mass is high or (ii) the planet has been present for a long
time. Both constraints are useful for interpreting observations. However, so
far the problem of determining the vortex emergence timescale has received
little attention111Hallam & Paardekooper (2020) studied the relation between
the planetary growth timescale and the emergence of vortices, but they also
varied disc viscosity and adopted a particular accretion history of the
planet., with most studies focusing on the long-term survival of already well-
developed vortices under various secondary instabilities (Lesur & Papaloizou,
2009), thermal relaxation (Fung & Ono, 2021; Rometsch et al., 2021) and other
processes.

The main goal of our present study is to determine the critical time for the
appearance of planet-driven vortices in inviscid discs and to explore the
dependence of this timescale on the planetary and disc parameters, which are
fixed in time. We do this using both direct hydrodynamic simulations and semi-
analytical theory based on the recent work of Cimerman & Rafikov (2021,
hereafter Paper I), in which we developed from first principles (and verified
numerically) a semi-analytical framework for calculating the vortensity jump
$\Delta\zeta$ across a planet-driven spiral shock in a barotropic disc, taking
disc and planet parameters as inputs. This framework as well as the linear
stability analysis of the planet-induced gaps against the RWI are the key
theoretical foundations of our present study. We also perform a thorough test
of our method against direct, high resolution, 2D, hydrodynamic simulations
for the range of relevant parameters, allowing us to verify its validity and
identify its limitations.

This paper is organised as follows (readers who are only interested in results
for the timescales and their applicability might skip to Section 6). We
describe the problem setup in Section 2, followed by a description of non-
linear simulations and their typical outcomes in Section 3. In Section 4 we
introduce the linear stability analysis for RWI at gap edges, explain our
method for obtaining the disc profile from its vortensity structure and show
typical results for unstable modes. Relevant timescales for different stages
of developing RWI are introduced and motivated, together with related
diagnostics in Section 5. In Section 6 we show results for these timescales
and provide power law fits, propose a heuristic theoretical explanation for
these dependencies and study the vortensity levels at the onset of RWI. We
discuss our findings, their applications and limitations in Section 7 before
summarizing in Section 8. Appendices contain further information and tests
regarding our methods.

##  2 Problem setup

Our setup for studying planet-disc interaction is the same as in Paper I.
Here, we recite its most important aspects.

We consider a planet of mass $M_{\mathrm{p}}$ orbiting a star of mass
$M_{\star}$ on a circular orbit with a semi-major axis $R_{\mathrm{p}}$ that
lies within a thin two-dimensional gas disc. We adopt polar coordinates
$(R,\phi)$. The initial (background) disc state, unperturbed by the planet has
a power law profile of the surface density

$\displaystyle\Sigma_{\mathrm{i}}(R)=\Sigma_{\mathrm{i}}(R_{\mathrm{p}})\left(\frac{R}{R_{\mathrm{p}}}\right)^{-p}=\Sigma_{\mathrm{p}}\left(\frac{R}{R_{\mathrm{p}}}\right)^{-p},$
(2)

where $p$ is a constant.222We have changed the subscript for the initial
conditions from ’0’ to ’i’ as compared to Paper I, to avoid confusion with
azimuthally averaged variables, i.e. their $m=0$ Fourier components. We assume
that the mass of the disc is small, $M_{\mathrm{d}}\ll M_{\star}$, such that
its self-gravity can be neglected.

We use a globally isothermal equation of state (EoS),
$P=c_{\mathrm{s}}^{2}\Sigma$, with the spacially constant speed of sound
$c_{\mathrm{s}}$, which is barotropic (i.e. $P=P(\Sigma)$ only). The choice of
this EoS instead of an often used non-barotropic locally isothermal EoS, for
which the sound speed follows a prescribed radial profile $c_{\mathrm{s}}(R)$,
is motivated in more detail in Section 7.2.1 and Paper I. The 2D approximation
is appropriate for thin discs, such that the aspect-ratio of the disc
$h=H/R=c_{\mathrm{s}}/(\Omega R)\ll 1$.

Very importantly, in this study (and unlike e.g. Hammer et al., 2017; McNally
et al., 2019; Hallam & Paardekooper, 2020) we keep both planetary
($M_{\mathrm{p}}$ and $R_{\mathrm{p}}$) and disc ($h_{\mathrm{p}}$) properties
fixed in time to highlight the key physical processes and to reduce the number
of relevant parameters. In Section 7.1 we comment on how this assumption can
be relaxed.

In general, the disc maintains a radial centrifugal balance accounting for the
radial pressure gradient:

$\displaystyle\Omega^{2}(R)=\Omega_{\mathrm{K}}^{2}(R)+\frac{c_{\mathrm{s}}^{2}}{R\Sigma(R)}\frac{\mathrm{d}\Sigma(R)}{\mathrm{d}R},$
(3)

where $\Omega_{\mathrm{K}}=\sqrt{GM_{\star}/R^{3}}$ is the Keplerian orbital
frequency. The initial, unperturbed (by either a planet or vortices) disc has
radial velocity $u_{\mathrm{i},R}(R)=0$ and azimuthal velocity
$u_{\mathrm{i},\phi}(R)=R\Omega_{\mathrm{i}}$, where $\Omega_{\mathrm{i}}(R)$
is given by equation (3) with $\Sigma(R)=\Sigma_{\mathrm{i}}(R)$.

Our fiducial disc model is the same as in Paper I. It has an aspect-ratio
$h_{\mathrm{p}}=0.05$ at $R_{\mathrm{p}}$ and a surface density slope $p=3/2$.
The latter results in an initial vortensity profile $\zeta_{\mathrm{i}}(R)$
that is almost constant for slightly sub-Keplerian discs, see Section 6 of
Paper I.

When presenting our results we adopt units where
$G=M_{\ast}=\Omega_{K}(R_{\mathrm{p}})=\Sigma_{\mathrm{p}}=1$. In these units
the planetary orbital time is $P_{\mathrm{p}}=2\pi$, which we use as a unit of
time.

We now describe the different methods, both numerical and analytical, used in
this work (their usage is illustrated in Fig. 7). We use the fiducial setup
described earlier to showcase our methods and typical results.

##  3 Methods: hydrodynamical simulations

Figure 1: Simulation results using fiducial disc parameters and an
intermediate-mass planet $M_{\mathrm{p}}/M_{\mathrm{th}}=0.25$. All
perturbations are w.r.t. the initial conditions
($\Sigma_{\mathrm{i}},\zeta_{\mathrm{i}}$). (a) Map of the relative surface
density perturbation. (b) Radial profile of the azimuthally averaged surface
density perturbation. (c) Radial profile of the azimuthally averaged
vortensity perturbation. (d) Map of the vortensity perturbation. The dashed
horizontal lines show the predicted shock locations $R=R_{\mathrm{p}}\pm
l_{\mathrm{sh}}$. At the outer gap edge, the RWI has developed and led to the
formation of three large vortices. These cause additional, weaker spiral waves
that are visible in panel (a). 2D maps and blue solid lines in panels (b) and
(c) show the disc state at $t=780P_{\mathrm{p}}$. Additional dot-dashed lines
in the middle panels correspond to $t=\\{500,640,920\\}P_{\mathrm{p}}$.
Figure 2:  Two-dimensional map of the non-axisymmetric part of the vortensity
perturbation, $\zeta-\left\langle\zeta\right\rangle_{\phi}$, for the fiducial
case. Panels show results for different times (labels in panels) during the
onset of instability. The shock front is visible as a thin curve. In panel (a)
the flow is stable and mostly axisymmetric. At the predicted distance
$l_{\mathrm{sh}}$ from the planet (black dashed lines), the shock starts
leading to roughly axisymmetric vortensity distribution — two thin rings
around $R\simeq 0.9$ and $R\simeq 1.1$ (locally perturbed around the shock).
Other panels show (b) the onset of RWI in its initial stage, (c) the
development of three vortices and (d) a single vortex as their merger product
at the outer gap edge. At the inner gap edge, we observe similar processes
starting later.

###  3.1 Setup and overview of simulations

To study vortex formation we perform global, non-linear hydrodynamic
simulations of planet-disc interaction using Athena++ (Stone et al.,
2020)333Athena++ is publicly available on GitHub.. The code uses a Godunov
scheme to solve the hydrodynamic equations in a conservative form:

$\displaystyle\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\mathbf{u})$
$\displaystyle=0,$ (4) $\displaystyle\frac{\partial(\rho\mathbf{u})}{\partial
t}+\nabla\cdot(\rho\mathbf{u}\otimes\mathbf{u}+P\mathbf{I})$
$\displaystyle=-\rho\nabla\Phi,$ (5)

where $P$ is the gas pressure, $\mathbf{I}$ the identity tensor, and $\Phi$ is
the total gravitational potential. Following Paper I, we employ a fourth-order
smoothed potential for the planet $\Phi_{\mathrm{p}}=\Phi_{\mathrm{p}}^{(4)}$
(Dong et al., 2011a).

The simulation domain extends over the radial range $0.2\leq
R/R_{\mathrm{p}}<4.0$, with logarithmic grid spacing. The full azimuthal range
$0\leq\phi\leq 2\pi$ is covered uniformly. In this work, we use a fixed
resolution of $N_{R}\times N_{\phi}=3448\times 7200$, corresponding to 58
cells per disc scale-height at $R_{\mathrm{p}}$ for the fiducial disc with
$h_{\mathrm{p}}=0.05$. For hotter discs (larger $h_{\mathrm{p}}$) this results
in higher effective resolution. We use Roe’s approximate Riemann solver,
linear spacial interpolation and second order accurate time-stepping. We
employ the orbital advection algorithm, which has recently been implemented in
Athena++. For more details regarding the numerical setup, see Paper I.

###  3.2 Typical simulation outcomes

We now describe typical outcomes of our simulations. Our fiducial model, which
we will use as a reference, has $M_{\mathrm{p}}=0.25M_{\mathrm{th}}$,
$h_{\mathrm{p}}=0.05$ and $p=1.5$ (i.e. $\zeta_{\mathrm{i}}$ close to a
constant). We illustrate our results with 2D snapshots of the surface density
and vortensity perturbations at a particular moment of time,
$t=780P_{\mathrm{p}}$, in panels (a) and (d) of Fig. 1. In panels (b) and (c)
of that figure we also show the azimuthally-averaged profiles of the $\Sigma$
and $\zeta$ perturbations at several moments of time to illustrate their
steady evolution. Additionally, we display maps of
$\zeta-\left\langle\zeta\right\rangle_{\phi}$, i.e. the non-axisymmetric part
of the vortensity, in Fig. 2, which makes it easier to identify vortices as
they develop.

For more than 700 orbits of the simulation, the vortensity perturbation,
driven by the spiral shocks, increases (close to) linearly with time,
maintaining an almost constant
$\partial_{t}\left\langle\zeta\right\rangle_{\phi}(R)$ and showing no sign of
saturation (see also Fig. 12 in Paper I). As expected, vortensity stays almost
unperturbed in the radial band $|R-R_{\mathrm{p}}|<l_{\mathrm{sh}}$ (shown
with horizontal dashed lines), where (Goodman & Rafikov, 2001)

$\displaystyle l_{\mathrm{sh}}\approx
0.86\,h_{\mathrm{p}}\,R_{\mathrm{p}}\left(M_{\mathrm{p}}/M_{\mathrm{th}}\right)^{-2/5}$
(6)

(for isothermal EoS) is the wave shocking length; after travelling this
distance the wave shocks, driving time evolution of $\zeta$ for
$|R-R_{\mathrm{p}}|>l_{\mathrm{sh}}$. This evolution, visible in Fig. 1c, is
associated with the formation of a gap that becomes deeper and deeper (see
Fig. 1b), leading to increasing radial gradients of $\Sigma$ and pressure.
This, in turn, leads to increased radial shear of the azimuthal flow, which
nevertheless remains stable (laminar) during this stage. This can be seen in
Fig. 2a, in which the only non-axisymmetric vortensity perturbation is
associated with the planetary spiral shock.

This picture begins to change at around $t=760P_{\mathrm{p}}$, when our
simulation starts exhibiting RWI at the outer gap edge
($R/R_{\mathrm{p}}\simeq 1.1$), which is clear from the non-axisymmetric
vortensity structures appearing in Fig. 2b. As the instability sets in, non-
axisymmetric features also become apparent in the surface density and
vortensity perturbations in Fig. 1a,d slightly later, at
$t=780P_{\mathrm{p}}$; the corresponding $m=3$ vortensity perturbation at the
outer gap edge is also clear in Fig. 2c, appearing as three vortices. These
vortices drive additional spiral density waves, visible in the inner and outer
disc, that are weaker than the planet-driven density waves (see also Fig. 1 in
Paper I, when vortices have not formed yet).

At $t=830$ orbits, the multiple vortices have coalesced to form one single
$m=1$ vortex. At the inner gap edge ($R/R_{\mathrm{p}}\simeq 0.9$), the RWI
sets in at around $t=880P_{\mathrm{p}}$, resulting in merging vortices. This
state is illustrated in Fig. 2d. By $t=1000P_{\mathrm{p}}$, both gap edges
feature a single vortex. The instability smooths out radial gradients of
vortensity (and related fluid variables), but does not eliminate their extrema
(see also Meheut et al., 2010).

##  4 Methods: linear analysis of RWI development

In addition to numerical simulations capturing the development of the RWI
directly, we also studied the instability using semi-analytical means. Below
we present our setup for linear stability analysis of an axisymmetric, 2D
barotropic and inviscid flow.

###  4.1 Linear RWI stability analysis

Consider a barotropic flow in an axisymmetric 2D disc in radial equilibrium
with $u_{0,R}=0$ everywhere, with prescribed radial profiles of surface
density and azimuthal velocity, $\Sigma_{0}(R)$ and
$u_{0,\phi}(R)=R\Omega_{0}(R)$ 444We note that this axisymmetric state is in
general different from the initial conditions we use for disc-planet
interactions, which are denoted $\Sigma_{0}$, $\Omega_{0}$ in Paper I and
$\Sigma_{\mathrm{i}}$, $\Omega_{\mathrm{i}}$ in this work..

Introducing small, non-axisymmetric perturbations to
${X\in[\Sigma,u_{R},u_{\phi}]}$ of the form 555The perturbations $\delta
X=X-X_{0}$ are perturbations to an axisymmetric state taken at some time $t$
after introduction of the planet and are different from the perturbations with
respect to the initial state (i.e. $X-X_{\mathrm{i}})$.

$\displaystyle\delta X(R,\phi,\omega,m,t)=\delta
X_{m}(R)\exp\left[\mathrm{i}(m\phi-\omega t)\right],$ (7)

($\left|{\delta X}\right|\ll\left|{X_{0}}\right|$) where $\delta X_{m}(R)$ is
a complex amplitude, $m$ is the azimuthal mode number and

$\displaystyle\omega=\omega_{\mathrm{R}}+\mathrm{i}\gamma$ (8)

is the complex mode frequency, with real part $\omega_{\mathrm{R}}$ and
imaginary part $\gamma$ (unstable modes have $\gamma>0$ and grow exponentially
with time).

As shown in Lovelace et al. (1999), the linearised perturbation equations
(ignoring the planetary potential) can be combined into a master equation

$\displaystyle\Psi_{m}^{\prime\prime}+B(R,\omega,m)\Psi_{m}^{\prime}+C(R,\omega,m)\Psi_{m}=0,$
(9)

where primes stand for radial derivatives and

$\displaystyle\Psi_{m}\equiv\frac{\delta P_{m}}{\Sigma_{0}},$ (10)

is the enthalpy perturbation ($\Psi_{m}$ is its Fourier component), equal to
$c_{\mathrm{s}}^{2}\delta\Sigma_{m}/\Sigma_{0}$ in the barotropic case
considered here. The coefficients in equation (9) are related to the disc
structure and are given by

$\displaystyle B(R,\omega,m)$
$\displaystyle\equiv\frac{1}{R}+\frac{\mathcal{F}^{\prime}}{\mathcal{F}}-\frac{\Omega_{0}^{\prime}}{\Omega_{0}},$
(11) $\displaystyle C(R,\omega,m)$
$\displaystyle\equiv-\frac{m^{2}}{R^{2}}-\frac{\kappa_{0}^{2}-\Delta\omega^{2}}{c_{0}^{2}}-2\frac{m}{R}\frac{\Omega_{0}}{\Delta\omega}\frac{\mathcal{F}^{\prime}}{\mathcal{F}},$
(12) $\displaystyle\mathcal{F}(R,\omega,m)$
$\displaystyle\equiv\frac{\Sigma_{0}\Omega_{0}}{\kappa_{0}^{2}-\Delta\omega^{2}}.$
(13)

Here $\Delta\omega\equiv\omega-m\Omega_{0}$ is the Doppler-shifted mode
frequency and

$\displaystyle\kappa_{0}^{2}=\frac{2\Omega_{0}}{R}\frac{\mathrm{d}(R^{2}\Omega_{0})}{\mathrm{d}R},$
(14)

is the square of the epicyclic frequency. Note that in an axisymmetric 2D disc
in radial force equilibrium (with $u_{R}=0$ everywhere) vortensity can be
written as

$\displaystyle\zeta_{0}(R)=\frac{1}{R\Sigma_{0}}\frac{\mathrm{d}(R^{2}\Omega_{0})}{\mathrm{d}R}=\frac{\kappa^{2}_{0}}{2\Sigma_{0}\Omega_{0}}.$
(15)

Our strategy for numerically solving equation (9) is described in Appendix A.

Extremal values in the radial vortensity profile are a necessary condition for
the development of the RWI in barotropic discs (Lovelace et al., 1999), and
Papaloizou & Lin (1989) argued that it is the vortensity minima that become
unstable. The edges of planet-induced gaps in protoplanetary discs can
naturally provide such extrema that may become linearly unstable (e.g. de Val-
Borro et al., 2007; Li et al., 2009). We can expose the background vortensity
$\zeta_{0}=\kappa_{0}^{2}/(2\Sigma_{0}\Omega_{0})$ in the coefficients
(11)-(13) by noting that

$\displaystyle\mathcal{F}^{-1}$
$\displaystyle=\frac{\kappa_{0}^{2}-\Delta\omega^{2}}{\Sigma_{0}\Omega_{0}}=2\zeta_{0}-\frac{\Delta\omega^{2}}{\Sigma_{0}\Omega_{0}}.$
(16)

For corotating modes with $|\Delta\omega|^{2}\ll|\kappa_{0}|^{2}$, one has
$\mathcal{F}^{-1}\simeq 2\zeta_{0}$ (Lovelace et al., 1999).

Defining a new function

$\displaystyle\Xi(R)\equiv\sqrt{\frac{R\mathcal{F}}{\Omega_{0}}}\Psi,$ (17)

one can also transform equation (9) into the form

$\displaystyle\Xi_{m}^{\prime\prime}-D(R,\omega,m)\Xi_{m}=0$ (18)

resembling the time-independent Schrödinger equation with complex $\Xi_{m}$
and

$\displaystyle D(R,\omega,m)\equiv\frac{B^{\prime}}{2}+\frac{B^{2}}{4}-C$ (19)

playing the role of the potential (Ono et al., 2016). Regions of
$\mathrm{Re}(D)<0$ are classically allowed for wave propagation, while regions
of $\mathrm{Re}(D)>0$ are classically forbidden regions, through which waves
can only tunnel. Sufficiently deep and wide troughs in $\mathrm{Re}(D)$ allow
for trapped modes.

Due to our barotropic equation of state, the Brunt-Väisälä frequency
associated with buoyancy vanishes ($N^{2}=0$), such that stability against
axisymmetric perturbations is dictated by the Rayleigh criterion (e.g.
Chandrasekhar, 1961), i.e. $\kappa_{0}^{2}<0$ is needed for instability. On
the other hand, equation (15) makes it clear that $\kappa_{0}^{2}<0$ requires
$\zeta_{0}<0$. We have checked that all disc models studied in this work
remain stable against axisymmetric perturbations (i.e. the minima of
vortensity never become negative), according to this criterion. In agreement
with previous works, we find that the RWI always occurs before the Rayleigh
criterion is violated (Li et al., 2000; Les & Lin, 2015).

Equations (9)-(19) make it clear that the linear analysis of RWI requires the
knowledge of radial profiles of $\Sigma_{0}$, $\Omega_{0}$. These could be
derived from simulations, but for a fully self-contained semi-analytical
analysis one would like to obtain $\Sigma_{0}(R)$, $\Omega_{0}(R)$ without
simulations. Next we propose an alternative semi-analytical procedure for
constructing $\Sigma_{0}(R)$, $\Omega_{0}(R)$ profiles.

For the rest of this work, we will drop the subscript on $\Sigma_{0}(R)$,
$\Omega_{0}(R)$ for brevity, with implicit understanding that $\Sigma(R)$,
$\Omega(R)$, etc. refer to the current azimuthally-averaged disc
characteristics.

###  4.2 Reconstruction of surface density and rotation profiles

In Paper I we developed a semi-analytical procedure for predicting the
planetary shock-driven evolution of the disc vortensity $\zeta(R,t)$ in time
and space. Now we describe how to use $\zeta(R,t)$ obtained using this
procedure and the method outlined in Lin & Papaloizou (2010) to reconstruct
the surface density and rotation profiles in the disc as a function of time.

####  4.2.1 Theoretical considerations

Knowing the vortensity jump at the planet-driven shock $\Delta\zeta(R)$ (Paper
I), we can calculate the associated rate of change of the disc vortensity
$S_{\zeta,\mathrm{sh}}$ according to the following formula:

$\displaystyle\left.\frac{\partial\zeta(R,t)}{\partial
t}\right|_{\mathrm{sh}}\equiv
S_{\zeta,\mathrm{sh}}=\Delta\zeta(R)\times\frac{|\Omega(R)-\Omega_{\mathrm{p}}|}{2\pi}.$
(20)

Assuming $\Omega(R)\simeq\Omega_{\mathrm{K}}$, approximating
$S_{\zeta,\mathrm{sh}}$ as constant in time (see Paper I, regarding the
validity of this assumption), and integrating over time we obtain the
vortensity profile at time $t$ after the introduction of the planetary
potential as

$\displaystyle\zeta(R,t)=\zeta(R,t=0)+tS_{\zeta,\mathrm{sh}}(R).$ (21)

Substituting $\Omega$ from the relation (3) into Equation (15) and assuming a
barotropic EoS one obtains

$\displaystyle\zeta=\frac{\Omega_{\mathrm{K}}^{2}+\frac{3c_{\mathrm{s}}^{2}}{R}\frac{\mathrm{d}\ln\Sigma}{\mathrm{d}R}+c_{\mathrm{s}}^{2}\frac{\mathrm{d}^{2}\ln\Sigma}{\mathrm{d}R^{2}}}{2\Sigma\left(\Omega_{\mathrm{K}}^{2}+\frac{c_{\mathrm{s}}^{2}}{R}\frac{\mathrm{d}\ln\Sigma}{\mathrm{d}R}\right)^{1/2}}.$
(22)

Following Lin & Papaloizou (2010), but using our globally isothermal EoS,
$c_{\mathrm{s}}=$ const, instead of their locally isothermal EoS, we can
rearrange this into a second order (non-linear) differential equation for
$\Sigma$ similar to their equation (18):

$\displaystyle\frac{1}{R^{3}}\frac{\mathrm{d}}{\mathrm{d}R}\left(R^{3}\frac{\mathrm{d}\ln\Sigma}{\mathrm{d}R}\right)+\frac{\Omega_{\mathrm{K}}^{2}}{c_{\mathrm{s}}^{2}}=\frac{2\Sigma}{c_{\mathrm{s}}^{2}}\zeta\left(\Omega_{\mathrm{K}}^{2}+\frac{c_{\mathrm{s}}^{2}}{R}\frac{\mathrm{d}\ln\Sigma}{\mathrm{d}R}\right)^{1/2},$
(23)

where vortensity enters only as a source term on the right-hand side.

This differential equation can be solved numerically in $R$, with the input
vortensity profile $\zeta(R,t)$ given by equation (21) at any $t$. To this
end, we use an iterative relaxation scheme and impose the unperturbed surface
density ($\Sigma=\Sigma_{\mathrm{i}}$) as a boundary condition far from the
planet (see Appendix B for more details).

Having found $\Sigma(R)$, we then retrieve the rotation profile $\Omega(R)$
using equation (3), thus fully determining the axisymmetric disc state (i.e.
$\zeta(R)\rightarrow[\Sigma(R),\Omega(R)])$. According to Paper I, our recipe
for constructing $\zeta$ is accurate for $0.05\leq
M_{\mathrm{p}}/M_{\mathrm{th}}<1$, implying that the method for constructing
time-dependent $\Sigma(R),\Omega(R)$ profiles for discs with gaps should work
well for intermediate-mass planets.

In this work we use a fixed planet orbital radius $R_{\mathrm{p}}$ and planet
mass $M_{\mathrm{p}}$ (resulting in time independent $S_{\zeta,\mathrm{sh}}$),
but we note that our 1D method can easily include planetary growth and a
prescribed planet migration by making these parameters (and
$S_{\zeta,\mathrm{sh}}$) functions of time (see also Section 7.1), in which
case the second term in (21) will no longer be linear in $t$.

####  4.2.2 Verification of $\Sigma(R),\Omega(R)$ retrieval via simulations

Figure 3: Radial profiles of (a) vortensity $\langle\zeta\rangle_{\phi}$ (in
units of $\Omega_{\mathrm{p}}/\Sigma_{\mathrm{p}}$), (b) surface density
(normalized by $\Sigma_{\mathrm{p}}$) and (c) its relative deviation from
$\Sigma_{\mathrm{i}}$, (d) angular frequency and (e) a measure of shear — the
second term in equation (25). Shown is the calculation for the fiducial disc
model with $M_{\mathrm{p}}=0.25M_{\mathrm{th}}$ at $t=700P_{\mathrm{p}}$. We
show results from the full hydro simulation (blue) and compare them to results
of application of our $\Sigma$, $\Omega$ retrieval procedure (Section 4.2.1)
with the vortensity profile from simulations (orange dashed) and from our
semi-analytical prescription (red dotted). The good agreement between the full
simulation results and our semi-analytical procedure is obvious.

We now check the performance of the $\Sigma(R),\Omega(R)$ retrieval procedure
described in Section 4.2 using simulations, by comparing the gap structures
obtained by both methods. In Fig. 3 we show the results for
$\Sigma(R),\Omega(R)$ at $t=700P_{\mathrm{p}}$ derived from simulations and
through our semi-analytical method for our fiducial disc model.

Panel (a) shows that the semi-analytical prediction $\zeta^{\mathrm{SA}}$ (red
dotted) provides an excellent match to the vortensity $\zeta^{\mathrm{sim}}$
(blue solid) obtained with Athena++, reliably reproducing the amplitude and
shape of the two vortensity rings. Slight disagreement is present in the inner
disc ($R\lesssim 0.7R_{\mathrm{p}}$), where $\zeta^{\mathrm{SA}}$ smoothly
tends to the unperturbed value, but the simulation shows a secondary bump and
trough. This is caused by the emergence of a secondary spiral arm in the
simulation that eventually shocks (Bae & Zhu, 2018; Miranda & Rafikov, 2019a).
Our semi-analytical model cannot capture this effect by design (see also Paper
I, ).

The azimuthally-averaged $\Sigma(R)$ and its relative deviation from the
initial $\Sigma_{\mathrm{i}}(R)$ profile, obtained by the two methods, are
shown in panels (b) and (c), respectively. As a test of our disc
reconstruction method, we also show the results obtained through the procedure
in Section 4.2 but with $\zeta^{\mathrm{sim}}(R,t)$ obtained from simulations
(orange dashed), rather than from equation (21). These show excellent
agreement with azimuthally averaged simulation results, proving the validity
of the $\Sigma(R),\Omega(R)$ retrieval technique (see also Lin & Papaloizou
2010). Note the characteristic double-trough profile of the gap, predicted in
Rafikov (2002b), with the two minima of $\Sigma$ corresponding to the density
wave shocking on each side of the planetary orbit.

When using $\zeta^{\mathrm{SA}}$ instead of $\zeta^{\mathrm{sim}}$ to solve
equation (23), the secondary surface density depression and enhancement in the
inner disc are missed by $\Sigma^{\mathrm{SA}}$, corresponding to the missing
bump and trough in $\zeta^{\mathrm{SA}}$. This leads to a smaller gradient at
the inner gap edge for the semi-analytical gap. Also, in the entire domain,
the semi-analytical solution for $\Sigma(R)$ shows a positive offset with
respect to simulation results. Nevertheless, the amplitude and width of the
primary bumps are matched reasonably well. In the outer disc, the agreement of
the surface density depression and enhancement is better due to the absence of
a secondary shock.

Despite these discrepancies in $\Sigma(R)$, the semi-analytical and numerical
radial profiles of the orbital frequency $\Omega(R)$ agree remarkably well,
see panel (d). For this parameter set, the outer gap edge shows stronger shear
(gradient of $\Omega$) and, quite importantly for the validity of our semi-
analytical analysis, $\Omega$ behaviour in the regions of strongest shear at
the gap edges is reproduced very well by our method.

According to the definition of vortensity, its deviation from the initial
value $\delta\zeta=\zeta-\zeta_{\mathrm{i}}$ depends on perturbations of both
$\Sigma$ and $\Omega$, i.e. on $\delta\Sigma=\Sigma-\Sigma_{\mathrm{i}}$ and
$\delta\Omega=\Omega-\Omega_{\mathrm{i}}$. One can easily show that as long as
$\Omega_{\mathrm{i}}\approx\Omega_{\mathrm{K}}$, we can express

$\displaystyle\frac{\delta\zeta}{\zeta_{\mathrm{i}}}$
$\displaystyle\approx\left(1+\frac{\delta\Sigma}{\Sigma_{\mathrm{i}}}\right)^{-1}\left[1+\frac{2}{\Omega_{\mathrm{K}}R}\frac{\mathrm{d}}{\mathrm{d}R}\left(R^{2}\delta\Omega\right)\right]-1$
(24)
$\displaystyle\approx-\frac{\delta\Sigma}{\Sigma_{\mathrm{i}}}+\frac{2}{\Omega_{\mathrm{K}}R}\frac{\mathrm{d}}{\mathrm{d}R}\left(R^{2}\delta\Omega\right),$
(25)

with the second line valid when $|\delta\Sigma|\ll\Sigma_{\mathrm{i}}$ and
$|\partial\delta\Omega/\partial R|\ll\Omega_{\mathrm{K}}/R$; equation (24)
holds even when these constraints are not met. The behaviour of the second
term in the right-hand side of (25) is illustrated in panel (e), while the
first one is shown in panel (c). One can see that even though
$|\delta\Omega|/\Omega_{\mathrm{K}}\lesssim 10^{-2}$ is rather small, see
panel (d), the term in panel (e) has a much higher amplitude $\sim 1$;
approximation (25) leads to substantial deviations near the peaks of $\zeta$
in this case. This term considerably exceeds
$\delta\Sigma/\Sigma_{\mathrm{i}}$ in magnitude, which implies that a proper
calculation of $\zeta$ evolution must account for the changes in both $\Sigma$
and $\Omega$.

To summarize, our $\Sigma$, $\Omega$ retrieval procedure shows good
performance when tested against simulations. In Appendix C we provide further
tests of the robustness of our semi-analytical reconstruction for various disc
parameters ($p$ and $h_{\mathrm{p}}$) and planetary masses, again finding good
agreement with simulations.

###  4.3 Illustration of our RWI analysis

We now show typical results obtained using linear stability analysis on gap
profiles produced both by simulations and by the semi-analytical method of
Section 4.2 for the fiducial disc model with
$M_{\mathrm{p}}=0.25M_{\mathrm{th}}$ and surface density slope $p=3/2$. For
this particular illustration, we use $\Sigma(R),\Omega(R)$ profiles derived
from the corresponding Athena++ simulation at $t=720P_{\mathrm{p}}$ (close in
time to Fig. 3). For comparison, in Appendix D we provide another illustration
of the RWI analysis, this time for a $p=0$ disc.

####  4.3.1 Typical mode structure

Figure 4: Top: Radial profile of the (normalized) unstable enthalpy eigenmode
for $m=3$, associated with the outer gap edge, for the fiducial disc
simulation at $t=720P_{\mathrm{p}}$. We show the real part, imaginary part and
absolute magnitude of $\Psi_{3}$. The frequency of this mode is
$\omega=2.3374+0.01653\mathrm{i}$. Bottom: Two-dimensional map of the
corresponding pressure perturbation.

We start by presenting the typical radial structure of RWI unstable modes,
obtained from linear stability analysis of the fiducial Athena++ model. In
Fig. 4 we show the radial profile of the normalized enthalpy perturbation
$\Psi_{m}/\max\left|{\Psi_{m}}\right|$ for $m=3$, since this $m$ dominates the
RWI in the fiducial run, see Fig. 1. We plot the (normalized) real part,
imaginary part and the absolute magnitude of $\Psi_{3}$ in the top panel. The
phase is arbitrarily set such that $\mathrm{Im}(\Psi_{m})=0$ at the location
of $\max|\Psi_{m}|$666This corresponds to an azimuthal rotation, under which
the axisymmetric problem is invariant.. One can see that $\Psi_{3}$ strongly
peaks at the vortensity minimum (as expected for trapped modes based on
previous studies, e.g. Papaloizou & Lin 1989, Li et al. 2000, Lin & Papaloizou
2010) at the outer gap edge, which becomes unstable to RWI earlier than the
inner edge. This is expected for $p=3/2$ disc, see Paper I for details. We
also find that spacial trapping of the modes depends on their azimuthal
wavenumber $m$, with higher-$m$ modes being less radially confined near the
vortensity minimum.

The $m=3$ mode is close to corotating with the background flow
($\omega_{R}\simeq m\Omega(R_{0})$) at the vortensity minimum at $R_{0}\simeq
1.17R_{\mathrm{p}}$ (see Fig. 3a) and has a growth rate corresponding to an
e-folding time of about 10 orbits, indicating growth on dynamical timescales.
To give an example of the 2D mode structure, we show the $\phi-R$ map of the
corresponding normalized pressure perturbation $\delta P_{m}/\max\left|{\delta
P_{m}}\right|$ in the bottom panel of Fig. 4. The pressure perturbations
($\delta P_{m}=\Sigma_{0}\Psi_{m}$) show $m=3$ regions of higher (red) and
lower (blue) pressure at the outer gap edge, characteristic of vortices.

####  4.3.2 Time-dependence of growth rates and most unstable modes

Figure 5: Radial profile of linearly unstable enthalpy eigenmode $\Psi_{3}(R)$
for $m=3$, associated with the outer gap edge, for the fiducial disc model at
$t=610P_{\mathrm{p}}$ (blue solid) and $740P_{\mathrm{p}}$ (orange dashed).
One can see that until the instability develops into the non-linear phase, the
localized part of the unstable mode does not change significantly, consistent
with the mode growing over long timescales while $\gamma$ changes.  Figure 6:
Instantaneous growth rates $\gamma_{m}$ of linearly unstable modes, associated
with the outer gap edge, found for different $m$ at different times for the
fiducial disc and $M_{\mathrm{p}}/M_{\mathrm{th}}=0.25$. We compare results of
linear RWI analysis (Section 4.1) using disc profiles obtained directly from
simulations (blue, downward tripods) and reconstructed using our semi-
analytical framework in Section 4.2 (orange, upward tripod). The vertical red
dotted line marks $t=750$ orbits, when growth rates begin falling off. At this
time we observe clearly $m=3$ vortices in the simulation, marking the
transition into the non-linear regime of the instability.

As the disc evolves under the action of the planetary perturbation, so do the
inputs $\Sigma(R),\Omega(R)$ for the modal equation (9). As a result, the RWI
eigenmodes and, most importantly, eigenfrequencies change in time.

Evolution of the eigenmodes is illustrated in Fig. 5, where we show the radial
structure of $m=3$ mode at $t=610P_{\mathrm{p}}$ (blue, solid) and
$t=740P_{\mathrm{p}}$ (orange, dashed) in our fiducial setup with
$M_{\mathrm{p}}=0.25M_{\mathrm{th}}$. One can see that the mode structure at
different times is rather similar, with the largest differences in
$\mathrm{Im}(\Psi_{3})$ and considerably smaller changes in
$\mathrm{Re}(\Psi_{3})$. We can also see that the radial behaviour of the mode
phase changes only slightly. These results imply that the radial structure of
the unstable modes change rather slowly during the time in which the
instability develops, a property that will later be used in Section 5. We find
similar results for all other $m$ we consider.

On the other hand, the growth rate of the mode changes significantly during
this interval of time. Indeed, the corresponding eigenfrequencies are
$\omega_{3}(t=610P_{\mathrm{p}})=2.352+0.00566\mathrm{i}$ and
$\omega_{3}(t=740P_{\mathrm{p}})=2.335+0.0184\mathrm{i}$, indicating
exponential growth over roughly 28 and 9 planet orbits at these moments of
time. As $|\mathrm{Re}(\omega)|\gg|\mathrm{Im}(\omega)|=\gamma$, the relative
changes of the growth rates are much larger than those of
$|\mathrm{Re}(\omega)|$. This is typical for all linearly unstable RWI modes.

Figure 6 provides a more detailed illustration of the rapid evolution of the
growth rates obtained via our RWI linear stability analysis, and does this for
a number of RWI modes ($m=1$ to $m=6$ in different panels). Orange upward
tripods show the growth rates obtained by using $\Sigma(R)$, $\Omega(R)$ given
by our semi-analytical reconstruction (Section 4.2) in equation (9). Blue
downward tripods show growth rates obtained when $\Sigma(R)$, $\Omega(R)$ are
taken directly from our fiducial simulation. Growth rates are computed every
10 orbits; note that the horizontal axis does not start at $t=0$, but at the
moment when unstable modes first appear (for $m=1$).

For both methods of computing the growth rates, unstable modes are found
earliest for the lowest azimuthal wave number $m=1$ (ground state). For all
$m$, growth rates increase close to linearly with time up until the non-linear
phase of instability begins. This is not surprising since the gap opening
process continues and the disc becomes more and more unstable (de Val-Borro et
al., 2007). The slope of $\gamma_{m}(t)$ increases with $m$, most notably from
$m=1$ to $m=3$. This means that the instantaneously most unstable mode will
change in time, with higher-$m$ modes with higher $\gamma$ dominating at later
time. For all $m$, we find an offset in time between $\gamma_{m}(t)$ obtained
by the two methods $t_{\mathrm{off}}\simeq 30P_{\mathrm{p}}$, with the semi-
analytical disc becoming unstable later. This offset does not depend strongly
on $m$ and the slope of $\gamma_{m}(t)$ is comparable between the two methods.
These results indicate that the semi-analytical method captures the essential
physics and is well suited to determine the timescale for development of the
RWI.

Note that the blue tripods, corresponding to $\Sigma$, $\Omega$ derived from a
simulation, begin falling off at around $t=750P_{\mathrm{p}}$, which is when
strong vortices have developed and the instability transitions into the non-
linear regime. This is because perturbations are then strong enough to make
the azimuthally averaged disc structure more stable against the RWI777We note
that linear stability analysis is not a suitable method at this stage due to
the large amplitude of perturbations.. On the other hand, the growth rates
obtained using the semi-analytical reconstruction do not show such a break and
keep increasing monotonically. This is to be expected since there is no
instability feedback built into our simple model that would stop $\delta\zeta$
from steadily growing near the planet.

##  5 Definition: characteristic timescales of the RWI

input: $(M_{\mathrm{p}}/M_{\mathrm{th}},h_{\mathrm{p}},p)$vortensity
$\zeta_{0}^{\mathrm{SA}}(R,t)$ from semi-analytical (SA) framework (Paper
I)non-linear simulations (Athena++)
(§3)$\left[\Sigma(R,\phi,t),\Omega(R,\phi,t)\right]^{\mathrm{sim}}$$\left[\Sigma_{0}(R,t),\Omega_{0}(R,t)\right]^{\mathrm{SA}}$linear
stability analysis (§4)azimuthal average
$\langle\cdot\rangle_{\phi}$$\left[\Sigma_{0}(R,t),\Omega_{0}(R,t)\right]^{\mathrm{sim}}$direct
Fourier analysis (§5)growth
rates$\gamma_{m}^{\mathrm{SA}}(t)$$\gamma_{m}^{\mathrm{sim}}(t)$Integrate
$\int\gamma_{m}(t^{\prime})\,\mathrm{d}t^{\prime}$ $\rightarrow$ amplification
factor $\mathcal{A}_{m}(t)$
(§5)$t^{\mathrm{SA}}_{\mathrm{nl}}$$t^{\mathrm{SA}}_{\mathrm{lin}}$$t^{\mathrm{sim}}_{\mathrm{lin}}$first
unstable
modes$t^{\mathrm{sim}}_{\mathrm{nl}}$$\hat{t}^{\mathrm{sim}}_{\mathrm{nl}}$developed
vortices solve Eq. (23) break in $\sigma_{m}(t)$ (Fig.
8)$\gamma_{m}>\gamma_{\mathrm{c}}$$\mathcal{A}_{m}>{\mathcal{A}_{m}}_{,\mathrm{c}}$
Figure 7: Flow chart showing the methods used in our study (rectangles) and
their outputs (rounded rectangles and circles). Teal fill marks purely semi-
analytical methods, orange fill marks full simulations and no fill marks
general methods. The timescales as final outputs are marked by circles. The
times when linear stability analysis first predicts unstable modes are
$t^{\mathrm{SA}}_{\mathrm{lin}}$ and $t^{\mathrm{sim}}_{\mathrm{lin}}$, marked
by the dashed circles on the left and grouped by the dashed rounded rectangle
on the left. The three timescales characterizing substantial RWI growth —
$t_{\mathrm{nl}}^{\mathrm{SA}}$, $t^{\mathrm{sim}}_{\mathrm{nl}}$,
$\hat{t}^{\mathrm{sim}}_{\mathrm{nl}}$ — are grouped in the dashed rounded
rectangle on the right. See Section 6 for details.

Determination of the timescale for the emergence of vortices at the edges of
planet-induced gaps is the main goal of this paper. To that end, we now define
a number of important timescales that characterize different stages of the
vortex development via the RWI. Given the variety of methods that we employ to
study this problem, some of the timescales also depend on the particular
method used in defining them, see Fig. 7 for illustration. In particular, we
use the superscript ’SA’ for the timescales based on the semi-analytical disc
reconstruction technique of Section 4.2, while the ones utilizing the
simulation data in any way are denoted ’sim’. We also use different timescales
to characterize the onset and non-linear development of the instability, as
described below.

###  5.1 Timescales based on linear RWI stability analysis

As the planet opens a deeper and deeper gap, the gap edge will at some point
become unstable, which results in RWI modes computed by our linear analysis
attaining $\gamma_{m}>0$. This allows us to define several characteristic
timescales, as follows.

$t^{\mathrm{SA}}_{\mathrm{lin}}$: onset of instability with semi-analytic
reconstruction
We define $t^{\mathrm{SA}}_{\mathrm{lin}}$ as the time after which linear
stability analysis first predicts unstable RWI modes with a growth rate
greater than a critical value $\gamma_{m}>\gamma_{\mathrm{crit}}=10^{-2}/2\pi$
(corresponding to exponential growth with an e-folding time of 100 planet
orbits). This particular definition assumes that the radial structure of the
disc used in our linear stability analysis has been obtained through the semi-
analytical reconstruction described in Section 4.2.

$t^{\mathrm{sim}}_{\mathrm{lin}}$: onset of instability with simulation input
The time $t^{\mathrm{sim}}_{\mathrm{lin}}$ is defined analogous to
$t^{\mathrm{SA}}_{\mathrm{lin}}$ but using azimuthal averages of 2D simulation
results as $\Sigma(R),\Omega(R)$ inputs for the RWI linear stability analysis.

Figure 8: Logarithmic measure of the relative amplitudes of the Fourier
components of the surface density perturbation $|\Sigma_{m}|/\Sigma$ measured
in our runs, integrated over a radial range around the outer gap edge
($\sigma_{m}$, see equation 30) and shown as a function of time for different
$m$. The break in the behaviour of $\sigma_{m}$ indicates the transition to
the non-linear phase of RWI in our simulations.

$t_{\mathrm{nl}}^{\mathrm{SA}}$: developed RWI with semi-analytic
reconstruction
While the above procedure predicts when linear instability first sets in, it
does not tell us when the instability has amplified to a point where
substantial vortices appear.

Results of Section 4.3.2 indicate that the radial structure of an unstable
mode does not vary much as the instability develops, while the growth rate
changes substantially. If we take two moments of time during this phase,
$t_{\mathrm{i}}$ and $t_{\mathrm{f}}$, then the mode amplitudes at these times
are related by the integral of instantaneous growth-rates over time:

$\displaystyle\ln\frac{\left|{\delta\Sigma_{m}(R,t_{\mathrm{f}})}\right|}{\left|{\delta\Sigma_{m}(R,t_{\mathrm{i}})}\right|}=\int_{t_{\mathrm{i}}}^{t_{\mathrm{f}}}\gamma_{m}\,\mathrm{d}t\equiv\Gamma_{m}.$
(26)

In other words,

$\displaystyle\left|{\delta\Sigma_{m}(R,t_{\mathrm{f}})}\right|=\exp\left(\Gamma_{m}\right)\left|{\delta\Sigma_{m}(R,t_{\mathrm{i}})}\right|=\mathcal{A}_{m}\left|{\delta\Sigma_{m}(R,t_{\mathrm{i}})}\right|,$
(27)

where we have defined the amplification factor
$\mathcal{A}_{m}\equiv\exp\left(\Gamma_{m}\right)$.

Using this logic, we can define a timescale $t_{\mathrm{nl}}^{\mathrm{SA}}$ as
the time over which an initial perturbation of order
$|\delta\Sigma|/\Sigma\simeq\mathcal{A}_{m,\mathrm{crit}}^{-1}$ grows to order
unity:

$\displaystyle\mathcal{A}_{m}(t_{\mathrm{nl}})\geq\mathcal{A}_{m,\mathrm{crit}},$
(28)

where $\mathcal{A}_{m,\mathrm{crit}}$ is some critical value (see next) and
the growth rates in (26) are computed using the inputs based on semi-
analytical reconstruction. By definition, we always have
$t^{\mathrm{SA}}_{\mathrm{lin}}<t^{\mathrm{SA}}_{\mathrm{nl}}$.

By comparing $t^{\mathrm{SA}}_{\mathrm{nl}}$ with the time when non-linear
simulations show substantial RWI vortices (see Section 5.2), we find good
agreement when setting $\mathcal{A}_{m,\mathrm{crit}}=10^{4}$, which gives an
estimate of initial deviations from axisymmetry ($\delta\Sigma/\Sigma\simeq
10^{-4}$). Due to the exponential growth of the instability, the timescale is
not very sensitive to this choice. Here we assume that the initial noise is of
similar amplitude for all mode numbers. In real systems, there might be
preferential forcing of specific modes. This might lead to a different offset
between $t^{\mathrm{SA}}_{\mathrm{lin}}$ and $t^{\mathrm{SA}}_{\mathrm{nl}}$
as $m$ varies.

$t^{\mathrm{sim}}_{\mathrm{nl}}$: developed RWI with simulation inputs
We also define another time $t^{\mathrm{sim}}_{\mathrm{nl}}$ for developed
RWI, analogous to $t^{\mathrm{SA}}_{\mathrm{nl}}$ but now using the
azimuthally-averaged simulation data when computing the growth rates
$\gamma_{m}$, which enter the equation (26).

###  5.2 $\hat{t}^{\mathrm{sim}}_{\mathrm{nl}}$: simulation-only timescale

Finally, we also diagnose the non-linear development of the instability using
a time $\hat{t}^{\mathrm{sim}}_{\mathrm{nl}}$ defined using only the
simulation results, i.e. not resorting to the RWI linear stability analysis as
we do for all other characteristic times defined in Section 5.1. This is done
by tracking the non-axisymmetric density structure in the simulations, the
amplitude of which is interpreted as indicating the RWI development.

In practice, we first perform Fourier decomposition of the surface density
perturbation of our 2D simulation data:

$\displaystyle\Sigma_{m}(R,t)\equiv\frac{1}{2\pi}\int_{0}^{2\pi}\Sigma(R,\phi,t)\exp(-\mathrm{i}m\phi)\,\mathrm{d}\phi,$
(29)

with the prefactor chosen such that
$\left\langle\Sigma\right\rangle_{\phi}=\Sigma_{0}$. We then define the
quantity

$\displaystyle\sigma_{m}\equiv\ln\left(\int_{R_{\mathrm{l}}}^{R_{\mathrm{u}}}\frac{\left|{\Sigma_{m}}\right|}{\Sigma_{0}}\,R\mathrm{d}R\right),$
(30)

where we perform the integral over a radial region
$R_{\mathrm{l}}<R<R_{\mathrm{u}}$ around the gap edge, where RWI unstable
modes peak in amplitude; at the outer gap edge we typically take
$R_{\mathrm{l}}=1.05R_{\mathrm{p}}$ and $R_{\mathrm{u}}=1.5R_{\mathrm{p}}$. A
similar metric has been used for example by Les & Lin (2015).

Once the instability develops, we expect $\sigma_{m}$ to become non-zero and
to grow linearly in time allowing us to diagnose the onset of RWI. However,
there is a complication: due to the non-axisymmetric perturbations induced by
the planet, $\left|{\Sigma_{m}}\right|\neq 0$ at all times, even before the
gap edge becomes unstable. This is illustrated in Fig. 8, where we plot
$\sigma_{m}(t)$ for the fiducial simulation for $1\leq m\leq 8$. One can see
that, indeed, $\sigma_{m}$ is always non-zero and shows slow variation early
on. The weak time dependence is caused by the effect of gap opening on the
planet-induced density waves888The radial locations of Lindblad resonances
depend on their order, such that local changes to the relative surface density
will affect the excitation of modes with different $m$ differently (e.g.
Petrovich & Rafikov, 2012)..

Nevertheless, we can still associate the onset of instability with the break
in the behaviour of $\sigma_{m}(t)$, which is clearly present in all panels.
Indeed, beginning at around $t=720P_{\mathrm{p}}$, small oscillations occur
for the three ($1\leq m\leq 3$) lowest mode numbers, which are most pronounced
for $m=2$ and $m=3$. Around $t=760P_{\mathrm{p}}$, there is a clear break in
$\sigma_{m}(t)$ behaviour of the $m=3$ mode (red dotted line vertical in panel
c), shortly after which all $1\leq m\leq 3$ show strong growth999Modes with
$m\geq 4$ show more noisy behaviour of $\sigma_{m}(t)$ instead of coherent
growth. of $\sigma_{m}(t)$, indicating (close to exponential) growth of
$\Sigma_{m}$ for around 50 orbits, after which $\sigma_{m}$ saturates. The
break thus corresponds to a time, which we call
$\hat{t}^{\mathrm{sim}}_{\mathrm{nl}}$, when RWI has grown substantially and
started dominating101010Note that $\hat{t}^{\mathrm{sim}}_{\mathrm{nl}}$
occurs slightly earlier than RWI saturation in simulations (which is not very
well defined). However the relative difference between the two moments of time
is typically $\lesssim 5\%$, so we define
$\hat{t}^{\mathrm{sim}}_{\mathrm{nl}}$ based on the initial break of
$\sigma_{m}(t)$. over the non-axisymmetric perturbation due to the planet.
This interpretation is confirmed by examining Figs. 1 & 2, which show that at
$t\approx\hat{t}^{\mathrm{sim}}_{\mathrm{nl}}=760P_{\mathrm{p}}$ vortices
become visible and comparable in amplitude to planetary perturbations in 2D
maps of $\Sigma$ and vortensity perturbations. As the determination of
$\hat{t}^{\mathrm{sim}}_{\mathrm{nl}}$ uses only the simulation data, it
allows for an independent verification of our semi-analytical calculations.

To summarize, we can define five different characteristic times, two of which
— $t^{\mathrm{SA}}_{\mathrm{lin}}$ and $t^{\mathrm{sim}}_{\mathrm{lin}}$ —
describe first appearance of unstable modes, while three others —
$t^{\mathrm{SA}}_{\mathrm{nl}}$, $t^{\mathrm{sim}}_{\mathrm{nl}}$, and
$\hat{t}^{\mathrm{sim}}_{\mathrm{nl}}$ — mark the development of the non-
linear phase of the RWI. These different metrics are compared in the following
section.

##  6 Results: timescales for instability

Table 1: Complete list of hydro simulations that were performed. We report planet mass $M_{\mathrm{p}}$ (in units of $M_{\mathrm{th}}$), disc parameters ($p$, $h_{\mathrm{p}}$), simulation time (in orbits) at which we observe the RWI at inner and outer gap edges, respectively. Last column indicates which gap edge shows vortices first. Our fiducial simulation is shown in boldface.  $M_{\mathrm{p}}/M_{\mathrm{th}}$ | $p$ | $h_{\mathrm{p}}$ | $\hat{t}^{\mathrm{sim}}_{\mathrm{inner}}/P_{\mathrm{p}}$ | $\hat{t}^{\mathrm{sim}}_{\mathrm{outer}}/P_{\mathrm{p}}$ | most unstable gap edge
---|---|---|---|---|---
0.1 | 1.5 | 0.05 | $>6800$ | $>6800$ | -
0.1 | 1.5 | 0.07 | $>4200$ | $>4200$ | -
0.1 | 1.5 | 0.10 | $>2970$ | $>2970$ | -
0.1 | 0 | 0.05 | $>5850$ | $>5850$ | -
0.1 | 0 | 0.07 | $>4750$ | $>4750$ | -
0.1 | 0 | 0.10 | $>2800$ | $>2800$ | -
0.25 | 1.5 | 0.05 | 880 | 760 | outer
0.25 | 1.5 | 0.07 | 690 | 540 | outer
0.25 | 1.5 | 0.10 | 520 | 380 | outer
0.25 | 0 | 0.05 | 710 | 1000 | inner
0.25 | 0 | 0.07 | 510 | 810 | inner
0.25 | 0 | 0.10 | 330 | 450 | inner
0.5 | 1.5 | 0.05 | 164 | 146 | outer
0.5 | 1.5 | 0.07 | 130 | 112 | outer
0.5 | 1.5 | 0.10 | 94 | 82 | outer
0.5 | 0 | 0.05 | 144 | 168 | inner
0.5 | 0 | 0.07 | 98 | 128 | inner
0.5 | 0 | 0.10 | 72 | 112 | inner

We now describe the results on the various instability timescales (defined in
previous section) obtained by applying the methods described in Sections 3 & 4
to a variety of disc-planet setups. We vary a number of parameters of the
problem, choosing three values each for the normalized planet mass111111The
lower limit on $(M_{\mathrm{p}}/M_{\mathrm{th}})$ is due to the requirement of
increasingly high resolutions to correctly capture vortensity generation by
low mass planets (see Paper I), coupled with increasing timescale for the
onset of RWI. $M_{\mathrm{p}}/M_{\mathrm{th}}\in\\{0.1,0.25,0.5\\}$, and the
disc scale-height $h_{\mathrm{p}}\in\\{0.05,0.07,0.1\\}$ and two values for
the surface density slope of the background disc $p\in\\{0,1.5\\}$. Covering
all possible combinations, this gives a total of 18 disc models listed in
Table 1.

For every disc model we run a direct simulation to get the simulation-based
timescales $t^{\mathrm{sim}}_{\mathrm{lin}}$,
$t^{\mathrm{sim}}_{\mathrm{nl}}$, and $\hat{t}^{\mathrm{sim}}_{\mathrm{nl}}$,
and use the semi-analytical reconstruction of Section 4.2 to obtain
$t^{\mathrm{SA}}_{\mathrm{lin}}$, $t^{\mathrm{SA}}_{\mathrm{nl}}$. For
$M_{\mathrm{p}}=0.1M_{\mathrm{th}}$, simulations could not be run long enough
to develop RWI, providing only a lower limit on the instability timescale
$\hat{t}^{\mathrm{sim}}_{\mathrm{nl}}$. To carry out our linear analysis
(Section 5.1) and to observe the development of RWI directly (Section 5.2), we
sample the disc state every 10 orbits in runs with
$M_{\mathrm{p}}=0.25M_{\mathrm{th}}$, and every 2 orbits in the highest
$M_{\mathrm{p}}=0.5M_{\mathrm{th}}$ runs.

###  6.1 First appearance of unstable modes

Figure 9:  Characteristic timescales for the appearance of first unstable RWI
modes. We show $t^{\mathrm{SA}}_{\mathrm{lin}}$ (crosses) for all parameter
sets and $t^{\mathrm{sim}}_{\mathrm{lin}}$ (circles) where available, see
Section 5.1 for definitions. We show results for (a) the fiducial disc
($p=1.5$) and (b) the constant surface density disc $(p=0)$. Colours
correspond to different values of $h_{\mathrm{p}}$. The dashed lines show a
power law fit for $t^{\mathrm{SA}}_{\mathrm{lin}}$ given by equation (31) with
parameters (33).

We start by presenting results for $t^{\mathrm{SA}}_{\mathrm{lin}}$ and
$t^{\mathrm{sim}}_{\mathrm{lin}}$ — the times when linearly unstable modes
with $\gamma_{m}>\gamma_{\mathrm{crit}}$ first appear (see Section 5.1) —
allowing us to compare RWI linear analysis with inputs based on simulations
and semi-analytical method. In Fig. 9, we plot
$t^{\mathrm{SA}}_{\mathrm{lin}}$ as crosses for different $h_{\mathrm{p}}$
(colours) as a function of normalized planet mass on a log-log scale. In
panels (a) and (b) we consider discs with $p=1.5$ and $p=0$, examining onset
of RWI at the outer and inner gap edge, respectively.

For all parameter sets except those including
$M_{\mathrm{p}}=0.1M_{\mathrm{th}}$, we also show
$t^{\mathrm{sim}}_{\mathrm{lin}}$ (circles). For the lowest planet mass
explored, $M_{\mathrm{p}}=0.1M_{\mathrm{th}}$, simulations do not show RWI
vortices at the end of run time indicated in Table 1, e.g. after
$t_{\mathrm{end}}=6800P_{\mathrm{p}}$ for the disc with fiducial
parameters121212This run required more than half a year of continuous
computation.. This is in full agreement with the estimates of
$t^{\mathrm{SA}}_{\mathrm{lin}}$ for these runs, which all exceed the run
times, see Fig. 9.

We can make several observations based on this figure. First, there is
generally good agreement between $t^{\mathrm{SA}}_{\mathrm{lin}}$ and
$t^{\mathrm{sim}}_{\mathrm{lin}}$ for any given disc-planet model, confirming
the robustness of our semi-analytic reconstruction technique. The agreement is
slightly worse for models with $p=0$, but this is to be expected because of
the secondary arm formation in the inner disc, see Appendix D.

Second, for both values of $p$, the behaviour of
$t^{\mathrm{SA}}_{\mathrm{lin}}$ and $t^{\mathrm{sim}}_{\mathrm{lin}}$ in Fig.
9 suggests a power law fit for the time $t_{\mathrm{lin}}$ when the first
unstable modes appear of the form (we will use this fit also in the next
section for $t_{\mathrm{nl}}$)

$\displaystyle\\{t_{\mathrm{lin}},t_{\mathrm{nl}}\\}=AP_{\mathrm{p}}\left(\frac{M_{\mathrm{p}}}{M_{\mathrm{th}}}\right)^{\alpha}h_{\mathrm{p}}^{\beta}$
(31)

or, equivalently,

$\displaystyle\log\frac{\\{t_{\mathrm{lin}},t_{\mathrm{nl}}\\}}{P_{\mathrm{p}}}=\log
A+\alpha\log\frac{M_{\mathrm{p}}}{M_{\mathrm{th}}}+\beta\log h_{\mathrm{p}}.$
(32)

To obtain the parameters of this fit, we perform a least-squares regression of
$t^{\mathrm{SA}}_{\mathrm{lin}}$ data (as we have them available for all
$M_{\mathrm{p}}$, even the lowest-mass case
$M_{\mathrm{p}}=0.1M_{\mathrm{th}}$) in log-log space over all 18 data points,
including $p=1.5$ and $p=0$. Our resultant fit of the form (31) is shown via
dashed lines in Fig. 9 and has parameters

$\displaystyle(\log A,\alpha,\beta)_{t_{\mathrm{lin}}}=(-0.63\pm 0.09,-3.23\pm
0.05,-1.13\pm 0.08).$ (33)

so that $A\approx 0.23$. The maximum relative deviation between this fit and
the $t^{\mathrm{SA}}_{\mathrm{lin}}$ data is 11% over all points, and we see a
good match for both surface density slopes.

The power law indices indicate that the planet mass and disc scale-height have
the strongest influence on $t_{\mathrm{lin}}$, with the former dominating, and
the surface density slope having insignificant effect.

###  6.2 Timescale for developed RWI (vortex formation)

Figure 10:  Timescales characterizing the non-linear development of RWI:
$t^{\mathrm{SA}}_{\mathrm{nl}}$ (pluses), $t^{\mathrm{sim}}_{\mathrm{nl}}$
(circles), and $\hat{t}^{\mathrm{sim}}_{\mathrm{nl}}$ (crosses), see Sections
5.1,5.2 for definitions. Panels again show (a) fiducial disc ($p=1.5$) and (b)
constant surface density disc $(p=0)$, colours indicate different
$h_{\mathrm{p}}$. The dashed lines show a power law fit (31) with parameters
(34) for all 18 data points for $t^{\mathrm{SA}}_{\mathrm{nl}}$.

Perhaps more important than $t_{\mathrm{lin}}$ from a practical standpoint is
the time-scale $t_{\mathrm{nl}}$ for well-developed RWI, indicating when
large-amplitude vortices appear and can be observed directly. We use the three
previously defined (Section 5) time scales — $t^{\mathrm{SA}}_{\mathrm{nl}}$,
$t^{\mathrm{sim}}_{\mathrm{nl}}$, and $\hat{t}^{\mathrm{sim}}_{\mathrm{nl}}$ —
as proxies for $t_{\mathrm{nl}}$ for all parameter sets, except those with
$M_{\mathrm{p}}/M_{\mathrm{th}}=0.1$ for which only
$t^{\mathrm{SA}}_{\mathrm{nl}}$ is available, since simulations have not
become unstable by the end of the run time.

In Fig. 10, we show these timescales in a way similar to Fig. 9, see the
legend. The agreement we find between $t^{\mathrm{sim}}_{\mathrm{nl}}$ (RWI
analysis using simulation data as input) and
$\hat{t}^{\mathrm{sim}}_{\mathrm{nl}}$ (timescale based directly on
simulations) is excellent for all simulations, with maximum relative
deviations between the timescales of 5%. This justifies our method of
integrating the linear growth rates (see Section 5.1) as a good way of
predicting the emergence of fully-developed vortices. We also find that the
semi-analytical estimate $t^{\mathrm{SA}}_{\mathrm{nl}}$ gives a maximum
relative deviation of 15% compared to $\hat{t}^{\mathrm{sim}}_{\mathrm{nl}}$,
once again, supporting our gap reconstruction technique of Section 4.2. We use
this relative deviation as a simple estimate for the uncertainty of
$t^{\mathrm{SA}}_{\mathrm{nl}}$ below.

Again, a least square fit for $t_{\mathrm{nl}}$ in the form (32) is
appropriate, and we obtain its parameters by linear regression of
$t^{\mathrm{SA}}_{\mathrm{nl}}$ data as

$\displaystyle(\log A,\alpha,\beta)_{t_{\mathrm{nl}}}=(0.21\pm 0.13,-2.68\pm
0.06,-0.86\pm 0.11).$ (34)

so that $A\approx 1.6$. With these parameters, equation (31) yields the vortex
emergence timescale in physical units as

$\displaystyle\tau$ $\displaystyle\approx 4\times
10^{3}\mbox{yr}\left(\frac{M_{\mathrm{p}}}{M_{\mathrm{th}}}\right)^{-2.7}\left(\frac{R_{\mathrm{p}}}{50\mbox{AU}}\right)^{1.5}\left(\frac{h_{\mathrm{p}}}{0.1}\right)^{-0.86}\left(\frac{M_{\star}}{M_{\odot}}\right)^{-0.5}$
(35) $\displaystyle\approx 3\times
10^{4}\mbox{yr}\left(\frac{M_{\mathrm{p}}}{0.5M_{\mathrm{J}}}\right)^{-2.7}\left(\frac{R_{\mathrm{p}}}{50\mbox{AU}}\right)^{1.5}\left(\frac{h_{\mathrm{p}}}{0.1}\right)^{7.2}\left(\frac{M_{\star}}{M_{\odot}}\right)^{2.2}.$
(36)

In agreement with the results of Paper I, we find that in a disc with constant
surface density ($p=0$), the inner gap edge is more unstable than the outer
one, opposite to the findings for the fiducial $p=3/2$ disc. We also see that
the RWI stability of the inner gap edge is affected to a certain degree by the
formation and shocking of a secondary (and higher order) spiral arms in the
inner disc (Bae & Zhu, 2018; Miranda & Rafikov, 2019a). This effect is not
captured by our gap reconstruction technique (Section 4.2), while being fully
accounted for in simulations. Nevertheless, our semi-analytical predictions
for the RWI development at the inner gap edge (e.g. in $p=0$ discs) match
$\hat{t}^{\mathrm{sim}}_{\mathrm{nl}}$ quite well, see Fig. 10.

###  6.3 Theoretical estimate

From our previous discussion (Section 4.1) and the results of Paper I it is
natural to expect that the timescale $\tau$ on which the disc structure
changes and the RWI sets in should depend on the vortensity production rate by
the planetary shocks $S_{\zeta,\mathrm{sh}}$. According to equation (20),
$S_{\zeta,\mathrm{sh}}$ is given by the product of the vortensity jump at the
shock $\Delta\zeta$ and the synodic frequency of a shocked gas parcel w.r.t.
the planet. To lowest order, we can estimate the synodic frequency using the
local approximation, i.e. expanding it to linear order in $R-R_{\mathrm{p}}$:

$\displaystyle|\Omega(R)-\Omega_{\mathrm{p}}|\simeq\frac{3\Omega_{\mathrm{p}}}{2}\left|\frac{R-R_{\mathrm{p}}}{R_{\mathrm{p}}}\right|.$
(37)

Figure 11 illustrates the performance of this approximation, with the local
expansion (37) shown in orange being compared to
$|\Omega-\Omega_{\mathrm{p}}|$ (solid blue). Vertical coloured lines
correspond to $R=R_{\mathrm{p}}\pm 3l_{\mathrm{sh}}$, computed for different
values of $h_{\mathrm{p}}$ and $M_{\mathrm{p}}/M_{\mathrm{th}}$: we know from
Paper I that $\Delta\zeta$ is significant only at radii satisfying
$l_{\mathrm{sh}}\lesssim|R-R_{\mathrm{p}}|\lesssim 3l_{\mathrm{sh}}$ (for
$M_{\mathrm{p}}\lesssim M_{\mathrm{th}}$). One can see that within this range
the local approximation (37) provides a decent fit for small
$l_{\mathrm{sh}}$, with order unity deviations becoming particularly
noticeable for highest $h_{\mathrm{p}}$ and lowest
$M_{\mathrm{p}}/M_{\mathrm{th}}$, mainly in the inner disc. This has important
implications (see Section 6.4), but for the sake of our simple argument we
will adopt the approximation (37) in what follows.

Figure 11: Relative orbital frequency of gas parcels w.r.t. the planet-driven
shock. We show the full expression (blue) and its linearised form (equation
(37), orange). Vertical lines indicate a radial separation
$|R-R_{\mathrm{p}}|=3l_{\mathrm{sh}}$ from the planet for different
$M_{\mathrm{p}}/M_{\mathrm{th}}$ (colours) and $h_{\mathrm{p}}=0.05$ and
$h_{\mathrm{p}}=0.1$ (dashed and dash-dotted lines, respectively).

Moreover, we have shown in Paper I that the vortensity jump at the shock
$\Delta\zeta\propto(M_{\mathrm{p}}/M_{\mathrm{th}})^{3}$, but is a weak
function of the disc aspect ratio $h_{\mathrm{p}}$. With this in mind, using
$|R-R_{\mathrm{p}}|\sim l_{\mathrm{sh}}$ and recalling that
$l_{\mathrm{sh}}\approx
R_{\mathrm{p}}h_{\mathrm{p}}(M_{\mathrm{p}}/M_{\mathrm{th}})^{-2/5}$ (Goodman
& Rafikov, 2001), equations (20) and (37) yield

$\displaystyle\frac{\partial\zeta}{\partial
t}=S_{\zeta,\mathrm{sh}}\propto\left(\frac{M_{\mathrm{p}}}{M_{\mathrm{th}}}\right)^{3}\times\Omega_{\mathrm{p}}\frac{l_{\mathrm{sh}}}{R_{\mathrm{p}}}\propto\Omega_{\mathrm{p}}h_{\mathrm{p}}\left(\frac{M_{\mathrm{p}}}{M_{\mathrm{th}}}\right)^{2.6}.$
(38)

Note that $S_{\zeta,\mathrm{sh}}$ depends on $h_{\mathrm{p}}$ both explicitly
and through $M_{\mathrm{th}}$, see equation (1).

Let us now make a simple assumption that RWI sets in and vortices appear at
the time $\tau$ corresponding to the moment when $\delta\zeta$ (e.g. its peak
or trough value, ignoring the radial structure for simplicity) reaches some
threshold value $\zeta_{\mathrm{th}}$; we discuss the validity of this
assumption next in Section 6.4. In other words, in this approach $\tau$ should
be determined from the condition

$\displaystyle\int^{\tau}_{0}S_{\zeta,\mathrm{sh}}(t)\mathrm{d}t\approx\zeta_{\mathrm{th}}.$
(39)

With $M_{\mathrm{p}}$, $h_{\mathrm{p}}$ and other parameters not changing in
time, as we assume in this work (but see Section 7.1), $S_{\zeta,\mathrm{sh}}$
remains constant. Then the integral in (39) is trivial, and the time for
vortices to appear is simply $\tau\approx
S_{\zeta,\mathrm{sh}}^{-1}\zeta_{\mathrm{th}}$. Given the expression (38), it
then follows that the characteristic timescale $\tau$ for the emergence of
vortices at the edges of a planetary gap should scale as

$\displaystyle\tau\propto
P_{\mathrm{p}}h_{\mathrm{p}}^{-1}\left(\frac{M_{\mathrm{p}}}{M_{\mathrm{th}}}\right)^{-2.6}.$
(40)

Given the simplicity of our argument, this estimate may apply to both
$t_{\mathrm{lin}}$ and $t_{\mathrm{nl}}$. Regarding the former, equation (33)
implies stronger $t_{\mathrm{lin}}$ dependence on planet mass,
$t_{\mathrm{lin}}\mathrel{\vbox{
\offinterlineskip\halign{\hfil$#$\cr\propto\cr\kern
2.0pt\cr\sim\cr\kern-2.0pt\cr}}}\left(M_{\mathrm{p}}/M_{\mathrm{th}}\right)^{-3.2}$,
than in (40). On the other hand, the value of $\alpha$ in the
$t_{\mathrm{nl}}$ fit (34) is close to $-2.6$ in equation (40). This improved
agreement for $t_{\mathrm{nl}}$ is likely a coincidence, since calculation of
$t^{\mathrm{SA}}_{\mathrm{nl}}$ (leading to $t_{\mathrm{nl}}$) takes into
account the evolution of $\gamma_{m}(t)$ during the linear stage of
instability, while the estimate (40) ignores such details.

Note also that our estimate (40) reproduces reasonably well the explicit
dependence of the instability timescale on $h_{\mathrm{p}}$. Indeed,
$t_{\mathrm{lin}}\mathrel{\vbox{
\offinterlineskip\halign{\hfil$#$\cr\propto\cr\kern
2.0pt\cr\sim\cr\kern-2.0pt\cr}}}h_{\mathrm{p}}^{-1.1}$ and
$t_{\mathrm{nl}}\mathrel{\vbox{
\offinterlineskip\halign{\hfil$#$\cr\propto\cr\kern
2.0pt\cr\sim\cr\kern-2.0pt\cr}}}h_{\mathrm{p}}^{-0.86}$ according to the fits
(33) & (34), which are reasonably close to $\tau\propto h_{\mathrm{p}}^{-1}$
in terms of the slope of scaling.

###  6.4 Vortensity threshold for RWI

Whether the onset of planet-induced RWI can be directly associated with
$\delta\zeta$ exceeding some threshold value $\zeta_{\mathrm{th}}$, as assumed
in the derivation of the scaling (40), is an important question. Existence of
a well-defined $\zeta_{\mathrm{th}}$ would be extremely useful since then one
could predict the emergence of vortices without running the time-consuming
linear RWI analysis (it would still be needed if one were interested in e.g.
the RWI growth rate); instead one would simply follow the evolution of the
vortensity profile which can be easily done using the methods of Paper I. And
it was shown in Ono et al. (2016) that for simple models of localized
vortensity perturbations in discs (e.g. due to a Gaussian bump, or a sharp
increase of $\Sigma$) some statements regarding the RWI triggering in terms of
the amplitude and radial scale of $\delta\Sigma(R)$ features can indeed be
made.

On the other hand, emergence of growing modes in equation (18) is determined
entirely by the behaviour of the potential $D$, given by equation (19). Even
for simple models considered in Ono et al. (2016) $D$ has a very complicated
radial structure, see Figs. 5 & 9 in that work. In the planetary case
$\zeta(R)$ is considerably more complex, see Fig. 3, which gets reflected in
even more sophisticated radial profile of $D$. Moreover, RWI triggering
depends not only on the amplitude but also on the width of vortensity
features, which scales with $l_{\mathrm{sh}}$ and changes as $h_{\mathrm{p}}$
or $M_{\mathrm{p}}$ are varied (also higher temperatures tend to promote RWI
growth, see Li et al., 2000).

Given these complications, to provide a direct test of our
$\zeta_{\mathrm{th}}$ assumption, we measured the extremal values of
$\delta\zeta$ (deviation from $\zeta_{\mathrm{i}}$) at time
$t^{\mathrm{SA}}_{\mathrm{lin}}$ (onset of RWI) in our semi-analytic131313We
verified that extrema of $\delta\zeta$ measured in simulations agree with
these values typically to $\lesssim 10\%$, which should be obvious from Figs.
3 & 14. calculations for different disc and planet parameters. The results are
shown as a function of $l_{\mathrm{sh}}$ in Fig. 12, where we display both the
peak ($\delta\zeta_{\mathrm{max}}$, top) and trough
($\delta\zeta_{\mathrm{min}}$, bottom) values of $\delta\zeta$ on the side of
the gap where the RWI sets in first (inner for $p=0$, outer for $p=3/2$). In
Fig. 13 we plot the same information but at time
$t_{\mathrm{nl}}^{\mathrm{SA}}$, when the RWI becomes non-linear and vortices
should appear.

Figure 12:  Characteristic extremal values of vortensity perturbation (in
units of $\Omega_{\mathrm{p}}/\Sigma_{\mathrm{p}}$) at the onset of linear RWI
(at time $t^{\mathrm{SA}}_{\mathrm{lin}}$), extracted for different disc-
planet parameter sets (indicated with various colours and symbols in the
legend). We show both the highest (top panel) and lowest (bottom panel) values
of the vortensity perturbation at the edge of the gap where the instability
develops first (which is determined by $p$ in this case). This calculation
uses our semi-analytical framework. Figure 13:  Same as Fig. 12, but now
showing the extrema of vortensity perturbation at time
$t_{\mathrm{nl}}^{\mathrm{SA}}$, when the RWI reaches its non-linear stage and
fully developed vortices appear in simulations.

There are several observations that we can make based on these plots. First,
extremal values of $\delta\zeta$ do depend on $l_{\mathrm{sh}}$ — the scale of
vortensity features; at any given $l_{\mathrm{sh}}$ there is a large scatter
in $\delta\zeta$ (a factor of several) as different parameters are varied, but
the envelope clearly increases with $l_{\mathrm{sh}}$. Second, extremal values
of $\delta\zeta$ follow well-defined branches when only one parameter of the
model is varied. In particular, we show in Appendix E that the data shown in
Fig. 12 can be matched using a simple formula (48) with three parameters
depending on $p$. Thus, despite the complexity of the radial behaviour of $D$,
the extrema of $\delta\zeta$ change in a systematic (and not chaotic) fashion.
Third, as one decreases $M_{\mathrm{p}}/M_{\mathrm{th}}$ while keeping
$h_{\mathrm{p}}$ fixed (varying colour of symbols), the extremal
$|\delta\zeta|$ increase, but much faster for $p=0$ than for $p=3/2$. Fourth,
as one increases $h_{\mathrm{p}}$ while keeping
$M_{\mathrm{p}}/M_{\mathrm{th}}$ fixed, the extrema of $|\delta\zeta|$
increase for $p=0$, while slightly decreasing for $p=3/2$. Fifth, the extremal
values of $\delta\zeta$ typically vary much less for $p=3/2$ (a factor of
$\sim 2$) than for $p=0$ (which vary by up to $\sim 6$).

Given these trends, it appears that the stability of the gap edge to RWI is
determined not only by the overall amplitude of the vortensity perturbation
but also by other factors, e.g. the radial scale and overall shape of
$\delta\zeta$. For that reason, it is also not surprising that the timescale
fit parameters (33) and (34) show some deviations from the predictions of our
simple estimate (40). What is remarkable, is that these deviations are rather
small and are almost independent of the surface density slope $p$, especially
given the clearly $p$-dependent variation of the critical $\delta\zeta$ values
(by a factor of several) and the diverse trends with $h_{\mathrm{p}}$ and
$M_{\mathrm{p}}$ in Figs. 12 & 13.

The larger amplitude of extremal $\delta\zeta$ for $p=0$ can be explained by
the fact that in this case vortices first appear at the inner edge of the
planetary gap. The rapid variation of $|\Omega-\Omega_{\mathrm{p}}|$ in the
inner disc (see Fig. 11) tends to broaden the $\delta\zeta(R)$ distribution,
requiring higher values of $\delta\zeta$ for RWI to set in when
$l_{\mathrm{sh}}$ is large. But this also explains why the instability
timescale stays roughly the same regardless of $p$: for large
$l_{\mathrm{sh}}$ the rate at which the fluid passes through the shock in the
inner disc can be considerably higher than predicted by the local
approximation (37), see Fig. 11. Our results imply that this effect almost
fully compensates for the increased $|\zeta_{\mathrm{th}}|$ in $p=0$ discs,
resulting in $t_{\mathrm{lin}}$ and $t_{\mathrm{nl}}$ being essentially
independent of $p$.

To summarize, we do not find that, in general, $\zeta_{\mathrm{sh}}$ is
confined within a well-defined narrow range, as assumed in Section 6.3. This
is approximately true for $p=3/2$, where we observe variations of
$\zeta_{\mathrm{th}}$ by a factor of only $\sim 2$, but for $p=0$ critical
$\zeta_{\mathrm{sh}}$ may vary much more. This is likely good news for
interpreting observations, since we expect higher $p$ (i.e. radially
decreasing $\Sigma$) to be more typical for PPDs. On the other hand, for any
$p$ we do observe that the characteristic values of $|\zeta_{\mathrm{th}}|$
are of order unity, which is reassuring and may explain the insensitivity of
both $t_{\mathrm{lin}}$ and $t_{\mathrm{nl}}$ to $p$, and their rough
agreement with the scaling (40), which we found in Sections 6.1 & 6.2. We
leave detailed exploration of the trends observed in Figs. 12 & 13 and 17 to
future work.

##  7 Discussion

Semi-analytic calculation of the timescales $t^{\mathrm{SA}}_{\mathrm{lin}}$
and $t_{\mathrm{nl}}^{\mathrm{SA}}$ (see Section 5.1) is a multi-step
procedure: we first need to compute the evolution of the vortensity profile,
then reconstruct $\Sigma,\Omega$ profiles (Section 4.2), after which we
perform the RWI stability analysis (Section 4.1) and integrate growth rates to
obtain $t_{\mathrm{nl}}^{\mathrm{SA}}$ (Section 5.1). Despite this complexity,
$t_{\mathrm{nl}}^{\mathrm{SA}}$ agrees remarkably well with the simulation-
based $\hat{t}^{\mathrm{sim}}_{\mathrm{nl}}$, providing strong support for the
robustness of our semi-analytical framework. We now discuss some ways in which
one can capitalize on this agreement, and put our results in context with
existing studies.

###  7.1 Applications of our results

Our results (34)-(36) for the planet-induced vortex development timescales can
be used for interpreting observations of protoplanetary discs showing evidence
of non-axisymmetric features (arcs, clumps, etc.), which may be interpreted as
vortices. Without going into too many details here, we note that equation (36)
can be used to constrain either the mass or the age of a planet responsible
for the appearance of an observed vortex in a protoplanetary disc. We provide
more details on this in Rafikov & Cimerman (in prep.).

Our timescale estimates (31)-(34) are also useful for planning future
numerical studies of RWI and vortices in protoplanetary discs. They allow one
to estimate the run time of the simulations producing vortices in advance,
facilitating making informed decisions about allocation of computing
resources. Our timescale fits (31)-(34) can also be used for benchmarking
numerical codes.

Our $\Sigma$, $\Omega$ reconstruction technique (see Section 4.2) based on Lin
& Papaloizou (2011) allows one to construct radial profile of the planet-
induced gap given the knowledge of the radial vortensity profile. Once
combined with our recipe for vortensity production (valid for
$M_{\mathrm{p}}\lesssim M_{\mathrm{th}}$) it provides a fully self-contained
semi-analytical method for predicting the radial profiles of shallow gaps at
arbitrary moments of time, which we thoroughly test against simulations (see
Section 4.2.2 and Appendix C). This is a powerful tool, complementary to the
existing methods for construction of profiles of deep gaps carved out by
massive ($M_{\mathrm{p}}\gtrsim M_{\mathrm{th}}$) planets (Crida et al., 2006;
Kanagawa et al., 2017; Ginzburg & Sari, 2018; Duffell, 2020), that can find
multiple uses.

In particular, radial profiles of $\Sigma$, which can be easily computed as a
function of time using our method, could be used to study dust trapping at
pressure bumps (e.g. Stammler et al., 2019). Especially at the outer gap edge,
where our semi-analytical model matches simulations well, we expect modelling
of dust accumulation to be quite accurate (in the inner disc multiple spirals
may lead to some discrepancies with simulations). This provides an efficient
and accurate way of studying e.g. the dependence of pebble isolation mass
(Paardekooper & Mellema, 2004; Lambrechts et al., 2014) on disc properties
such as $h_{\mathrm{p}}$ and $p$ (in the inviscid limit).

Numerical studies of RWI development, dust accumulation, etc. near the low
mass (sub-thermal) planets typically need to be run for thousands of orbits to
properly capture the disc structure at late times. Our semi-analytical method
provides a useful shortcut: one can use it to construct the approximate disc
state (including the forming gap) around the planet at any moment of time
after the planet has been introduced in the simulation. This state would
provide the initial condition for running the simulation, avoiding the initial
‘burn-in’ stage needed to reach this situation in the simulation.

We should also note that our semi-analytical method can naturally account for
various processes that have been neglected in this study — planet migration,
history of its mass accumulation (Hammer et al., 2017; Hallam & Paardekooper,
2020), evolution of the disc state on long timescales, etc. To enable this,
one simply needs to generalize equation (20) by allowing $\Omega_{\mathrm{p}}$
and $\Delta\zeta(R)$ to be functions of time set by the evolution of the
physical characteristics of the problem at hand, e.g. planetary semi-major
axis $R_{\mathrm{p}}$, mass $M_{\mathrm{p}}$, etc. This equation can then be
easily integrated, obtaining a more general $\zeta(R,t)$ than the solution
(21) used in this work141414This implicitly assumes that the time-variation of
$\Delta\zeta$ is slow compared to the dynamical time (which is typically
satisfied, e.g. since our $t_{\mathrm{lin}}\gg P_{\mathrm{p}}$, see Figure 9),
such that the radial force balance, i.e. equation (3), holds.. Using this
$\zeta(R,t)$ to reconstruct the evolving $\Sigma$, $\Omega$ profiles (see
Section 4.2) and running the RWI linear stability analysis (see Sections 4.1 &
4.3) with these inputs, one can determine the onset of RWI in these more
general situations.

For illustration, let us sketch how our heuristic determination of the vortex
emergence timescale $\tau$ (see Section 6.3) would change if $M_{\mathrm{p}}$
were not constant but was actually growing in time due to accretion (Hallam &
Paardekooper, 2020). With time-dependent $M_{\mathrm{p}}(t)$ the condition
(39) would become
$\zeta_{\mathrm{th}}\approx\int^{\tau}_{0}S_{\zeta,\mathrm{sh}}(t)\mathrm{d}t\propto\int^{\tau}_{0}M_{\mathrm{p}}^{2.6}(t)\mathrm{d}t$,
see equation (38), providing a new relation for determining $\tau$ instead of
equation (40). Because of the steep dependence on $M_{\mathrm{p}}(t)$ in this
integral we expect that the determination of $\tau$ would be most sensitive to
the details of the $M_{\mathrm{p}}(t)$ behaviour at late time, when the planet
is close to reaching its final mass.

###  7.2 Potential limitations of this work

In this work, we made several simplifying assumptions to highlight the most
important physical processes leading to vortex formation. In particular, we
neglected the presence of dust in the disc and its back-reaction on the fluid
motions and ignored the disc self-gravity (Lin, 2012b). Unlike Hammer et al.
(2017) and Hallam & Paardekooper (2020), we keep the planetary mass fixed in
time since accounting for the planetary accretion history would necessarily
introduce some ad hoc assumptions about the $M_{\mathrm{p}}(t)$ dependence.
Similarly, we fix the semi-major axis of the planetary orbit not allowing it
to migrate and set planetary eccentricity to zero. We also allow only one
planet to be present in the disc. As mentioned in Section 7.1, our semi-
analytical approach outlined in Section 4.2 can be easily generalized to
account for the history of planetary mass accretion and migration, multiple
planets (Garrido-Deutelmoser et al., 2022), and so on. We restrict the disc to
be strictly two-dimensional. However, in their study of linear RWI in a 3D
disc, Lin (2012a, b) found that its growth rates can be accurately predicted
from the 2D problem alone. Thus, introducing the third dimension is unlikely
to change our main conclusions.

####  7.2.1 Non-barotropic thermodynamics

In line with our previous work (Paper I), all our calculations assume a
globally isothermal EoS, which is a special case of a barotropic EoS
($P=P(\Sigma)$ only). In real discs there are various heating and cooling
processes that would cause departures from the barotropic setup. Adoption of
more sophisticated thermodynamic assumptions, even as simple as the
introduction of $\beta$-cooling on a timescale
$\beta\Omega_{\mathrm{K}}^{-1}$, would have several consequences.

First, the excitation of planetary density waves would change, see Miranda &
Rafikov (2020a) who explored the dependence of planetary torques on $\beta$.
Second, (linear) thermal relaxation can dramatically modify damping of the
spiral waves as shown in Miranda & Rafikov (2020a, b). This must affect radial
distribution of the vortensity production around the planet and the resulting
gap profiles. Third, baroclinic effects may start affecting vortensity
evolution. Fourth, irreversible heating of the disc by planet-driven spiral
shocks (e.g. Rafikov, 2016; Ziampras et al., 2020) can modify the local
temperature profile near the planet, further impacting density wave
propagation and potentially driving additional baroclinic effects. Fifth, a
specific form of the EoS affects the vortensity generation at the shock front,
e.g. compare Lin & Papaloizou (2010) and Paper I. Finally, even for identical
surface density perturbations (i.e. neglecting the aforementioned effects),
the linear growth of the RWI is affected by thermal relaxation as shown by
e.g. Les & Lin (2015) and Huang & Yu (2022).

Thermal relaxation also affects the long-term survival of vortices (Fung &
Ono, 2021; Rometsch et al., 2021), but this issue is beyond the scope of our
study, which focuses on vortex generation.

####  7.2.2 The effect of shear viscosity

In line with Paper I, we assumed the disc to be inviscid. This assumption is
supported by many observations suggesting that viscosity is likely low in most
PPDs (Pinte et al., 2016; Rafikov, 2017; Flaherty et al., 2020). Nevertheless,
it is still important to assess the consequences of disc viscosity being non-
zero.

Miranda & Rafikov (2020a) have shown that, unless the disc is very viscous
(with effective $\alpha\gtrsim 0.01$), viscous stresses do not affect density
wave dissipation, which is determined by the combination of non-linear and
radiative damping. However, another effect of viscosity is to diffusively
smooth out any features in the radial vortensity distribution. In particular,
it introduces the diffusive term in the vortensity evolution equation, which
will tend to smear out the sharp peaks and troughs of $\zeta$ produced in the
vicinity of the planet, see Figs. 3 & 14. As a result, viscosity would slow
down the growth of $\zeta$ and could delay the onset of RWI as has been shown
in e.g. Hallam & Paardekooper (2020). Thus, our inviscid calculations provide
us with the lower limit for the time to reach instability, and $\tau$ will be
longer in sufficiently viscous discs.

High viscosity may suppress the RWI entirely (Hammer et al., 2017), by
preventing $\zeta$ from reaching the threshold value $\zeta_{\mathrm{th}}$
necessary for instability to set in. We can provide a simple heuristic
estimate of when this might happen. Without invoking the explicit form of the
viscous term in the vortensity evolution equation, we can model it as a
diffusion process151515It should be remembered that the real stress in the
disc, for example caused by the MRI, may respond to planetary torques
differently from a simple shear viscosity model (e.g. Zhu et al., 2013). with
a characteristic diffusion coefficient equal to the kinematic viscosity $\nu$
(similar to the momentum equation), i.e.
$\partial_{t}\zeta|_{\nu}\sim\nu\nabla^{2}\zeta$. Since $\zeta$ varies on
scales $\sim l_{\mathrm{sh}}$, this term should be of order
$\nu\zeta/l_{\mathrm{sh}}^{2}$. Viscous diffusion would stop the planet-driven
growth of vortensity at some equilibrium value $\zeta_{\mathrm{eq}}$ at which
$\partial_{t}\zeta|_{\nu}$ would match $S_{\zeta,\mathrm{sh}}$ given by
equation (38). In other words,

$\displaystyle\Omega_{\mathrm{p}}h_{\mathrm{p}}\left(\frac{M_{\mathrm{p}}}{M_{\mathrm{th}}}\right)^{2.6}\propto\nu\frac{\zeta_{\mathrm{eq}}}{l_{\mathrm{sh}}^{2}},$
(41)

which, using equation (6) and $\nu=\alpha
h_{\mathrm{p}}^{2}\Omega_{\mathrm{p}}R_{\mathrm{p}}^{2}$, can be written as

$\displaystyle\zeta_{\mathrm{eq}}\propto\alpha^{-1}h_{\mathrm{p}}\left(\frac{M_{\mathrm{p}}}{M_{\mathrm{th}}}\right)^{1.8}.$
(42)

RWI would be suppressed altogether and vortices would not emerge if
$\zeta_{\mathrm{eq}}<\zeta_{\mathrm{th}}$. In particular, equation (42)
implies that for this to be the case planetary mass must satisfy

$\displaystyle\frac{M_{\mathrm{p}}}{M_{\mathrm{th}}}<B\left(\frac{\alpha}{h_{\mathrm{p}}}\right)^{5/9},$
(43)

where $B$ is a (dimensionless) constant which can be calibrated using
simulations. This relation is necessarily approximate, with true power law
exponents possibly deviating somewhat from our predictions, which again can be
checked using simulations. We note that McNally et al. (2019) found a
different relation between $M_{\mathrm{p}}$ and $\alpha$ (with the opposite
sign of the exponent), which separates evolution with and without vortices.
The reason for this difference is that (unlike us) they also considered planet
migration, which can lead to non-trivial feedback effects in the presence of
viscosity, see Rafikov (2002b). We also cannot directly compare our relation
(43) with the results of Hallam & Paardekooper (2020), since planet mass was
varied in their study. In light of our results, formation of vortices in
viscous discs clearly warrants further investigation.

###  7.3 Comparison with previous works

A number of past studies explored the emergence of vortices at the edges of
planetary gaps, and some of their approaches share similarities with our work.
For example, in their study of high mass planets (their lowest mass is
$M_{\mathrm{p}}=0.8M_{\mathrm{th}}$), de Val-Borro et al. (2007) investigated
linear stability to RWI of a simulated disc perturbed by a planet in a time-
dependent manner, like we do when determining
$t^{\mathrm{sim}}_{\mathrm{lin}}$ and $t^{\mathrm{sim}}_{\mathrm{nl}}$. They
found an increase in the growth rate of unstable modes with time (see their
Fig. 5) as the planet opens a deeper and deeper gap, similar to what we show
in our Fig. 6. Adopting a constant $\Sigma_{\mathrm{i}}$ (i.e. $p=0$) disc,
they found higher growth rates at the inner gap edge, in line with our
findings. While they speculated that this is caused by the proximity to domain
boundaries, our analysis shows that this is a real effect due to the
accelerated steepening of the inner spiral arm for $p=0$ (Paper I).

Lin & Papaloizou (2010) carried out a similar analysis, also focusing on
higher mass planets ($M_{\mathrm{p}}\geq M_{\mathrm{th}}$). Their work was the
first to introduce $\Sigma,\Omega$ reconstruction from vortensity $\zeta$, and
they found good performance of this technique using $\zeta$ profiles from
simulations. We extend this method even further by employing a semi-analytic
prescription for $\zeta$ evolution from Paper I, see Section 4.2.1, and
confirm its accuracy for a range of relevant parameters, see Section 4.2.2 and
Appendix C.

A number of studies (e.g. Li et al., 2009; Yu et al., 2010; McNally et al.,
2019) numerically studied emergence of vortices in viscous discs with
migrating planets. Because of these additional physical ingredients, which
introduce new effects (e.g. gas redistribution associated with migration,
Rafikov, 2002a) we cannot provide a direct comparison of our results with
their findings. This also precludes us from comparing our calculations of the
vortex emergence timescale $\tau$ with the works of Hammer et al. (2017) and
Hallam & Paardekooper (2020), since they considered time-varying planet mass
and viscous discs (although we discuss ways in which the impact of these
additional physical ingredients might be understood in Sections 7.1 and 7.2.2,
respectively).

Finally, we note the similarity of planet-driven vortex production in PPDs
with the origin of vortices found by Coleman et al. (2022) in their inviscid
2D simulations of accretion disc boundary layers. The only key difference is
that in their case density waves are driven not by planets but by the acoustic
instability in the boundary layer (Belyaev & Rafikov, 2012; Belyaev et al.,
2012, 2013), while all subsequent processes — production of vortensity at the
shock fronts, triggering of RWI, mergers of the resultant vortices — are
essentially identical.

##  8 Summary

We studied the stability of gaps carved by (sub-thermal mass, fixed in time)
planets to RWI in globally isothermal, inviscid, 2D protoplanetary discs. Our
primary goals were to study the disc evolution leading to instability and to
determine the timescale on which RWI sets in and vortices form at the edges of
planet-driven gaps. We used a two-stage approach to reach these goals. First,
we developed a closed-form semi-analytical approach for computing the radial
profile of a gap induced by a sub-thermal mass planet at any moment of time
given a set of disc and planetary parameters (Section 4.2). This was achieved
by coupling the semi-analytical calculation of the vortensity production at
the planetary shock (Paper I) with the $\Sigma,\Omega$-reconstruction
technique (Lin & Papaloizou, 2010). Second, we carried out a linear RWI
stability analysis (Section 4.1) on these gap profiles to determine the onset
of instability. Along the way, our results have been verified at all levels
against direct hydro simulations. Below we briefly summarize our main
findings.

  * •

We showed (Section 4.2.2 and Appendix C) that our semi-analytical method
reliably reproduces gap profiles when compared to expensive 2D simulations for
a mass range $M_{\mathrm{p}}=(0.1-0.5)M_{\mathrm{th}}$ (with somewhat reduced
agreement in the inner disc caused by the emergence and dissipation of the
secondary spiral arm). This allowed us to study the RWI stability of gap edges
in a fast and efficient manner.

  * •

We calculated the timescales for the onset of linear RWI and for the
development of its non-linear phase. We provided useful fitting formulae for
these timescales (Sections 6.1,6.2) valid in the $M_{\mathrm{p}}\lesssim
M_{\mathrm{th}}$ regime and proposed an approximate, heuristic explanation of
these scalings (Section 6.3).

  * •

Consistent with previous studies, we found that hotter discs (keeping
$M_{\mathrm{p}}/M_{\mathrm{th}}$ constant) and higher mass planets accelerate
the development of RWI.

  * •

In agreement with Paper I, we found that the disc surface density slope $p$
controls, which side of the gap develops RWI first: inner gap edge for
constant-$\Sigma$ discs ($p=0$), but outer edge for discs with radially
decreasing $\Sigma$ ($p=3/2$).

  * •

We found that at the point when RWI starts developing, the planet-induced
vortensity deviation $\delta\zeta$ shows considerable variation as disc and
planetary parameters are changed (Section 6.4). While we observe certain
trends in the behaviour of this characteristic $\delta\zeta$, its high
sensitivity to the value of the surface density slope $p$ precludes us from
drawing simple and universal conclusions.

  * •

We have shown how our semi-analytical approach can be extended beyond our
simple setup, for example, to account for planetary accretion and migration
(Section 7.1), and assessed the impact of non-zero disc viscosity on our
results (Section 7.2.2).

  * •

Our findings can be used for interpreting observations of PPDs with vortex-
like structures that may be caused by planets, as well as for setting up and
testing simulations of disc-planet interaction.

In Rafikov & Cimerman (in prep.) we apply the results of this study to
constrain masses and ages of putative planets in protoplanetary discs with
vortex-like structures.

## Acknowledgements

Software: NumPy (Harris et al., 2020), SciPy (Virtanen et al., 2020), IPython
(Perez & Granger, 2007), Matplotlib (Hunter, 2007), Athena++ (Stone et al.,
2020). We thank the referee Takayuki Muto for a constructive report which
helped us in clarifying several points made in this work. N.P.C. would like to
thank Robin Croft for helpful discussions regarding the relaxation solver and
all developers of Athena++ for making their code publicly available. N.P.C. is
funded by an Isaac Newton Studentship and a Science and Technology Facilities
Council (STFC) studentship. R.R.R. acknowledges financial support through the
NASA grant 15-XRP15-2-0139, Ambrose Monell Foundation, and STFC grant
ST/T00049X/1. A large part of the long term simulations were performed on the
HPC cluster FAWCETT at DAMTP, University of Cambridge. Part of this work was
performed using resources provided by the Cambridge Service for Data Driven
Discovery (CSD3) operated by the University of Cambridge Research Computing
Service (www.csd3.cam.ac.uk), provided by Dell EMC and Intel using Tier-2
funding from the Engineering and Physical Sciences Research Council (capital
grant EP/P020259/1), and DiRAC funding from the Science and Technology
Facilities Council (www.dirac.ac.uk).

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Andrews (2020) Andrews S. M., 2020, ARA&A, 58, 483
  * Bae & Zhu (2018) Bae J., Zhu Z., 2018, ApJ, 859, 118
  * Bae et al. (2017) Bae J., Zhu Z., Hartmann L., 2017, ApJ, 850, 201
  * Baruteau et al. (2019) Baruteau C., et al., 2019, MNRAS, 486, 304
  * Belyaev & Rafikov (2012) Belyaev M. A., Rafikov R. R., 2012, ApJ, 752, 115
  * Belyaev et al. (2012) Belyaev M. A., Rafikov R. R., Stone J. M., 2012, ApJ, 760, 22
  * Belyaev et al. (2013) Belyaev M. A., Rafikov R. R., Stone J. M., 2013, ApJ, 770, 67
  * Boehler et al. (2021) Boehler Y., et al., 2021, A&A, 650, A59
  * Chandrasekhar (1961) Chandrasekhar S., 1961, Hydrodynamic and Hydromagnetic Stability. Dover Books on Physics, Dover Publications
  * Cimerman & Rafikov (2021) Cimerman N. P., Rafikov R. R., 2021, MNRAS, 508, 2329
  * Coleman et al. (2022) Coleman M. S. B., Rafikov R. R., Philippov A. A., 2022, MNRAS, 509, 440
  * Crida et al. (2006) Crida A., Morbidelli A., Masset F., 2006, Icarus, 181, 587
  * Dong et al. (2011a) Dong R., Rafikov R. R., Stone J. M., Petrovich C., 2011a, ApJ, 741, 56
  * Dong et al. (2011b) Dong R., Rafikov R. R., Stone J. M., 2011b, ApJ, 741, 57
  * Dong et al. (2017) Dong R., Li S., Chiang E., Li H., 2017, ApJ, 843, 127
  * Dong et al. (2018) Dong R., et al., 2018, ApJ, 860, 124
  * Duffell (2020) Duffell P. C., 2020, ApJ, 889, 16
  * Flaherty et al. (2020) Flaherty K., et al., 2020, ApJ, 895, 109
  * Fung & Ono (2021) Fung J., Ono T., 2021, ApJ, 922, 13
  * Garrido-Deutelmoser et al. (2022) Garrido-Deutelmoser J., Petrovich C., Krapp L., Kratter K. M., Dong R., 2022, ApJ, 932, 41
  * Ginzburg & Sari (2018) Ginzburg S., Sari R., 2018, MNRAS, 479, 1986
  * Goodman & Rafikov (2001) Goodman J., Rafikov R. R., 2001, ApJ, 552, 793
  * Hallam & Paardekooper (2020) Hallam P. D., Paardekooper S. J., 2020, MNRAS, 491, 5759
  * Hammer et al. (2017) Hammer M., Kratter K. M., Lin M.-K., 2017, MNRAS, 466, 3533
  * Hammer et al. (2021) Hammer M., Lin M.-K., Kratter K. M., Pinilla P., 2021, MNRAS, 504, 3963
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357–362
  * Huang & Yu (2022) Huang S., Yu C., 2022, arXiv e-prints, p. arXiv:2205.11969
  * Hunter (2007) Hunter J. D., 2007, Computing in Science Engineering, 9, 90
  * Kanagawa et al. (2017) Kanagawa K. D., Tanaka H., Muto T., Tanigawa T., 2017, PASJ, 69, 97
  * Kevlahan (1997) Kevlahan N. K. R., 1997, Journal of Fluid Mechanics, 341, 371
  * Kraus et al. (2017) Kraus S., et al., 2017, ApJ, 848, L11
  * Lambrechts et al. (2014) Lambrechts M., Johansen A., Morbidelli A., 2014, A&A, 572, A35
  * Les & Lin (2015) Les R., Lin M.-K., 2015, MNRAS, 450, 1503
  * Lesur & Papaloizou (2009) Lesur G., Papaloizou J. C. B., 2009, A&A, 498, 1
  * Li et al. (2000) Li H., Finn J. M., Lovelace R. V. E., Colgate S. A., 2000, ApJ, 533, 1023
  * Li et al. (2009) Li H., Lubow S. H., Li S., Lin D. N. C., 2009, ApJ, 690, L52
  * Lin (2012a) Lin M.-K., 2012a, MNRAS, 426, 3211
  * Lin (2012b) Lin M.-K., 2012b, ApJ, 754, 21
  * Lin & Papaloizou (2010) Lin M.-K., Papaloizou J. C. B., 2010, MNRAS, 405, 1473
  * Lin & Papaloizou (2011) Lin M.-K., Papaloizou J. C. B., 2011, MNRAS, 415, 1445
  * Lovelace et al. (1999) Lovelace R. V. E., Li H., Colgate S. A., Nelson A. F., 1999, ApJ, 513, 805
  * Marr & Dong (2022) Marr M., Dong R., 2022, arXiv e-prints, p. arXiv:2203.11953
  * McNally et al. (2019) McNally C. P., Nelson R. P., Paardekooper S.-J., Benítez-Llambay P., 2019, MNRAS, 484, 728
  * Meheut et al. (2010) Meheut H., Casse F., Varniere P., Tagger M., 2010, A&A, 516, A31
  * Meheut et al. (2012) Meheut H., Meliani Z., Varniere P., Benz W., 2012, A&A, 545, A134
  * Miranda & Rafikov (2019a) Miranda R., Rafikov R. R., 2019a, ApJ, 875, 37
  * Miranda & Rafikov (2019b) Miranda R., Rafikov R. R., 2019b, ApJ, 878, L9
  * Miranda & Rafikov (2020a) Miranda R., Rafikov R. R., 2020a, ApJ, 892, 65
  * Miranda & Rafikov (2020b) Miranda R., Rafikov R. R., 2020b, ApJ, 904, 121
  * Muller (1956) Muller D. E., 1956, Mathematical Tables and Other Aids to Computation, 10, 208
  * Muto et al. (2010) Muto T., Suzuki T. K., Inutsuka S.-i., 2010, ApJ, 724, 448
  * Ono et al. (2016) Ono T., Muto T., Takeuchi T., Nomura H., 2016, ApJ, 823, 84
  * Ono et al. (2018) Ono T., Muto T., Tomida K., Zhu Z., 2018, ApJ, 864, 70
  * Paardekooper & Mellema (2004) Paardekooper S. J., Mellema G., 2004, A&A, 425, L9
  * Papaloizou & Lin (1989) Papaloizou J. C. B., Lin D. N. C., 1989, ApJ, 344, 645
  * Perez & Granger (2007) Perez F., Granger B. E., 2007, Computing in Science Engineering, 9, 21
  * Pérez et al. (2018) Pérez L. M., et al., 2018, ApJ, 869, L50
  * Petrovich & Rafikov (2012) Petrovich C., Rafikov R. R., 2012, ApJ, 758, 33
  * Pinte et al. (2016) Pinte C., Dent W. R. F., Ménard F., Hales A., Hill T., Cortes P., de Gregorio-Monsalvo I., 2016, ApJ, 816, 25
  * Rafikov (2002a) Rafikov R. R., 2002a, ApJ, 569, 997
  * Rafikov (2002b) Rafikov R. R., 2002b, ApJ, 572, 566
  * Rafikov (2016) Rafikov R. R., 2016, ApJ, 831, 122
  * Rafikov (2017) Rafikov R. R., 2017, ApJ, 837, 163
  * Rometsch et al. (2021) Rometsch T., Ziampras A., Kley W., Béthune W., 2021, A&A, 656, A130
  * Stammler et al. (2019) Stammler S. M., Drążkowska J., Birnstiel T., Klahr H., Dullemond C. P., Andrews S. M., 2019, ApJ, 884, L5
  * Stone et al. (2020) Stone J. M., Tomida K., White C. J., Felker K. G., 2020, ApJS, 249, 4
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Yu et al. (2010) Yu C., Li H., Li S., Lubow S. H., Lin D. N. C., 2010, ApJ, 712, 198
  * Zhang et al. (2018) Zhang S., et al., 2018, ApJ, 869, L47
  * Zhu & Stone (2014) Zhu Z., Stone J. M., 2014, ApJ, 795, 53
  * Zhu et al. (2013) Zhu Z., Stone J. M., Rafikov R. R., 2013, ApJ, 768, 143
  * Zhu et al. (2014) Zhu Z., Stone J. M., Rafikov R. R., Bai X.-n., 2014, ApJ, 785, 122
  * Zhu et al. (2015) Zhu Z., Dong R., Stone J. M., Rafikov R. R., 2015, ApJ, 813, 88
  * Ziampras et al. (2020) Ziampras A., Ataiee S., Kley W., Dullemond C. P., Baruteau C., 2020, A&A, 633, A29
  * de Val-Borro et al. (2007) de Val-Borro M., Artymowicz P., D’Angelo G., Peplinski A., 2007, A&A, 471, 1043
  * van der Marel et al. (2013) van der Marel N., et al., 2013, Science, 340, 1199
  * van der Marel et al. (2016) van der Marel N., Cazzoletti P., Pinilla P., Garufi A., 2016, ApJ, 832, 178

##  Appendix A Search for RWI modes

We discretize the differential operators in equation (9) using $N_{R}$ cells
in radius, which allows us to write this equation as

$\displaystyle W(R,\omega,m)\Psi_{m}=0,$ (44)

where $W$ is a complex, tri-diagonal square matrix of size $N_{R}\times
N_{R}$. Existence of non-trivial solutions $\Psi_{m}$ thus requires $\det W=0$
or, in other words, that $W$ has at least one eigenvalue $\lambda_{i}=0$. The
problem is then reduced to finding the values of $\omega$ giving $\det W=0$.
For most parameter sets, we solve this problem on a subdomain of the same
logarithmic radial grid that is used in the hydrodynamical models, ranging
from $R_{\mathrm{min}}=0.3R_{\mathrm{p}}$ to
$R_{\mathrm{max}}=3.25R_{\mathrm{p}}$. For some parameter sets (i.e. small
$M_{\mathrm{p}}$ and large $h_{\mathrm{p}}$), this radial domain is increased
to ensure that boundaries are far enough from the gap edges.

In order to allow for wave propagation across boundaries, we use outgoing wave
boundary conditions under the WKB approximation (Ono et al., 2016). While
low-$m$ modes ($m$ $\lesssim 6$) are very localized and insensitive to
boundary conditions, modes of higher $m$ ($\gtrsim 6$) are more global, such
that BCs become important (e.g. Lin & Papaloizou, 2010).

We use an eigenvalue search strategy similar to Li et al. (2000), de Val-Borro
et al. (2007) and Ono et al. (2016), which we have implemented using Python3.
Knowing the radial location of the vortensity minimum $R_{0}$, we use the
result of previous studies that trapped, unstable modes are corotation modes,
i.e. we start our search around $\omega_{\mathrm{R}}\simeq m\Omega(R_{0})$
(and $|\Delta\omega|^{2}\ll\Omega^{2}$). We thus restrict our search to a
region around corotation, typically
$0.9<\omega_{\mathrm{R}}/(m\Omega(R_{0}))<1.1$ and begin by looking for modes
with low growth rates $0.01<2\pi\gamma_{m}<0.4$. Poles and roots of $\det W$
are found via contour integration in the complex $\omega_{\mathrm{R}}-\gamma$
plane as described in de Val-Borro et al. (2007). When a contour contains a
root, we use Muller’s method (Muller, 1956) for locating it more precisely.
Once a mode is found, $\Psi_{m}$ is obtained by solving the eigenvalue problem
(44) using the SciPy routine scipy.sparse.linalg.eigs for sparse matrices
(since $W$ is tri-diagonal). Especially for modes with low $\gamma_{m}$, the
solver sometimes finds spurious unphysical modes, that do not represent
trapped modes. Thus, our method requires some human interaction, but could be
extended to deal with such spurious solutions. Having found an unstable mode
for the disc structure at a time $t$, we can use the complex frequency as an
initial guess in Muller’s method for the next search at time $t+\Delta t$, as
the disc evolves slowly (Li et al., 2000). This saves the more expensive
contour integration time.

We have verified our method by confirming that we find the same modes and
growth rates as Ono et al. (2016) for their fiducial ‘Gaussian Bump’ setup. As
another test, we have run several simulations in which we removed planet after
the disc reached an unstable state. Re-initializing the disc with the
azimuthally averaged $\langle\Sigma(R)\rangle_{\phi}$,
$\langle\Omega(R)\rangle_{\phi}$ from the last frame with the planet, we
monitored the development of RWI with a well-defined (constant in time) growth
rate. In parallel, we also carried out linear RWI analysis on the same
background disc state. We found very good agreement (at the level of a few
percent) of both the growth rates and the radial structure of the mode
measured in simulations with linear theory. This experiment not only verifies
our linear RWI analysis but also provides a very useful test of the code, in
particular, of the orbital advection module, now implemented in Athena++,
which has not been used by Ono et al. (2018).

Figure 14:  Same as Fig. 3a,c,d,e, but for different sets of parameters and
times (what is different compared to the fiducial set
$M_{\mathrm{p}}=0.25M_{\mathrm{th}};h_{\mathrm{p}}=0.05;p=3/2$ is indicated
above each column), comparing near-gap $\zeta$, $\Sigma$, and $\Omega$
profiles from simulations and the semi-analytical method. See Appendix C for
details.

##  Appendix B $\Sigma,\Omega$ reconstruction

To solve their version of equation (23), Lin & Papaloizou (2010) adopted a
shooting method. However, we found that it fails if the radial region of
integration is too large, as only exponentially growing or decaying solutions
are found.

For this reason, we developed a more robust iterative method for solving this
non-linear differential equation: we used a discretized version of equation
(23),

$\displaystyle\frac{\Sigma_{i+1}-2\Sigma_{i}+\Sigma_{i-1}}{w^{2}}=F_{i}(\Sigma,\zeta,R),$
(45)

where $w\equiv R_{i+1}-R_{i}$ is a constant grid spacing and subscript $i$
indicates evaluation at $R_{i}$ and

$\displaystyle F_{i}(\Sigma,\zeta,R)=$
$\displaystyle\frac{2\zeta\Sigma_{i}^{2}}{c_{\mathrm{s}}^{2}}\left(\Omega_{\mathrm{K}}^{2}+\frac{c_{\mathrm{s}}^{2}\left(\Sigma_{i+1}-\Sigma_{i-1}\right)}{2wR\Sigma_{i}}\right)^{1/2}$
(46)
$\displaystyle-\frac{\Omega_{\mathrm{K}}^{2}}{c_{\mathrm{s}}^{2}}\Sigma_{i}+\frac{\left(\Sigma_{i+1}-\Sigma_{i-1}\right){}^{2}}{4w^{2}\Sigma_{i}}-\frac{3\left(\Sigma_{i+1}-\Sigma_{i-1}\right)}{2wR}.$

This can be rearranged to give an iterative scheme for $\Sigma_{i}$:

$\displaystyle\Sigma_{i}^{(n+1)}=(1-C)\Sigma_{i}^{(n)}+\frac{C}{2}\left(-w^{2}F_{i}+\Sigma_{i+1}+\Sigma_{i-1}\right)^{(n)},$
(47)

where the superscript indicates the number of iteration. We typically chose
$C=2/3$. If required, adopting a multi-grid method and making $C$ a function
of iteration number for successive over-relaxation (SOR) could accelerate
convergence.

##  Appendix C Further tests of the semi-analytical reconstruction

Here we describe further tests of our semi-analytical gap reconstruction
technique (Section 4.2), for disc and planetary parameters different from the
fiducial setup.

###  C.1 Variation of the surface density slope $p$

First, we test reconstruction in a constant surface density disc ($p=0$). In
column A of Fig. 14, we show the disc structure at $t=650P_{\mathrm{p}}$ for
our standard $M_{\mathrm{p}}=0.25M_{\mathrm{th}}$ and $h_{\mathrm{p}}=0.05$
but a constant $\Sigma_{\mathrm{i}}$, i.e. $p=0$. As pointed out in Paper I,
this disc model has a non-zero background vortensity gradient (compared to the
fiducial $p=3/2$ disc), such that $\zeta_{\mathrm{i}}$ is no longer radially
constant. The most important effect of the density slope $p$ is that it
controls the part of the disc (inner or outer), in which the vortensity jump
at the shock $\Delta\zeta$ is greater, as the non-linearity of the spiral wake
is modified (e.g. ; , ). Accordingly, the inner gap edge experiences faster
vortensity evolution (and becomes unstable earlier) than the outer one in a
$p=0$ disc, opposite to the $p=3/2$ case covered in Section 4.2.1. This is
what we see in Fig. 14a.

As in the fiducial case, vortensity and rotation profiles show good agreement
between the semi-analytical method and simulations. Panel (b) reveals that the
semi-analytical method matches $\Sigma$ at the inner gap edge well, but over-
predicts it at the outer gap edge. This is due to a consistently stronger
perturbation in $\zeta^{\mathrm{SA}}$ in this region. Both methods show a
stronger asymmetry in the $\Sigma$ peaks at the inner and outer gap edges as
compared to the fiducial case, with the inner peak dominating. The radial
width of the gap is similar to the fiducial case, as expected, since
$l_{\mathrm{sh}}$ is unchanged.

Overall, the level of accuracy of our semi-analytical reconstruction for $p=0$
disc is similar to that in Fig. 3. We note, however, that for the lowest
planet mass case, $M_{\mathrm{p}}=0.1M_{\mathrm{th}}$, with
$h_{\mathrm{p}}=0.1$ and $p=0$, we find poorer agreement between simulations
and semi-analytical method in the inner disc. We attribute this to the fact
that in our semi-analytical model, we neglect the advection term in the
conservation equation for vortensity, as we found it to be negligible in Paper
I. However, for this parameter set, the shocking distance $l_{\mathrm{sh}}$ is
large, such that gap opening in the inner disc occurs in a region with
initially large vortensity gradient. As a result, the vortensity advection
term ($v_{R}\partial_{R}\zeta$) might become important as disc fluid gets
radially redistributed. This argument is supported by the fact that we do not
see such disagreement in the $p=3/2$ disc with radially constant
$\zeta_{\mathrm{i}}$ (with other parameters kept the same).

###  C.2 Variation of the disc scale-height $h_{\mathrm{p}}$

We next change the disc scale-height at the planet location $h_{\mathrm{p}}$
(and $c_{\mathrm{s}}$), considering hotter discs. In column B of Fig. 14 we
compare our gap reconstruction with simulations at $t=300P_{\mathrm{p}}$ for
fiducial values of the planet mass and surface density slope
($M_{\mathrm{p}}=0.25M_{\mathrm{th}}$ $p=1.5$) but an increased
$h_{\mathrm{p}}=0.1$. Again, we find good agreement between our semi-
analytical method and simulations.

In line with our findings in Paper I (e.g. Fig. 15 therein), the radial scale
(width) of vortensity and density perturbation induced by the spiral shocks
increases roughly linearly with $h_{\mathrm{p}}$ since $l_{\mathrm{sh}}\propto
h_{\mathrm{p}}$ for fixed $M_{\mathrm{p}}/M_{\mathrm{th}}$. While the
vortensity jump at the shock $\Delta\zeta$ remains almost unchanged, the rate
of vortensity production scales with the relative (synodic) orbital period of
gas parcels with respect to the shock, see equation (20). This means that
vortensity perturbations increase more rapidly in hotter discs (as
$l_{\mathrm{sh}}$ is larger), leading to higher surface density perturbations
at the same time after introduction of the planet: for $h_{\mathrm{p}}=0.1$,
$\delta\Sigma/\Sigma_{0}\sim 0.1$ at the gap edges is reached after about half
the time it takes in the $h_{\mathrm{p}}=0.05$ disc.

###  C.3 Variation of the planet mass $M_{\mathrm{p}}/M_{\mathrm{th}}$

The amplitude of the vortensity jump at the shock strongly depends on the
normalized planet mass (close to
$\Delta\zeta\propto\left(M_{\mathrm{p}}/M_{\mathrm{th}}\right)^{3}$, Paper I,
). Also, according to the equation (6), the width of the gap (which is
$\propto l_{\mathrm{sh}}$) decreases as $M_{\mathrm{p}}/M_{\mathrm{th}}$
increases. Both these scalings conspire to produce steeper gradients in
surface density and pressure, and thus increased shear, as
$M_{\mathrm{p}}/M_{\mathrm{th}}$ increases.

In columns C and D of Fig. 14 we show results for the fiducial disc parameters
with a lower mass planet $M_{\mathrm{p}}=0.1M_{\mathrm{th}}$ at
$t=6800P_{\mathrm{p}}$ and a higher mass planet
$M_{\mathrm{p}}=0.5M_{\mathrm{th}}$ at $t=60P_{\mathrm{p}}$. Comparing the two
illustrates the aforementioned trends. Note that the significance of the
secondary gap in the inner disc (which causes vortensity to deviate from our
reconstructed $\zeta$-profile, see a peak in $\zeta$ derived from simulations
at $R=0.6$ in Fig. 14i) relative to the primary gap is decreasing as
$M_{\mathrm{p}}/M_{\mathrm{th}}$ increases. This is due to the fact that the
density wave produced by a higher mass planet is more non-linear and loses its
angular momentum to the background flow more efficiently. This leads to a less
efficient formation of the secondary spiral arm (e.g. Paper I, ).

Note that regardless of the disc and planetary parameter choices we always
find the second term in the equation (25), shown in panels (d), (h), (l), (p),
to greatly exceed the first term in that equation,
$\delta\Sigma/\Sigma_{\mathrm{i}}$ (panels (b), (f), (j), (n)). This
emphasizes once again the importance of properly accounting for the small
variations of $\Omega$ when computing $\zeta(R)$.

##  Appendix D RWI analysis for a $p=0$ disc

Figure 15: Same as the top panel Fig. 4 but for an inner gap edge mode in the
$p=0$ disc at $t=650P_{\mathrm{p}}$. The frequency of this mode is
$\omega=3.874+0.01534\mathrm{i}$. Figure 16: Growth rates $\gamma$ of unstable
modes, associated with the inner gap edge, found at different times for the
constant surface density disc $(p=0$) and
$M_{\mathrm{p}}/M_{\mathrm{th}}=0.25$. Compare with Figure 6.

In a constant surface density $(p=0)$ disc, the inner gap edge becomes
unstable first as the vortensity evolves faster there, see column A of Fig.
14. Figure 15 illustrates this by showing the radial profile of the dominant
$m=2$ eigenmode in $p=0$ disc at $t=650P_{\mathrm{p}}$. One can see that it
peaks at $R\approx 0.85R_{\mathrm{p}}$, i.e. at the inner gap edge, meaning
that this side of the gap turns RWI-unstable first.

Due to the presence of the secondary spiral shock in the inner disc, not
captured by our semi-analytical model, we might then also expect greater
deviations of the growth rates from simulations in $p=0$ disc. Fig. 16
confirms this expectation by showing growth rates from our linear stability
analysis for $m=1$ and $m=2$ modes. While the general behaviour of the growth
rates (their overall increase with time) is similar to Fig. 6, the time offset
between the rates derived using the semi-analytical reconstruction (yellow)
and azimuthally-averaged simulation data is larger in the $p=0$ disc, around
$100P_{\mathrm{p}}$.

##  Appendix E Vortensity threshold for RWI onset

Figure 17: Data for $\delta\zeta_{\mathrm{th}}^{\mathrm{lin}}$ from Fig. 12 rescaled via multiplication by $h_{\mathrm{p}}/h_{0}$ ($h_{0}=0.05$). One can see the points collapsing on two well-defined branches for different values of $p$. A fit (48) to these data with parameters from Table 2 is shown by black dotted (for $p=0$) and dashed (for $p=3/2$) curves. Table 2: Fitting parameters for the threshold vortensity perturbation at $t=t^{\mathrm{SA}}_{\mathrm{lin}}$ (see equation 48 and Fig. 17).  | $p$ | $a_{1}$ | $a_{2}$ | $a_{3}$
---|---|---|---|---
$\delta\zeta_{\mathrm{min}}$ | 3/2 | 0.035 | -1.3 | 0.8
| 0 | -0.05 | -36 | 2.2
$\delta\zeta_{\mathrm{max}}$ | 3/2 | -0.43 | 4.4 | 0.48
| 0 | 0.43 | 153 | 2.1

We found that multiplication by $h_{\mathrm{p}}$ causes
$\delta\zeta_{\mathrm{th}}^{\mathrm{lin}}$ data from Fig. 12 to fall onto two
very well-defined branches as a function of $l_{\mathrm{sh}}$, distinguished
by the corresponding $p$. This is shown in Fig. 17. There we also show a fit
of the form

$\displaystyle\delta\zeta_{\mathrm{th}}^{\mathrm{lin}}=\left[a_{1}+a_{2}\,\left(\frac{l_{\mathrm{sh}}}{R_{\mathrm{p}}}\right)^{a_{3}}\right]\left(\frac{h_{\mathrm{p}}}{h_{0}}\right)^{-1},$
(48)

where $h_{0}=0.05$ and $(a_{1},a_{2},a_{3})$ are the fit parameters obtained
via a least-square fit, done separately for the two values of $p$ and for the
minima and maxima of vortensity at the inner (outer) gap edge for $p=0$
($3/2$), giving us four sets of these coefficients. These parameter sets are
listed in Table 2.

The correlations found in Fig. 17 clearly merit further investigation,
especially for other $p$, to determine the dependence of $(a_{1},a_{2},a_{3})$
on $p$. For now, we refrain from attaching any significance to these trends
beyond the statements made in Section 6.4.