11institutetext:  Institut für Astronomie und Astrophysik, Universität
Tübingen, Auf der Morgenstelle 10, 72076 Tübingen, Germany 22institutetext:
Astronomy Unit, School of Physical and Chemical Sciences, Queen Mary
University of London, London E1 4NS, UK
22email<EMAIL_ADDRESS><EMAIL_ADDRESS>
# Hydrodynamic turbulence in disks with embedded planets

Alexandros Ziampras Hydrodynamic turbulence in disks with embedded
planetsHydrodynamic turbulence in disks with embedded planetsHydrodynamic
turbulence in disks with embedded planetsHydrodynamic turbulence in disks with
embedded planets Wilhelm Kley W. Kley is deceased and is included as a co-
author for his significant contribution and guidance over the course of this
project.Hydrodynamic turbulence in disks with embedded planetsHydrodynamic
turbulence in disks with embedded planets Richard P. Nelson Hydrodynamic
turbulence in disks with embedded planetsHydrodynamic turbulence in disks with
embedded planets

The vertical shear instability (VSI) is a source of hydrodynamic turbulence
that can drive vigorous vertical mixing and moderate levels of accretion in
protoplanetary disks, and it could be observable in the near future. With
high-resolution three-dimensional numerical hydrodynamics simulations, we
modeled the behavior of the VSI in protoplanetary disks with and without
embedded planets. We then measured its accretion and mixing capabilities by
comparing the full Reynolds stress, which includes the contribution of
nonaxisymmetric features, such as spiral arms and vortices, to the Reynolds
stress due to the azimuthally averaged velocity field, which can be attributed
to good approximation to the VSI. We verified that the VSI can contribute to
the accretion stress and showed that, depending on disk conditions, an
embedded planet can coexist with or suppress VSI turbulent stress.
Specifically, the presence of spiral shocks launched by a planet or planet-
generated vortices can interfere with the VSI near the planet’s vicinity, with
the instability recovering at large enough distances from the planet or
vortex. Our results suggest that observations of VSI signatures are unlikely
in disks that contain massive, nonaxisymmetric features.

###### Key Words.:

planet formation, accretion disks, numerical hydrodynamics

##  1 Introduction

Turbulence is thought to be a key mechanism that affects the evolution of
circumstellar disks, as well as the planet formation process. From
facilitating accretion of gas onto the central star through sustained outward
angular momentum transport (Lynden-Bell & Pringle 1974; Shakura & Sunyaev
1973; Balbus & Papaloizou 1999), to the vertical stirring of dust grains
observable in millimeter observations with the ALMA telescope (e.g., Andrews
et al. 2018; Dullemond et al. 2018; Villenave et al. 2022), turbulence plays a
significant role in setting the environment in which planets form and grow.

In the context of planet–disk interaction, understanding turbulence is
fundamentally important in several aspects. In the immediate vicinity of an
embedded planet, turbulence can regulate the accretion of gas and dust
particles onto the planet, influencing its growth (Picogna et al. 2018). In
the case of massive planets, which can strongly affect the disk background by
opening gaps and forming pressure bumps in the gas (e.g., Kley & Nelson 2012),
turbulence dictates the formation and decay of vortices at the planet’s gap
edge (e.g., Li et al. 2005; de Val-Borro et al. 2007; Rometsch et al. 2021),
as well as the radial diffusion of dust that accumulates on the resulting
pressure bumps (Dullemond et al. 2018) that can lead to observable rings in
continuum emission (Huang et al. 2018). As a result, understanding how
turbulence operates on both large and small scales is necessary in studying
planet–disk interaction.

A variety of (magneto-)hydrodynamical mechanisms have been proposed as a means
of generating turbulence, depending on the mass, radial, and vertical
structure, as well as the ionization fraction of the disk (for a comprehensive
review, see Lyra & Umurhan 2019). Of particular interest to our work is the
vertical shear instability (VSI, Nelson et al. 2013), which is expected to be
active in the “dead zone” of the disk, where the magnetorotational instability
(MRI, Balbus & Hawley 1991) is suppressed due to low gas ionization fractions
(Gammie 1996) or nonideal MHD effects (Cui & Bai 2020) and therefore it cannot
operate efficiently.

The VSI, originally conceived by Goldreich & Schubert (1967) and Fricke (1968)
in the context of rotating stars, and having been brought into focus in the
protoplanetary disk context by Nelson et al. (2013), has received extensive
attention in the last decade as it arises naturally in typical protoplanetary
disks. It requires the presence of a vertical shear (i.e., a height-dependent
rotation profile) and relatively short cooling timescales (Lin & Youdin 2015),
both of which are satisfied at distances of $\gtrsim$10 au for typical disk
parameters in models with realistic disk thermodynamics (Stoll & Kley 2014;
Flock et al. 2017) and up to a height of approximately 3 gas scale heights
when accounting for dust–gas coupling (Pfeil & Klahr 2021). Its ability to
produce accretion rates with a turbulent $\alpha\sim 10^{-4}$–$10^{-3}$
(Nelson et al. 2013; Stoll & Kley 2014; Manger & Klahr 2018; Flores-Rivera et
al. 2020) coupled with the strongly anisotropic vertical stirring it generates
(Stoll et al. 2017b) deem the VSI a mechanism compatible with observed levels
of turbulence in protoplanetary disks (e.g., Dullemond et al. 2018; Flaherty
et al. 2020).

Regarding an embedded planet’s influence on the VSI, Stoll et al. (2017a) have
shown that VSI-active disks behave very similarly to traditionally viscous
ones in terms of the gap structure sculpted by and the migration torques
acting on a planet of up to $100\,\text{M}_{\oplus}$. They have also shown
that, even in this high-mass case, the disk maintains an accretion turbulence
parameter $\alpha\sim 10^{-4}$. The vertical stirring efficiency of the VSI,
however, is not addressed quantitatively, apart from an apparent damping of
the gas vertical velocity in the outer disk for
$M_{\text{p}}=100\,\text{M}_{\oplus}$, attributed in part to vortex activity.
This possible quenching of the VSI has important implications for any hopes of
observing the kinematic signatures of the VSI in molecular line emission, as
its signature vertical motion is extremely faint and very likely not
observable even by modern ALMA standards in disks without an azimuthal
structure (Barraza-Alfaro et al. 2021).

Our goal is to study the interaction between an embedded planet or vortex and
the VSI and isolate the contribution of the latter, so as to quantify the
extent to which the instability can coexist alongside nonaxisymmetric features
such as planet-generated vortices and spirals. In doing so, we aim to gauge
whether an observation of the VSI in disks with an azimuthal structure could
be realized in the near future using CO kinematics with ALMA.

We describe our physical setup and numerical methods in Sect. 2. We present
our results in Sect. 3, and discuss their implications in Sect. 4. We finally
list our conclusions in Sect. 5.

##  2 Physics and numerics

In this section we first provide a brief rundown of the physical model and
parameters we used and motivate our choices. We then describe our numerical
setup and the methods used to analyze our results, highlighting the specifics
of the averaging procedure for our time- and space-dependent datasets.

###  2.1 Physics

We solve the inviscid, time-dependent three-dimensional (3D) hydrodynamics
equations (Tassoul 1978) for an ideal gas with density $\rho$, pressure $P$
and velocity field $\bm{u}$ orbiting around a star with mass
$M_{\star}=1\,\text{M}_{\odot}$

$\displaystyle\frac{\partial{\rho}}{\partial{t}}$
$\displaystyle+\bm{u}\cdot\nabla\rho=-\rho\nabla\cdot\bm{u}$ (1)
$\displaystyle\frac{\partial{\bm{u}}}{\partial{t}}$
$\displaystyle+(\bm{u}\cdot\nabla)\bm{u}=-\frac{1}{\rho}\nabla P-\nabla\Phi$
(2) $\displaystyle\frac{\partial{e}}{\partial{t}}$
$\displaystyle+\bm{u}\cdot\nabla e=-\gamma e\nabla\cdot\bm{u}+S,$ (3)

where $\Phi=-\text{G}M_{\star}/r$ is the gravitational potential of the star
at distance $r$, $e=P/(\gamma-1)$ is the thermal energy density of the gas,
and $S$ encapsulates any additional terms that can affect the gas temperature
$T$. The gas has an adiabatic index $\gamma=\nicefrac{{7}}{{5}}$ and a mean
molecular weight $\mu=2.35$ such that $e=\rho c_{\text{v}}T$, with
$c_{\text{v}}=\frac{\mathcal{R}}{\mu(\gamma-1)}$ being its heat capacity at
constant volume. The gravitational constant and gas constant are indicated
with G and $\mathcal{R}$, respectively. Finally, the isothermal sound speed of
the gas is given by $c_{\text{s}}=\sqrt{P/\rho}$ and relates to the adiabatic
sound speed through $c_{\text{s}}^{\text{ad}}=\sqrt{\gamma}c_{\text{s}}$.

We now adopt a cylindrical coordinate system $\\{R,\phi,z\\}$ where
$r=\sqrt{R^{2}+z^{2}}$. We can derive a hydrodynamic equilibrium state for an
axisymmetric ($\nicefrac{{\partial}}{{\partial\phi}}=0$), nonaccreting disk
($\bm{u}=u_{\phi}\hat{\phi}$) with a vertically isothermal temperature profile
($\nicefrac{{\partial T}}{{\partial z}}=0$). If we further assume that the gas
density at the midplane $\rho_{\text{mid}}$ as well as the temperature both
follow power law profiles in the radial direction $R$ with

$\rho_{\text{mid}}(R)=\rho_{z=0}=\rho_{0}\left(\frac{R}{R_{0}}\right)^{p},\qquad
T(R)=T_{0}\left(\frac{R}{R_{0}}\right)^{q},$ (4)

we can write the density and velocity profiles in equilibrium (see for example
Nelson et al. 2013)

$\displaystyle\rho^{\text{eq}}(R,z)$
$\displaystyle=\rho_{\text{mid}}(R)\,\exp\left[\frac{\text{G}M_{\star}}{c_{\text{s}}^{2}}\left(\frac{1}{r}-\frac{1}{R}\right)\right],$
(5) $\displaystyle u_{\phi}^{\text{eq}}(R,z)$
$\displaystyle=R\Omega_{\text{K}}\left[1+(p+q)\left(\frac{H}{R}\right)^{2}+q\left(1-\frac{R}{r}\right)\right]^{\nicefrac{{1}}{{2}}},$
(6)

where $\Omega_{\text{K}}(R)=\sqrt{\text{G}M_{\star}/R^{3}}$ is the Keplerian
angular velocity and $H(R)=c_{\text{s}}/\Omega_{\text{K}}$ is the pressure
scale height. Through the latter we can also define the aspect ratio $h=H/R$.

We choose to relax the temperature to its initial profile $T_{t=0}$ using a
parametrized $\beta$-cooling approach (e.g., Gammie 2001), where the cooling
timescale is given by $t_{\text{cool}}=\beta/\Omega_{\text{K}}$ and the
corresponding source term is

$S_{\text{relax}}=-\rho
c_{\text{v}}\frac{T-T_{\text{init}}}{t_{\text{cool}}}\Rightarrow\frac{\partial{T}}{\partial{t}}=-\frac{T-T_{\text{init}}}{\beta}\Omega_{\text{K}}.$
(7)

For simplicity, we adopt a constant $\beta=0.01$ in our models. This would
correspond to a reference radius $R_{0}\sim 50$–100 au for typical disk
parameters, placing our models in the radial range observable by ALMA (e.g.,
Andrews et al. 2018), but our models are effectively compatible with locally
isothermal models in the literature ($\beta\rightarrow 0$; e.g., Stoll et al.
2017b; Manger & Klahr 2018; Barraza-Alfaro et al. 2021). Our choice of $\beta$
also results in a disk where the cooling criterion of Lin & Youdin (2015) is
globally satisfied, and as a result the VSI operates at nearly full capacity
(Manger et al. 2021).

The disk is initialized in an equilibrium state according to Eqs. (4), (5) and
(6) with $p=-\nicefrac{{3}}{{2}}$, $q=-1$, and a constant $h=0.05$ for most
models, with some models using $h=0.1$ instead. This corresponds to a
nonflared disk with $T(R)=\\{120.7,12.6\\}$ K at $\\{5.2,50\\}$ au for a
solar-mass star and a surface density profile $\Sigma(R)\propto
R^{-\nicefrac{{1}}{{2}}}$. Our models are scale-free in terms of the reference
radius $R_{0}$ and surface density $\Sigma_{0}=\Sigma(R=0)$, since we do not
consider the back-reaction of gas onto the star. While our choice of
parameters is not directly comparable to the temperature conditions at the
ALMA range ($q\approx-\nicefrac{{1}}{{2}}$, $h\approx 0.1$, see for example
Flock et al. 2020), it allows us to construct models where the generated VSI
turbulence corresponds to a constant effective $\alpha$ parameter (Shakura &
Sunyaev 1973) on the order of a few $10^{-4}$ similar to Flock et al. (2020),
while also maintaining compatibility with previous studies (e.g., Nelson et
al. 2013; Stoll et al. 2017b; Manger et al. 2020). This becomes especially
important when addressing the observability of VSI signatures, for which we
compare our findings to those of Barraza-Alfaro et al. (2021).

###  2.2 Numerics

We use the numerical hydrodynamics code PLUTO 4.3 (Mignone et al. 2007) in a
spherical polar $\\{r,\theta,\phi\\}$ geometry. PLUTO employs a finite-volume
approach with a Riemann solver to integrate Eqs. (1)–(3) in time. Our models
cover the full azimuthal range $\phi\in[0,2\pi)$, a radial range
$r\in[0.4,2.5]\,R_{0}$ and a vertical extent $z\in[-4,4]\,H$, for a total of 8
scale heights from end to end vertically. After conducting a resolution study,
we opted for a grid with 16 cells per scale height in the vertical and radial
directions and 5 scale heights in the azimuthal direction, which translates to
$N_{r}\times N_{\theta}\times N_{\phi}$ = $600\times 128\times 600$ cells for
$h=0.05$. We found that this cell count adequately resolves both the growth
and saturation phases of the VSI, while also keeping simulation runtimes
reasonably feasible given their computational cost. A resolution analysis on
our simulated VSI growth rates is provided in Appendix B.

Whenever possible, we use the FARGO module introduced by Masset (2000) and
implemented into PLUTO by Mignone et al. (2012). This yields a necessary
speedup factor of $\sim$10 by subtracting the Keplerian background flow before
solving the Riemann problem across all cell interfaces in the domain. We
choose the HLLc Riemann solver (Harten 1983) with the Van Leer flux limiter
(Van Leer 1974) and second-order reconstruction and time advancement schemes
(options LINEAR and RK2 respectively), but note that more accurate schemes
such as third-order-accurate methods or less diffusive Riemann solvers and
flux limiters yielded functionally identical results compared to our fiducial
2D model (discussed in Sect. 3.1).

We apply wave-damping boundary conditions in the radial domain
$R\in[0.4,0.5]\cup[2.1,2.5]\,R_{0}$ following de Val-Borro et al. (2006),
using a damping timescale of ten orbits at the corresponding boundary edge.
This suppresses the reflection of planet-generated spiral wakes while also
providing a mass reservoir for the accreting, VSI-active disk. To further
prevent waves of any kind from reentering the domain, we utilize “strict
outflow” boundary conditions at the radial and polar directions, meaning that
we allow outflow but not inflow of material at the boundary edge. In the polar
direction, we extrapolate the density and pressure into the ghost cells
assuming vertical hydrostatic equilibrium.

In models with an embedded planet, the latter is treated as a point-like
gravitational potential with mass $M_{\text{p}}=3\times 10^{-4}\,M_{\star}$
($100\,\text{M}_{\oplus}$ for a solar-mass star), smoothed over a distance of
half its Hill sphere following the polynomial prescription of Klahr & Kley
(2006). We neglect gas back-reaction on both the planet and star. While in the
linear regime we only expect the planet to launch spiral arms (Ogilvie & Lubow
2002), planets that are massive enough can open one or more gaps in their
vicinity (Rafikov 2002). The transition from one regime to the other can be
defined using the disk thermal mass (Zhu et al. 2015)

$M_{\text{th}}=\left.\frac{c_{\text{s}}^{3}}{\text{G}\Omega_{\text{K}}}\right|_{R_{\text{p}}}\approx
1\,M_{\text{J}}\leavevmode\nobreak\
\left(\frac{h_{\text{p}}}{0.1}\right)^{3}\frac{M_{\star}}{\text{M}_{\odot}},$
(8)

which, for our choices of aspect ratio $h=\\{0.05,0.1\\}$ evaluates to
$\\{0.125,1\\}\,M_{\text{J}}$. As a result, with our planet mass choice we can
probe both ends of the planet–disk interaction spectrum, from the quasi-linear
($h=0.1$, $M_{\text{p}}=0.3\,M_{\text{th}}$) to the nonlinear regime
($h=0.05$, $M_{\text{p}}=2.4\,M_{\text{th}}$).

Figure 1: Timeline for our 3D models with embedded planets. Models in the
blue, orange and green phases are discussed in Sects. 3.1, 3.2 and 3.3
respectively. The planet finishes growing by $t=t_{0}+100\,P_{0}$ (purple
dotted line).

###  2.3 Model timeline

Taking advantage of the axisymmetric nature of the VSI we first construct 2D
models of the VSI on the $\\{r,\theta\\}$ plane. These models, while
axisymmetric, involve all 3 velocity components. We then carry out a
resolution analysis, in order to measure the turbulent viscosity $\alpha$ in a
way compatible with previous 2D studies and verify that relevant VSI features
in our 3D models are adequately resolved. While the VSI typically saturates
within 30 local orbits (Nelson et al. 2013) which would translate to a
required simulation time of $t_{\text{2D}}\approx 2.5^{3/2}\,P_{0}$ $\approx$
120 orbits at $R=R_{0}$, we continue until $t_{\text{2D}}=1000\,P_{0}$ to
ensure the disk reaches a quasi-steady state while also allowing a reasonably
long timeframe of 500 $P_{0}$ (between 500–1000 $P_{0}$, with an output
frequency of $1/P_{0}$) which we use to compute our time-averaged Reynolds
stress once this quasi-steady state is reached. The 2D models are initialized
at dynamical equilibrium through Eqs. (5) & (6), with all velocity components
perturbed by $1\%\,c_{\text{s}}$ to encourage the growth of VSI modes. We
discuss our 2D results in Sect. 3.1.

We then expand the disk into 3D by including the azimuthal direction. Similar
to Flock et al. (2020), we duplicate the last state of our 2D models along the
$\phi$-axis while also seeding the velocity with a further $1\%\,c_{\text{s}}$
perturbation to break the axisymmetry and allow the formation of azimuthal
structure. The 3D models are then evolved for a further $200\,P_{0}$, and
related results are presented in Sect. 3.2.

Finally, a planet is inserted at $R=R_{0}$ in each of the now quasi-steady 3D
disks and grows to its final mass over 100 orbits using the formula in de Val-
Borro et al. (2006). The planet–disk models are then evolved for another
$400\leavevmode\nobreak\ P_{0}$, until all relevant planet-generated features
(spirals, gaps, vortices) have fully developed and the planet–VSI interplay
can be assessed. Results from this phase of the timeline are discussed in
Sect. 3.3.

A schematic of our model timeline is shown in Fig. 1. The 2D, 3D and
planet–disk phases are color-coded in blue, orange and green, respectively.

##  3 Results

In this section we present the results of our simulations. We begin by
profiling the VSI, using an axisymmetric model to measure the growth rate and
resulting stress levels in a saturated state. We then highlight the impact of
3D effects such as spiral arms and vortices on VSI activity. Finally, we
analyze the effect of planet-generated features (spirals, vortices, rings and
gaps) on VSI activity and attempt to disentangle the individual contributions
of the planet and the VSI to the levels of turbulence we observe in our
simulations.

###  3.1 Reference model

As has been shown before by numerous studies (e.g., Nelson et al. 2013; Stoll
& Kley 2014; Flock et al. 2017), the VSI is characterized by quasi-
axisymmetric, vertically elongated flows of gas that span the vertical extent
of the disk while covering a proportionally much narrower distance in the
radial direction. The result is the formation of channels of vertical gas
motion between the two disk surfaces that turn over near each surface, before
heading back toward the midplane. This is easily discernible through the
vertical component of the gas velocity field $u_{z}$, which shows sheets of
up- and downward motion along the $z$ direction throughout the disk (Fig. 2).

Figure 2: Radial ($u_{R}$), vertical ($u_{z}$), and azimuthal ($u_{\phi}$) gas
velocity components in our fiducial axisymmetric model, highlighting the
vertically elongated sheets of motion generated by the VSI. The Keplerian
background $u_{\text{K}}$ has been subtracted from $u_{\phi}$, and all
components have been normalized to the local sound speed.

Consistent with previous studies, we compute a VSI-driven accretion stress
$\alpha_{\text{acc}}\approx 1\text{--}3\times 10^{-4}$ and a vertical
turbulent mixing parameter $\alpha_{\text{mix}}\approx 0.1\text{--}0.3$ for
our choice of disk parameters. To measure $\alpha$ we adopt the method in
Balbus & Papaloizou (1999), which is similar to that of Stoll et al. (2017b)
albeit also accounting for the background velocity field (rather than assuming
a stationary state with $\bar{u}_{R}=0$) and further adjusting the result by a
factor of $\nicefrac{{2}}{{3}}$. Our method is described in detail in Appendix
A, and our results for $\alpha_{\text{acc}}$ and $\alpha_{\text{mix}}$ as
functions of $R$ and $z$ respectively are shown in Fig. 3.

Figure 3: Radial and vertical profiles of the turbulent accretion parameter
$\alpha_{\text{acc}}$ and the turbulent mixing parameter $\alpha_{\text{mix}}$
in our fiducial axisymmetric model. We find a typical
$\alpha_{\text{acc}}\approx 1.5\times 10^{-4}$ and $\alpha_{\text{mix}}\approx
0.15\text{--}0.3$, in agreement with previous results. The faint curve on the
upper panel refers to $\alpha_{\text{acc}}$ without radial smoothing through
Eq. (14).

###  3.2 3D effects

While the VSI is primarily axisymmetric, including the azimuthal direction in
fully 3D models can allow the growth of small-scale vortices at sufficient
resolution (Richard et al. 2016; Manger & Klahr 2018) or trigger the Rossby-
wave instability (RWI) at the radial location where the disk transitions from
quiet to VSI-turbulent due to the formation of a gap edge (Flock et al. 2020).
Such vortices launch spiral arms as they decay (Rometsch et al. 2021),
generating a Reynolds stress that can contribute to the total accretion stress
in the disk (Goodman & Rafikov 2001).

Figure 4: Highlights of our 3D fiducial model. Left: surface density
deviations in the disk compared to the azimuthally averaged surface density
$\bar{\Sigma}$, showing azimuthal structure in the form of spirals and
vortices. Middle: vertical velocity of the gas at the disk midplane; quasi-
axisymmetric stripes reveal VSI activity. Right: A map of the gas vortensity
deviations compared to the initial (Keplerian) vortensity $\varpi_{0}$,
highlighting two vortices orbiting at $R\approx 1.15\,R_{0}$.

Figure 4 shows a snapshot of the gas surface density (left) and vertical
velocity component (middle) for our fiducial 3D model that started from an
axisymmetric, fully VSI-saturated state and has evolved for 200 orbits at
$R=R_{0}$. The gas shows azimuthal structure, with spiral arms permeating the
disk and two vortices orbiting at $R\approx 1.15\,R_{0}$. Nevertheless, the
disk is still VSI-active, as indicated by the largely axisymmetric velocity
component $u_{z}$ with $u_{z}^{\text{mid}}\approx 0.2\,c_{\text{s}}$, which is
comparable to what we found in our axisymmetric VSI-saturated models in Fig.
2. To further highlight the presence of the vortices, we also plot the
vortensity (Masset & Benítez-Llambay 2016)

$\varpi=\left[\int\limits_{-\infty}^{+\infty}\frac{\rho}{(\nabla\times\bm{u})\cdot\hat{z}}\text{d}z\right]^{-1},$
(9)

on the right panel of the same figure.

In our 3D simulations we also observed the development of a parasitic
instability in the inner boundary, likely caused by the damping of VSI modes
on the interface between the active domain and the wave damping zone. This
instability induced strong vertical motion near the interface of the inner
damping zone ($R\approx 0.4\,R_{0}$) and contributed to the formation of
several small-scale vortices near the inner radial boundary. Similar to the
vortices shown in Fig. 4, they launched spirals that spanned the entire
domain, similar to what was observed by Barraza-Alfaro et al. (2021), likely
contributing to the total Reynolds stress. We note that, in the planet–disk
phase of our simulations, the planet overshadows the activity of such
vortices. We provide more details on this phenomenon in Appendix C.

Figure 5: Turbulent parameters $\alpha_{\text{acc}}$ and $\alpha_{\text{mix}}$
in our fiducial 3D model. Our 2D results from Fig. 3 are shown in blue curves.
The total stress in the disk is shown with solid curves, while dashed curves
isolate the contribution of the VSI. While $\alpha_{\text{acc}}^{\text{tot}}$
increases slightly, $\alpha_{\text{mix}}$ decreases such that
$\alpha_{\text{mix}}^{\text{tot}}\approx\alpha_{\text{mix}}^{\text{VSI}}$.
This happens due to the stronger accretion stress supplied by spiral arms that
in fact weaken the VSI, which is the primary driver of vertical motion. The
bottom panel is calculated at the radial location shown by the vertical line
in the top panel (here, $R=1.4\,R_{0}$).

We now measure turbulent stress by repeating the procedure in Sect. 3.1, this
time computing $\alpha$ using the full 3D gas data. As highlighted above and
shown by the solid orange lines in Fig. 5, this method captures the accretion
stress generated by the spiral arms seen in our 3D models and results in
slightly higher $\alpha_{\text{acc}}$ compared to our 2D models. However, the
disk shows a weaker vertical mixing parameter $\alpha_{\text{mix}}$ in the
region $R\sim R_{0}$, which is likely more strongly affected by the vortex
present there, as well as in the region $R\lesssim 0.7\,R_{0}$, which is
dominated by the spirals launched by the vortices that formed near the inner
boundary.

By assuming that, approximately, the total stress is
$\alpha^{\text{tot}}=\alpha^{\text{VSI}}+\alpha^{\text{spiral}}$, we attempt
to isolate the contribution of the VSI to the total turbulent stress by
recomputing $\alpha$ after first averaging our snapshots along the $\phi$ axis
as explained in Appendix A. This results in spiral-free profiles, and the
computed stress (dashed orange curves in Fig. 5) can be directly compared to
the 2D models in Sect. 3.1 (blue curves) and contrasted against the fully 3D
$\alpha$ stress (solid orange curves). We find that regions with a low
$\alpha_{\text{mix}}^{\text{tot}}$ coincide with those with low
$\alpha_{\text{mix}}^{\text{VSI}}$ and $\alpha_{\text{acc}}^{\text{VSI}}$,
even though $\alpha_{\text{acc}}^{\text{tot}}$ stays consistently high due to
spiral-driven Reynolds stress. In other words, while the spiral-rich 3D disk
shows a slightly higher accretion stress (as has been shown by previous 3D
studies such as Manger & Klahr 2018), the contribution of the VSI to this
total stress is in fact lower due to the spiral arms interfering with VSI
modes, and vertical mixing in this 3D disk is lower than what an axisymmetric
model would suggest. The latter supports the idea that the main driver of
vertical mixing ($\delta u_{z}$) is the VSI, while the spiral arms launched by
these relatively weak vortices mainly affect the radial and azimuthal
components of gas velocity ($\delta u_{R}$, $\delta u_{\phi}$) instead,
driving accretion.

Our argument, however, also implies that the spiral arms in our 3D runs should
be responsible for part of the accretion stress observed. Larson (1990) showed
that spiral shocks (such as those present in our models) can produce an
effective $\alpha_{\text{acc}}^{\text{spiral}}\approx 0.013h\sqrt{h+0.08}$, or
$\approx 2\times 10^{-4}$ for $h=0.05$. This estimate is compatible with our
results for $\alpha_{\text{acc}}^{\text{tot}}$, and supports the idea that the
observed accretion stress is driven by a combination of spiral arms and VSI
modes in the disk.

###  3.3 Planet–VSI interaction

In the previous section we showed how spiral arms and vortex activity can
inhibit the VSI by interfering with its vertical mixing capabilities while
maintaining a high accretion stress. Here, we extend our models by embedding
planets in our disks. While Stoll et al. (2017a) have already examined the
planet–VSI interplay as a function of planet mass, qualitatively showing that
VSI activity is weaker near the planet for sufficiently massive planets, we
intend to disentangle and quantify the contribution of the VSI to the total
accretion and mixing Reynolds components compared to that of the planet and
the features it generates (spirals, gaps, vortices).

Here, we focus on two models representative of the sub- and superthermal
planet mass regimes by first showing the qualitative effect of the planet on
the vertical velocity component of the disk, before calculating both
$\alpha^{\text{tot}}$ and $\alpha^{\text{VSI}}$ in order to highlight the
effect of the planet on the VSI. We then collect our estimates on
$\alpha_{\text{acc}}$ and $\alpha_{\text{mix}}$ with and without the
contribution of nonaxisymmetric features for all 3D runs in Table 1.

Table 1: Estimates of turbulent $\alpha$ in our models. Planet mass is quoted in $M_{\text{th}}$. Values for models with embedded planets serve as upper limits, especially for the last model where the VSI is only present in in the inner disk (see Fig. 9). $h$ | $M_{\text{p}}$ | $\alpha_{\text{acc}}^{\text{tot}}$ | $\alpha_{\text{mix}}^{\text{tot}}$ | $\alpha_{\text{acc}}^{\text{VSI}}$ | $\alpha_{\text{mix}}^{\text{VSI}}$
---|---|---|---|---|---
0.10 | — | $4.51\times 10^{-4}$ | $0.115$ | $2.53\times 10^{-4}$ | $0.103$
0.10 | 0.3 | $8.39\times 10^{-4}$ | $0.058$ | $1.40\times 10^{-4}$ | $0.045$
0.05 | — | $1.82\times 10^{-4}$ | $0.128$ | $1.06\times 10^{-4}$ | $0.117$
0.05 | 2.4 | $1.32\times 10^{-3}$ | $0.056$ | $4.21\times 10^{-5}$ | $0.039$

####  3.3.1 Quasi-linear regime: $h=0.1\rightarrow
M_{\text{p}}=0.3\,M_{\text{th}}$

Figure 6 shows an example of the typical disk structure around a subthermal-
mass planet with $M_{\text{p}}=0.3\,M_{\text{th}}$, taken after 100 planetary
orbits in our model with $h=0.1$, $M_{\text{p}}=100\,\text{M}_{\oplus}$. While
the planet’s spiral arms are visible (top panel), the disk is otherwise weakly
affected by the planet and the VSI remains operational, seemingly undeterred
at least qualitatively (bottom panel).

Similar to our analysis in Sect. 3.2, we now calculate the accretion stress
$\alpha_{\text{acc}}$ and vertical mixing parameter $\alpha_{\text{mix}}$ by
first considering the full 3D structure of the disk (which includes the
effects of spiral arms and results in a total stress $\alpha^{\text{tot}}$),
and then comparing our findings to an estimate of the VSI-driven stress
$\alpha^{\text{VSI}}$ by first averaging all quantities along the $\phi$ axis.

Our results, shown in Fig. 7, show that the presence of the planet does indeed
affect both the total Reynolds stress in the disk, as well as the isolated
contribution of the VSI to said stress. The total accretion stress
$\alpha_{\text{acc}}^{\text{tot}}$ (top panel, solid orange curve) increases
when compared to the 3D run without an embedded planet (solid blue curve, see
also the comparison between 2D and 3D for $h=0.05$ in Fig. 5) due to the
spiral arms launched by said planet, in line with our expectations from Sect.
3.2. At the same time, spiral activity interferes with the vertical velocity
field, weakening VSI-driven stress $\alpha_{\text{acc}}^{\text{VSI}}$ (dashed
orange curve) by a factor of $\approx 2$ compared to its contribution in the
3D run without a planet (dashed blue curve). The same effect can be observed
for the vertical mixing efficiency through $\alpha_{\text{mix}}$ on the right
panel of Fig. 7, as the VSI is the main driver of vertical motion (and thus
$\alpha_{\text{mix}}$).

Figure 6: Surface density deviations and vertical velocity maps at the
midplane for our disk with an embedded subthermal-mass planet ($h=0.1$,
$M_{\text{p}}=0.3\,M_{\text{th}}$) at $R=R_{0}$ after 500 planetary orbits.
While the planet’s primary spiral is prominent in this color range, the VSI
continues to operate as shown by the quasi-axisymmetric stripes of vertical
motion in the bottom panel.

While the embedded planet can also excite vertical motion via buoyancy waves
(Zhu et al. 2015) we could not find decisive evidence of their existence, most
likely due to the quite low $\beta=0.01$. An additional contribution to
$\alpha_{\text{acc}}$, also suggested by Stoll et al. (2017a), would very
likely exist due to the spiral wave instability (SWI, Bae et al. 2016a, b)
producing an accretion stress with an equivalent $\alpha_{\text{acc}}\sim
10^{-4}$ for this low planet mass.

All in all, in the quasilinear, subthermal-mass regime, we conclude that the
VSI can coexist with a planet, albeit in a slightly weaker state. A model with
$h=0.05$ and $M_{\text{p}}=13\,\text{M}_{\oplus}$, yielding
$0.3\,M_{\text{th}}$, resulted in very similar behavior, further highlighting
that $M_{\text{th}}$ is the relevant quantity to gauge the planet’s effect on
the disk regardless of disk model. In the next section, we investigate a
scenario with a superthermal-mass planet.

Figure 7: Turbulent parameters $\alpha_{\text{acc}}$ and $\alpha_{\text{mix}}$
presented similar to Fig. 5, where the blue lines instead refer to the
corresponding 3D model before the planet is introduced. As with Fig. 5,
$\alpha_{\text{acc}}$ increases but the overall mixing driven by the VSI is
weaker by a factor of 2.

####  3.3.2 Nonlinear regime: $h=0.05\rightarrow
M_{\text{p}}=2.4\,M_{\text{th}}$

Contrary to the disk model with $h=0.1$ discussed in the previous paragraph, a
planet with $M_{\text{p}}=100\,\text{M}_{\oplus}$ amounts to
$2.4\,M_{\text{th}}$ in a disk with $h=0.05$. This superthermal-mass planet,
shown in Fig. 8, shows significantly more features and affects the disk more
strongly and on a larger radial scale. Here, the disk shows much more
pronounced spirals as well as two gaps, with the center of the secondary gap
agreeing very well with the position predicted by Zhang et al. (2018) for a
low-viscosity protoplanetary disk ($\alpha_{\text{acc}}\lesssim 5\times
10^{-5}$). The conditions for the RWI are satisfied at the outer edge of both
gaps, giving rise to two large-scale vortices. In particular, the massive
vortex orbiting the planet’s primary gap encloses a gas mass of
$3.3\,M_{\text{p}}=1\,\mathrm{M}_{\mathrm{J}}$ and survives over the entire
duration of the simulation, continuously launching spiral arms that are
comparable in magnitude to those formed by the planet into the outer disk.

In this regime, we see both qualitatively (Fig. 8) and quantitatively (Fig. 9)
that the VSI is significantly weaker, if not completely quenched, in the outer
disk ($R\gtrsim R_{0}$, see in particular Fig. 8–middle or Fig. 9–top). This
is a result of the combined activity of the planet’s strong spiral shocks near
the planet, the massive vortex at the outer edge of the primary gap, and the
presence of both planet- and vortex-driven spiral arms in the entire outer
disk.

Figure 8: Two-dimensional heatmaps of relevant quantities similar to Fig. 4,
showcasing the planet’s effect on the disk after 500 orbits. The planet carves
a deep gap and forms an additional, secondary ring (top) which lead to the
formation of two massive vortices (top, bottom). VSI activity is significantly
weaker in the outer disk but present in the inner disk (middle). The initial
surface density and vortensity profiles are labeled with $\Sigma_{0}$ and
$\varpi_{0}$. Figure 9: Turbulence parameters $\alpha_{\text{acc}}$ and
$\alpha_{\text{mix}}$ similar to Fig. 7, but for the superthermal-mass case
with $h=0.05$, $M_{\text{p}}=2.4\,M_{\text{th}}$. The surface density snapshot
of Fig. 8 is also shown on panel _a_ , to highlight features present at
different radii. The parameter $\alpha_{\text{acc}}$ (panel _b_) is dominated
by the planet’s spiral arms, but the VSI maintains a contribution in the inner
disk (see also panel _c_). In the outer disk, VSI activity is weaker by 1–3
orders of magnitude compared to its fiducial value (see also panel _d_). A
gray band marks the planet’s gap region, defined as $R_{\text{p}}\pm
2.5\,R_{\text{Hill}}$.

We further show that the VSI remains active in the inner disk ($R\lesssim
R_{0}$), consistent with Stoll et al. (2017a), even though a vortex is present
at the edge of the planet’s secondary gap (see Fig. 9). This likely relates to
the stronger vertical shear in the inner disk, which is proportional to the
local $\Omega_{\text{K}}$, as well as the faster growth rate of the
instability in the inner disk ($\propto\Omega_{\text{K}}^{-1}$) when compared
to the timescale over which spiral arms, which corotate with the planet
($\propto\Omega_{\text{p}}^{-1}$), propagate radially through the gas. The
quicker growth of the VSI in the inner disk, combined with the smaller, weaker
vortex there compared to the massive vortex in the outer disk, allows the VSI
to recover to levels comparable to those found in 2D or 3D runs without an
embedded planet (see top and middle panels of Fig. 9).

Overall, a superthermal-mass planet in a nearly inviscid disk generates strong
nonaxisymmetric features (spiral shocks, massive vortices and subsequent
spirals) that can directly quench the VSI in their vicinity and especially so
in the outer disk ($R\gtrsim R_{\text{p}}$). We therefore conclude that the
VSI is very unlikely to coexist with massive planets ($M_{\text{p}}\gtrsim
M_{\text{th}}$), consistent with the qualitative findings of Stoll et al.
(2017a) and the synthetic observations by Barraza-Alfaro et al. (2022).

###  3.4 Vortex activity

To further elucidate the effect of vortices and their spiral arms that was
discussed in the previous section, we construct an additional 3D model based
on our fiducial setup ($h=0.05$), with the goal being to artificially spawn a
massive vortex into the disk and measure the turbulent stress in its vicinity
similar to Sects. 3.2 and 3.3. We initialize the disk with a circular density
gap centered at $R=R_{0}$, and place a blob of gas with total mass
$1\,\mathrm{M}_{\mathrm{J}}$ at the outer gap edge using a Gaussian
distribution centered at $R=1.25\,R_{0}$, $\phi=\pi$ (see top panel of Fig.
10)

$\rho(R,\phi)=\rho^{\text{eq}}\,\times\,f_{\text{gap}}(R)\,\times\,f_{\text{vortex}}(R,\phi),$
(10)

with $f_{\text{gap}}$ and $f_{\text{vortex}}$ motivated by our embedded-planet
models in Sect. 3.3:

$\begin{split}f_{\text{gap}}(R)=1-0.95&\exp\left[-\frac{(R/R_{0}-1)^{2}}{0.05}\right]\\\
f_{\text{vortex}}(R,\phi)=1+3.5&\exp\left[-0.5\left(\frac{R/R_{0}-1.25}{0.1}\right)^{2}\right]\\\
\times&\exp\left[-0.5\left(\frac{\phi-\pi}{0.35}\right)^{2}\right].\\\
\end{split}$ (11)

Figure 10: Surface density, vertical velocity and vortensity snapshots for our
run with an embedded vortex after 400 orbits at $R_{0}$. To eliminate the
stripes of vortensity deviations (see Fig. 4) and highlight the $\varpi$
minimum about the vortex, we normalize $\varpi$ to its azimuthal median
$\langle{\varpi}\rangle$.

The gas mass is initialized with Keplerian vorticity and quickly shears into
small gas clumps spread along the gap edge due to the radially-varying
$\Omega(R)$ profile. At the same time the RWI, the conditions for which are
satisfied at the gap edge, causes the formation of small-scale vortices at the
location of each clump, which over time coalesce into a single vortex. The
result is the formation of a massive vortex that orbits at the outer gap edge,
launching spiral arms through the otherwise VSI-active disk.

Figure 11: Turbulent $\alpha_{\text{acc}}$ and $\alpha_{\text{mix}}$ in our
vortex runs, similar to Fig. 5. VSI activity is significantly lower near the
radial vicinity of the vortex ($R\approx 1.3\,R_{0}$).

Fig. 11 shows the results of our data reduction for these models. Consistent
with our findings in Sects. 3.2 and 3.3, we find that $\alpha_{\text{acc}}$ is
slightly stronger in the model with a vortex orbiting the gap edge due to the
spirals generated by the vortex. At the same time, $\alpha_{\text{mix}}$ drops
near the vortex, signaling weaker VSI activity due to the disturbance of the
gas velocity field caused by the vortex.We also highlight that the drop in gas
density inside the gap does not seem to affect $\alpha$ in any noticeable way.

##  4 Discussion

In this section we justify our assumptions regarding our cooling prescription
and discuss our findings with regard to observations of the VSI. We also
briefly relate our findings to previous studies.

###  4.1 Radiative effects

Treatment of radiative effects—be it through a local, opacity-dependent
parametrized cooling prescription (e.g., Pfeil & Klahr 2021) or a full
treatment of radiation transport via flux-limited diffusion (FLD, Stoll & Kley
2014; Flock et al. 2020)—will certainly influence the planet’s gap-opening
capabilities (Ziampras et al. 2020b; Miranda & Rafikov 2020). Nevertheless,
the position of spirals will remain largely unaffected and the resulting
changes in the overall temperature profile will be minuscule in the regime
where the VSI is expected be active ($R\gtrsim 10$ au, Flock et al. 2020) due
to the significantly weaker spiral shock heating at these distances from the
star (Rafikov 2016; Ziampras et al. 2020a). For this reason, and since our
topic of interest is the planet–VSI interaction in a fully VSI-active disk
where radiative effects have been shown to not affect VSI activity
significantly (Stoll & Kley 2014; Manger et al. 2021), we do not carry out FLD
models in 3D.

###  4.2 On the observability of the VSI

Barraza-Alfaro et al. (2022) showed that the observational signature of the
VSI in a disk with an embedded massive planet is significantly weaker by
simulating ALMA gas kinematics observations. Stoll et al. (2017a) similarly
showed that, using $u_{\text{z}}$ as a tracer for VSI activity, the
instability should be weaker in the outer disk for
$M_{\text{p}}=100\,\text{M}_{\oplus}$, or $2.4\,M_{\text{th}}$ in their
models, but did not measure the vertical mixing parameter
$\alpha_{\text{mix}}$ in their simulations and instead focused on accretion
stress ($\alpha_{\text{acc}}$).

Our models, for which we provide quantitative values of $\alpha_{\text{mix}}$,
support the findings of both studies and enrich them further by providing
insight into the effect of a large-scale vortex on $\alpha_{\text{mix}}$. This
agreement also supports the robustness of our method of isolating the
contribution of the VSI to the total stress by simply azimuthally averaging
our data before computing the turbulent stress. It is also worth noting that
our vertical profiles of $\alpha_{\text{mix}}(z)$ are approximately constant,
which means that the VSI does not recover at large heights (see for example
Fig. 11) even though features such as vortices are significantly denser near
the midplane, with no sign of their activity near the disk surfaces.

This raises, however, the question whether VSI signatures could be observable
in protoplanetary disks with azimuthal structures such as spiral arms,
vortices, or planets. The findings of Barraza-Alfaro et al. (2021) suggest
that the vertical motion induced by the VSI could be observable in the near
future in disks where the VSI can operate unhindered and determine the
velocity field of the gas. By combining our results with those of Stoll et al.
(2017a) and Barraza-Alfaro et al. (2022), we conclude that it might be
possible, albeit challenging, to observe VSI activity in disks with
subthermal-mass planets or vortices and for typical disk parameters for ALMA
($h\sim 0.1$, $q\approx-0.5$, such that $\alpha^{\text{VSI}}\sim 10^{-4}$) in
the future, should the instability be operating in the disk. However, we also
deduce that it is unlikely that VSI signatures could be detected in the
exterior regions of a disk containing a planet whose mass exceeds the thermal
mass, since the planet-induced strong spiral shocks quench the activity of the
instability.

###  4.3 On the lifetime of vortices in our models

In our model with $M_{\text{p}}=2.4\,M_{\text{th}}$ a planet opened two gaps,
the outer edge of which grows RWI-unstable and spawns large-scale vortices. In
particular, the vortices that form outside of the planet’s primary gap grow to
masses of up to $3.3\,M_{\text{p}}=1\,\mathrm{M}_{\mathrm{J}}$. Given the
nearly-locally isothermal disk conditions ($\beta=0.01$) and the lack of VSI
activity in their vicinity resulting in a nearly inviscid disk with
$\alpha\sim\max(\alpha_{\text{num}},10^{-7})\ll 10^{-4}$, such vortices are
expected to survive for tens of thousands of planetary orbits (Rometsch et al.
2021).

However, the vortex growing on the outer edge of the secondary gap is smaller
in comparison. This vortex, also visible in Fig. 8, is not massive enough to
strongly affect the vertical velocity field in its vicinity (bottom panel of
Fig. 8), and as a result it is effectively embedded in a VSI-active region
with $\alpha\approx 5\times 10^{-5}$ (see top panel of Fig. 9), where it
decays over $\approx 10^{3}$ local orbits at $R\approx 0.7\,R_{0}$. This
lifetime is consistent with the estimates of Rometsch et al. (2021) for a disk
with $\alpha=10^{-4}$, $\beta=0.01$ and $h=0.05$, suggesting that the
turbulence driven by the VSI is compatible with 2D laminar-$\alpha$ models.
This last point complements the findings of Stoll et al. (2017a), who showed
that $\alpha$\- and VSI-driven disks behave very similarly in terms of planet
gap opening and torques exerted on embedded planets.

###  4.4 On resolution and grid requirements for planet–VSI modeling

Recent work by Flores-Rivera et al. (2020) has shown that the VSI is subject
to the Kelvin–Helmholtz instability on the interface between its vertically-
shearing corrugation modes, and Cui & Latter (2022) have further explored the
nonlinear saturation of the VSI via a parametric instability. Both studies
highlight the necessity of high-resolution modeling ($\gtrsim 200$–300 cells
per scale height) in order to fully resolve the interplay between all
mechanisms present in VSI-active disks. Since we are mainly interested in the
stress generated by the VSI and the planet–VSI interaction, which can be
resolved adequately with our grid resolution, and for computational reasons,
we chose not to run models at a higher resolution.

In addition, Svanberg et al. (2022) have shown that the VSI saturates into
coherent, large-scale, outwardly-traveling inertial waves. In this framework,
the planet can be seen as a potential that disrupts this outward propagation
of VSI waves beyond its semimajor axis. While our data supports this idea,
with VSI stress plummeting in the outer disk in our models with
$M_{\text{p}}=M_{\text{th}}$ (Sect. 3.3.2), it is difficult to explore this
scenario without actually carrying out simulations with a large enough radial
range that can accommodate the several radial zones that the VSI partitions
the disk into.

Based on the above, a much higher grid resolution as well as larger radial
domain would be necessary to resolve the global structure of the VSI and the
instabilities that can attack it. The interplay between the VSI and the
mechanisms covered in the above studies in disks with embedded planets is
worth exploring in the future.

##  5 Conclusions

We executed 3D numerical hydrodynamics simulations to study the vertical shear
instability (VSI) and its interaction with nonaxisymmetric features such as
embedded planets or vortices in a protoplanetary disk. We then isolated the
contribution of the VSI to the total radial accretion and vertical mixing
turbulence parameters ($\alpha_{\text{acc}}$ and $\alpha_{\text{mix}}$
respectively) and quantified to what extent a feature such as a planet or
vortex can interfere with VSI activity.

To achieve this, we first computed the total accretion and mixing parameters
using the full 3D structure of a snapshot ($\alpha^{\text{tot}}$), which
includes the contribution of nonaxisymmetric features. We then compared this
to the stress calculated using azimuthally averaged datasets
($\alpha^{\text{VSI}}$), taking advantage of the inherently axisymmetric
nature of the VSI.

By first comparing a typical 2D axisymmetric VSI model to a fully 3D one, we
found that the total accretion stress in 3D is slightly higher than in the 2D
case, consistent with previous studies (Manger & Klahr 2018). We attribute
this additional stress to the spiral arms that form in the 3D model, and
highlight that the VSI-induced stress is in fact lower, as can be seen in both
$\alpha_{\text{acc}}$ and $\alpha_{\text{mix}}$. This is not surprising, as
spiral arms can contribute strongly to accretion (Goodman & Rafikov 2001)
while the VSI remains the main driver of vertical motion in the disk. As a
result, we find that
$\alpha_{\text{acc}}^{\text{VSI}}<\alpha_{\text{acc}}^{\text{2D}}<\alpha_{\text{acc}}^{\text{tot}}$
and
$\alpha_{\text{mix}}^{\text{VSI}}\approx\alpha_{\text{mix}}^{\text{tot}}<\alpha_{\text{mix}}^{\text{2D}}$.

In models with embedded planets, the spiral wakes excited by the planet can
strongly affect the velocity field in the disk, generating an effective
accretion stress. The result, similar to the 3D case without a planet, is
slightly weaker vertical mixing while the radial accretion stress becomes
higher, in line with the expectations from theory (Goodman & Rafikov 2001).

In addition, due to the nearly inviscid disk conditions, a planet with
$M_{\text{p}}\gtrsim M_{\text{th}}$ will open a gap in the gas (Crida et al.
2006) while also generating strong zonal flows near its vicinity (Bi et al.
2021), as well as large-scale vortices at its gap edge (Paardekooper et al.
2010; Rometsch et al. 2021) that launch waves of similar magnitude as they
decay over time. All of these processes can, to different extents, disrupt the
velocity field in such a way that the contribution of the VSI to the total
accretion stress is significantly weaker in the outer disk for a massive
planet. In the inner disk, due to the absence of massive vortices, the
stronger vertical shear rate, and the shorter growth timescale of the VSI
($\propto\Omega_{\text{K}}^{-1}$) in contrast to the “disturbance timescale”
by the planet’s spirals as they propagate radially inward
($\propto\Omega_{\text{p}}^{-1}$), the VSI remains active at levels similar to
those expected by axisymmetric models, without embedded planets.

Additional mechanisms that a planet excites and can drive vertical motion,
such as the spiral wave instability (SWI, Bae et al. 2016a, b) or buoyancy
waves (Zhu et al. 2015; McNally et al. 2020), could further hinder the
contribution of the VSI to the total accretion and mixing turbulent
parameters. In our case, the short cooling timescale of $\beta=0.01$ that we
impose should minimize the excitation of buoyancy waves, while our grid
resolution might be too low to reliably invoke the SWI as an alternative
driver of turbulence. A study on the effect of either mechanism on the
activity of the VSI is unfortunately out of the scope of our project, but is
worth investigating in the future.

Vortices can, in principle, form without the necessity of the presence of a
planet in VSI-active disks (Richard et al. 2016; Manger & Klahr 2018; Flock et
al. 2020). A large-scale vortex can, similar to a planet, launch spiral arms
that corotate with its orbit. Such spirals will also drive an accretion stress
with $\alpha_{\text{acc}}\sim 10^{-4}$ for our disk parameters (Larson 1990),
and can subject the gas to vertical motion through one or more of the
mechanisms outlined above. This results in a weaker VSI stress at the
immediate vicinity of the vortex, that recovers at sufficient distances
($\approx 4H$ from the eye of the vortex for our models). While this would
very likely render the VSI undetectable and the mm dust layer considerably
flatter at the vicinity of the vortex, the effect is weaker than in the case
of an equally massive planet.

In summary, our findings suggest that the VSI will be severely hindered in
disks with massive embedded planets, and considerably weaker in the vicinity
of large-scale vortices. This has implications for the height of the mm dust
layer (Dullemond et al. 2022), which is commonly the target of observations
using the ALMA array (e.g., Andrews et al. 2018), as well as for the
observability of the instability using ALMA detections of molecular line
emission (Barraza-Alfaro et al. 2021, 2022).

###### Acknowledgements.

AZ and RPN are grateful for the guidance and friendship of Willy Kley, and
dedicate this paper to his memory. We thank the anonymous referee for their
helpful suggestions and comments. The authors acknowledge support by the High
Performance and Cloud Computing Group at the Zentrum für Datenverarbeitung of
the University of Tübingen, and the state of Baden-Württemberg through bwHPC.
This research utilized Queen Mary’s Apocrita HPC facility, supported by QMUL
Research-IT (http://doi.org/10.5281/zenodo.438045). This work was performed
using the DiRAC Data Intensive service at Leicester, operated by the
University of Leicester IT Services, which forms part of the STFC DiRAC HPC
Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding
via STFC capital grants ST/K000373/1 and ST/R002363/1 and STFC DiRAC
Operations grant ST/R001014/1. DiRAC is part of the National e-Infrastructure.
AZ and RPN are supported by STFC grant ST/P000592/1, and RPN is supported by
the Leverhulme Trust through grant RPG-2018-418. All plots in this paper were
made with the Python library matplotlib (Hunter 2007).

## References

  * Andrews et al. (2018) Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41
  * Bae et al. (2016a) Bae, J., Nelson, R. P., & Hartmann, L. 2016a, ApJ, 833, 126
  * Bae et al. (2016b) Bae, J., Nelson, R. P., Hartmann, L., & Richard, S. 2016b, ApJ, 829, 13
  * Balbus & Hawley (1991) Balbus, S. A. & Hawley, J. F. 1991, ApJ, 376, 214
  * Balbus & Papaloizou (1999) Balbus, S. A. & Papaloizou, J. C. B. 1999, ApJ, 521, 650
  * Barraza-Alfaro et al. (2022) Barraza-Alfaro, M., Flock, M., & Doe, J. 2022, A&A
  * Barraza-Alfaro et al. (2021) Barraza-Alfaro, M., Flock, M., Marino, S., & Pérez, S. 2021, arXiv e-prints, arXiv:2106.01159
  * Bi et al. (2021) Bi, J., Lin, M.-K., & Dong, R. 2021, ApJ, 912, 107
  * Crida et al. (2006) Crida, A., Morbidelli, A., & Masset, F. 2006, Icarus, 181, 587
  * Cui & Bai (2020) Cui, C. & Bai, X.-N. 2020, ApJ, 891, 30
  * Cui & Latter (2022) Cui, C. & Latter, H. N. 2022, MNRAS, 512, 1639
  * de Val-Borro et al. (2007) de Val-Borro, M., Artymowicz, P., D’Angelo, G., & Peplinski, A. 2007, A&A, 471, 1043
  * de Val-Borro et al. (2006) de Val-Borro, M., Edgar, R. G., Artymowicz, P., et al. 2006, MNRAS, 370, 529
  * Dullemond et al. (2018) Dullemond, C. P., Birnstiel, T., Huang, J., et al. 2018, ApJ, 869, L46
  * Dullemond et al. (2022) Dullemond, C. P., Ziampras, A., Ostertag, D., & Dominik, C. 2022, arXiv e-prints, arXiv:2210.13413
  * Flaherty et al. (2020) Flaherty, K., Hughes, A. M., Simon, J. B., et al. 2020, ApJ, 895, 109
  * Flock et al. (2017) Flock, M., Nelson, R. P., Turner, N. J., et al. 2017, ApJ, 850, 131
  * Flock et al. (2020) Flock, M., Turner, N. J., Nelson, R. P., et al. 2020, ApJ, 897, 155
  * Flores-Rivera et al. (2020) Flores-Rivera, L., Flock, M., & Nakatani, R. 2020, A&A, 644, A50
  * Fricke (1968) Fricke, K. 1968, ZAp, 68, 317
  * Gammie (1996) Gammie, C. F. 1996, ApJ, 457, 355
  * Gammie (2001) Gammie, C. F. 2001, ApJ, 553, 174
  * Goldreich & Schubert (1967) Goldreich, P. & Schubert, G. 1967, ApJ, 150, 571
  * Goodman & Rafikov (2001) Goodman, J. & Rafikov, R. R. 2001, ApJ, 552, 793
  * Harten (1983) Harten, A. 1983, Journal of Computational Physics, 49, 357
  * Huang et al. (2018) Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2018, ApJ, 869, L42
  * Hunter (2007) Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90
  * Klahr & Kley (2006) Klahr, H. & Kley, W. 2006, A&A, 445, 747
  * Kley & Nelson (2012) Kley, W. & Nelson, R. P. 2012, ARA&A, 50, 211
  * Larson (1990) Larson, R. B. 1990, MNRAS, 243, 588
  * Li et al. (2005) Li, H., Li, S., Koller, J., et al. 2005, ApJ, 624, 1003
  * Lin & Youdin (2015) Lin, M.-K. & Youdin, A. N. 2015, ApJ, 811, 17
  * Lynden-Bell & Pringle (1974) Lynden-Bell, D. & Pringle, J. E. 1974, MNRAS, 168, 603
  * Lyra & Umurhan (2019) Lyra, W. & Umurhan, O. M. 2019, PASP, 131, 072001
  * Manger & Klahr (2018) Manger, N. & Klahr, H. 2018, MNRAS, 480, 2125
  * Manger et al. (2020) Manger, N., Klahr, H., Kley, W., & Flock, M. 2020, MNRAS, 499, 1841
  * Manger et al. (2021) Manger, N., Pfeil, T., & Klahr, H. 2021, MNRAS, 508, 5402
  * Masset (2000) Masset, F. 2000, A&AS, 141, 165
  * Masset & Benítez-Llambay (2016) Masset, F. S. & Benítez-Llambay, P. 2016, ApJ, 817, 19
  * McNally et al. (2020) McNally, C. P., Nelson, R. P., Paardekooper, S.-J., Benítez-Llambay, P., & Gressel, O. 2020, MNRAS, 493, 4382
  * Mignone et al. (2007) Mignone, A., Bodo, G., Massaglia, S., et al. 2007, The Astrophysical Journal Supplement Series, 170, 228
  * Mignone et al. (2012) Mignone, A., Flock, M., Stute, M., Kolb, S. M., & Muscianisi, G. 2012, A&A, 545, A152
  * Miranda & Rafikov (2020) Miranda, R. & Rafikov, R. R. 2020, ApJ, 904, 121
  * Nelson et al. (2013) Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610
  * Ogilvie & Lubow (2002) Ogilvie, G. I. & Lubow, S. H. 2002, MNRAS, 330, 950
  * Paardekooper et al. (2010) Paardekooper, S.-J., Lesur, G., & Papaloizou, J. C. B. 2010, ApJ, 725, 146
  * Pfeil & Klahr (2021) Pfeil, T. & Klahr, H. 2021, ApJ, 915, 130
  * Picogna et al. (2018) Picogna, G., Stoll, M. H. R., & Kley, W. 2018, A&A, 616, A116
  * Rafikov (2002) Rafikov, R. R. 2002, ApJ, 569, 997
  * Rafikov (2016) Rafikov, R. R. 2016, ApJ, 831, 122
  * Richard et al. (2016) Richard, S., Nelson, R. P., & Umurhan, O. M. 2016, MNRAS, 456, 3571
  * Rometsch et al. (2021) Rometsch, T., Ziampras, A., Kley, W., & Béthune, W. 2021, A&A, None, None
  * Shakura & Sunyaev (1973) Shakura, N. I. & Sunyaev, R. A. 1973, A&A, 500, 33
  * Stoll & Kley (2014) Stoll, M. H. R. & Kley, W. 2014, A&A, 572, A77
  * Stoll et al. (2017a) Stoll, M. H. R., Kley, W., & Picogna, G. 2017a, A&A, 599, L6
  * Stoll et al. (2017b) Stoll, M. H. R., Picogna, G., & Kley, W. 2017b, A&A, 604, A28
  * Svanberg et al. (2022) Svanberg, E., Cui, C., & Latter, H. N. 2022, MNRAS, 514, 4581
  * Tassoul (1978) Tassoul, J.-L. 1978, Theory of rotating stars
  * Van Leer (1974) Van Leer, B. 1974, Journal of computational physics, 14, 361
  * Villenave et al. (2022) Villenave, M., Stapelfeldt, K. R., Duchêne, G., et al. 2022, ApJ, 930, 11
  * Zhang et al. (2018) Zhang, S., Zhu, Z., Huang, J., et al. 2018, ApJ, 869, L47
  * Zhu et al. (2015) Zhu, Z., Dong, R., Stone, J. M., & Rafikov, R. R. 2015, ApJ, 813, 88
  * Ziampras et al. (2020a) Ziampras, A., Ataiee, S., Kley, W., Dullemond, C. P., & Baruteau, C. 2020a, A&A, 633, A29
  * Ziampras et al. (2020b) Ziampras, A., Kley, W., & Dullemond, C. P. 2020b, A&A, 637, A50

##  Appendix A Measuring turbulent stress

To measure turbulent stress we followed a four-step process based on Balbus &
Papaloizou (1999). In the section below, we denote a quantity $q$ averaged
with respect to $y$ as $\langle{q}\rangle_{y}$.

We begin by first computing an axisymmetric, time-averaged disk state for all
quantities $\rho$, $\bm{u}$ and $P$ using 150–500 snapshots depending on
model, taken at intervals of $1/P_{0}$. For a quantity $q$, we have

$\langle{q}\rangle_{\phi
t}=\frac{1}{N_{\text{files}}}\sum\limits_{i=i_{\text{start}}}^{i_{\text{end}}}\langle{q^{i}}\rangle_{\phi},\quad\langle{q^{i}}\rangle_{\phi}=\frac{1}{2\pi}\int\limits_{0}^{2\pi}q^{i}\,\text{d}\phi,$
(12)

where $i_{\text{start}}$ and $i_{\text{end}}$ denote the indices of the first
and last snapshots used, and
$N_{\text{files}}=i_{\text{end}}-i_{\text{start}}+1$ is the number of
snapshots.

We then compute the $R\phi$ and $z\phi$ components of the Reynolds stress
tensor $\overline{\textbf{T}}(\bm{r})$

$T_{x\phi}^{i}(\bm{r})=\rho^{i}\cdot\left[u_{x}^{i}-\langle{u_{x}}\rangle_{\phi
t}\right]\cdot\left[u_{\phi}^{i}-\langle{u_{\phi}}\rangle_{\phi
t}\right],\qquad x\in\\{R,z\\},$ (13)

and obtain two 3D data cubes $\langle{T_{R\phi}}\rangle_{t}$ and
$\langle{T_{z\phi}}\rangle_{t}$ by averaging over the same snapshots that we
used to compute the mean, background field. We then compute the azimuthal
average and apply a radial smoothing over a window of a local scale height
$H(R)$ on both datasets:

$T^{\prime}_{x\phi}(R,z)=\frac{1}{\Delta R}\int\limits_{R-\Delta
R/2}^{R+\Delta
R/2}\langle{T_{x\phi}}\rangle_{t\phi}\,\text{d}R^{\prime},\qquad\Delta
R=H(R).$ (14)

This is equivalent to a radial rolling average, and can be easily implemented
with a first-degree Savitzky–Golay filter (available as
scipy.signal.savgol_filter in Python).

To compute $\alpha_{\text{acc}}(R)$ we integrate $T^{\prime}_{R\phi}$
vertically and finally obtain

$\alpha_{\text{acc}}(R)=-W_{R\phi}\cdot\left(R\frac{\text{d}{\Omega}}{\text{d}{R}}\Sigma
c_{\text{s}}H\right)^{-1},\qquad
W_{R\phi}=\int\limits_{z_{\text{min}}}^{z_{\text{max}}}T^{\prime}_{R\phi}\,\text{d}z,$
(15)

where $\Sigma$ is computed using the time- and azimuthally-averaged
$\langle{\rho}\rangle_{t\phi}$. We note that, since the mean azimuthal
velocity and disk temperature do not change significantly due to the VSI
($\Delta u_{\phi}/{u_{\phi}}_{0}\sim 10^{-2}$, $\Delta T/T_{0}\sim 10^{-6}$),
one can also instead use the initial, analytical profiles of $\Omega$,
$c_{\text{s}}$ and $H$ in these calculations. However, since the disk is
constantly accreting (i.e., losing mass in our models), the mean density
$\langle{\rho}\rangle_{t\phi}$ is more appropriate than the initial profile
$\rho^{\text{eq}}$ (see Eq. (5)) when calculating $\Sigma$.

Finally, we compute $\alpha_{\text{mix}}$ by replacing the viscous stress
tensor component in Eq. (10) of Stoll et al. (2017b) with $T_{z\phi}^{\prime}$
and find

$\alpha_{\text{mix}}=\frac{T_{z\phi}^{\prime}}{\langle{P}\rangle_{\phi
t}}\cdot\left(\frac{|q|h}{2}\frac{z}{H}\right)^{-1}.$ (16)

We note that this only provides an estimate of $\alpha_{\text{mix}}$, as it is
not backed by the accretion model by Balbus & Papaloizou (1999). A more
appropriate method of calculating $\alpha_{\text{mix}}$ would involve
scattering dust particles in a numerical model and measuring their scale
height as they are lifted by the gas (e.g., Dullemond et al. 2022).
Nevertheless, our method provides compatible results with those of Dullemond
et al. (2022) for our models.

To isolate the contribution of the VSI, we simply replace $\rho^{i}$ and
$u^{i}$ in Eq. (13) with their azimuthal averages
$\langle{\rho^{i}}\rangle_{\phi}$ and $\langle{u^{i}}\rangle_{\phi}$
respectively. This means that azimuthal averaging in Eq. (14) is not
necessary, as $\langle{T_{x\phi}}\rangle_{t}$ is already axisymmetric in this
case.

##  Appendix B Resolution study

Previous studies have shown that a resolution of 16 cells per scale height
(cps) is enough to resolve the growth and saturation of the VSI and recover a
value of $\alpha\sim 10^{-4}$. Nevertheless, to verify that the VSI is
resolved appropriately in our models, we conduct a short resolution study
where we measure the growth rate (Fig. 12) of the instability using 2D
axisymmetric models.

After a local linear stability analysis we found that the fastest-growing VSI
mode has a wavelength of approximately $0.1\,H$, which requires a radial
resolution of 20 cps in order to be resolved over 2 cells. To maintain
compatibility with previous studies, we then run models at 8, 16, 32 and 48
cps in both $r$ and $\theta$ directions.

We find that a resolution of 16 cps captures a growth rate of $\approx
0.58/\text{orbit}$, similar to that reported by Stoll & Kley (2014), with
higher-resolution models showing very similar results (see Fig. 12).

Figure 12: Growth rates of the VSI as a function of resolution in cells per
scale height (cps) in our axisymmetric models. The two shaded bands denote a
growth rate of 0.58/orbit and 0.13/orbit for the steep and shallow bands
respectively, in agreement with Stoll & Kley (2014).

##  Appendix C Issues with the inner boundary

In our 3D runs shown in Sect. 3.2 we noted the presence of a parasitic
instability in the inner radial boundary, likely owing to a problematic
interface between the domain and the wave damping zone. As shown in Fig. 13.1,
nonaxisymmetric features in the form of over- and under-densities form on the
inner boundary edge at $R=0.4\,R_{0}$ and launch spirals akin to the vortices
discussed in Sects 3.2 and 3.4. These spirals permeate the active domain,
contributing slightly to the total Reynolds stress and likely inhibiting VSI
activity to an extent. This phenomenon has been seen before (e.g., Barraza-
Alfaro et al. 2021, see Fig. 3 therein), and might relate to the weaker
vertical velocities observed near the inner boundary for low-mass planets in
Stoll et al. (2017a).

In an attempt to understand the origin of these nonaxisymmetric features, we
calculate the potential vorticity at the midplane

$\zeta_{\text{mid}}=\frac{(\nabla\times\bm{u}_{\text{mid}})\cdot\hat{z}}{\rho_{\text{mid}}}$
(17)

and plot it for various timestamps in Fig. 14.2. We find that vortex-like
features start forming as early as 10 $P_{0}$ into the simulation at both the
inner boundary at $R=0.4\,R_{0}$ but also on the interface between the damping
zone and the active domain at $0.5\,R_{0}$. While these features grow in size
in the absence of a planet (see panels _a_ – _e_ in Figs. 13.1 & 14.2), the
planet’s inner spiral arms break up these overdensities, restabilizing the
inner boundary to an extent and weakening the spiral arms they launch in the
process (see panels _f_ – _h_ in the same figures). This would then explain—in
part—why the VSI can remain active closer to the inner radial boundary in both
our runs and in those by Stoll et al. (2017a).

Figure 13: Time evolution of deviations in the midplane density
$\rho_{\text{mid}}$ with respect to its azimuthal average, showing the
development of over- and under-dense areas near the inner radial boundary
edge. The time shown on each panel refers to the time elapsed since the disk
has been expanded to 3D (see Fig. 1), and the planet grows between panels _e_
– _g_. The nonaxisymmetric features grow over time in panels _a_ – _e_ , but
are dampened within a few orbits after the planet is introduced in the domain
(panel _f_). Panel _h_ shows that these features no longer grow in the
presence of the embedded planet. Vertical dotted lines mark the end of the
wave-damping zone. Figure 14: Time evolution of the deviations in the
potential vorticity in the midplane according to Eq. (17) with respect to its
azimuthal average. The figure is annotated similar to Fig. 13.1. In
particular, panel _a_ shows that nonaxisymmetric features appear at the inner
radial boundary edge at $R=0.4\,R_{0}$ as well as at the edge of the wave-
damping zone at $0.5\,R_{0}$.