© The Authors 2023

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 ≳10au 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 α ~ 10⁻⁴−10⁻³ (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 M_⊕. They have also shown that, even in this high-mass case, the disk maintains an accretion turbulence parameter α ~ 10⁻⁴. 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_p = 100 M_⊕, 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 ρ, pressure P and velocity field u orbiting around a star with mass M_⋆ = 1 M_⊙ $$\frac{\partial \rho }{\partial t}+u\cdot \nabla \rho =-\rho \nabla \cdot u$$(1) $$\frac{\partial u}{\partial t}+\left(u\cdot \nabla \right)u=-\frac{1}{\rho }\nabla P-\nabla \text{Φ}$$(2) $$\frac{\partial e}{\partial t}+u\cdot \nabla e=-\gamma e\nabla \cdot u+S,$$(3)

where Φ = −GM_⋆/r is the gravitational potential of the star at distance r, e = P/(γ − 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 =\raisebox{1ex}{$7$}\!\left/ \!\raisebox{-1ex}{$5$}\right.$ and a mean molecular weight µ = 2.35 such that e = ρc_vT, with $c_{\text{v}}=\frac{ℛ}{\mu \left(\gamma -1\right)}$ being its heat capacity at constant volume. The gravitational constant and gas constant are indicated with G and ℛ, 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, ϕ, z} where $r=\sqrt{R^2+z^2}$. We can derive a hydrodynamic equilibrium state for an axisymmetric $\left(\raisebox{1ex}{$\partial $}\!\left/ \!\raisebox{-1ex}{$\partial \varphi $}\right.=0\right)$, nonaccreting disk $\left(u=u_{\varphi }\widehat{\varphi }\right)$ with a vertically isothermal temperature profile $\left(\raisebox{1ex}{$\partial T$}\!\left/ \!\raisebox{-1ex}{${\partial }_z$}\right.=0\right)$. If we further assume that the gas density at the midplane ρ_mid as well as the temperature both follow power law profiles in the radial direction R with $$\begin{array}{cc}{\rho }_{\text{mid}}\left(R\right)={\rho }_{z=0}={\rho }_0{\left(\frac{R}{R_0}\right)}^p,& T\left(R\right)\end{array}=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) $${\rho }^{\text{eq}}\left(R,z\right)={\rho }_{\text{mid}}\left(R\right)\exp \left[\frac{G{M}_{\star }}{c_{\text{s}}^2}\left(\frac{1}{r}-\frac{1}{R}\right)\right],$$(5) $$u_{\varphi }^{\text{eq}}\left(R,z\right)=R{\text{Ω}}_{\text{K}}{\left[1+\left(p+q\right){\left(\frac{H}{R}\right)}^2+q\left(1-\frac{R}{r}\right)\right]}^{\mathrm{½}},$$(6)

where ${\text{Ω}}_{\text{K}}\left(R\right)=\sqrt{G{M}_{\star }/R^3}$ is the Keplerian angular velocity and $H\left(R\right)=c_{\text{s}}/{\text{Ω}}_{\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 β-cooling approach (e.g., Gammie 2001), where the cooling timescale is given by t_cool = β/Ω_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 }{\text{Ω}}_{\text{K}}.$$(7)

For simplicity, we adopt a constant β = 0.01 in our models. This would correspond to a reference radius R₀ ~ 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 (β → 0; e.g., Stoll et al. 2017b; Manger & Klahr 2018; Barraza-Alfaro et al. 2021). Our choice of β 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)–(6) with $p=-\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$, 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 $\text{Σ}\left(R\right)\propto R^{-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}$. Our models are scale-free in terms of the reference radius R₀ and surface density Σ₀ = Σ(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 -\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$, h ≈ 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 α parameter (Shakura & Sunyaev 1973) on the order of a few 10⁻⁴ 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,θ,ϕ} 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 ϕ ϵ [0,2π), a radial range r ϵ [0.4,2.5] R₀ and a vertical extent z ϵ [−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 × N_θ × N_ϕ = 600 × 128 × 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 ~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 ϵ [0.4,0.5] ∪ [2.1,2.5] R₀ 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_p = 3 × 10⁻⁴ M_⋆ (100 M_⊕ 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}}=\frac{c_{\text{s}}^3}{G{\text{Ω}}_{\text{K}}}|{}_{R_{\text{p}}}\approx 1M_{\text{J}}{\left(\frac{h_{\text{p}}}{0.1}\right)}^3\frac{M_{\star }}{M_{\odot }},$$(8)

which, for our choices of aspect ratio h = {0.05,0.1} evaluates to {0.125,1} M_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_p = 0.3 M_th) to the nonlinear regime (h = 0.05, M_p = 2.4 M_th).

2.3 Model timeline

Taking advantage of the axisymmetric nature of the VSI we first construct 2D models of the VSI on the {r, θ} plane. These models, while axisymmetric, involve all 3 velocity components. We then carry out a resolution analysis, in order to measure the turbulent viscosity a 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_2D≈2.5^3/2 P₀ ≈ 120 orbits at R = R₀, we continue until t_2D = 1000 P₀ to ensure the disk reaches a quasi-steady state while also allowing a reasonably long timeframe of 500 P₀ (between 500–1000 P₀, with an output frequency of 1/P₀) 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) and (6), with all velocity components perturbed by 1% c_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 ϕ-axis while also seeding the velocity with a further 1% c_s perturbation to break the axisymmetry and allow the formation of azimuthal structure. The 3D models are then evolved for a further 200 P₀, and related results are presented in Sect. 3.2.

Finally, a planet is inserted at R = R₀ 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 P₀ , 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.

Fig. 1 Timeline for our 3D models with embedded planets. Models in the blue, orange and green phases are discussed in Sects. 3.1–3.3 respectively. The planet finishes growing by t = t0 + 100 P0 (purple dotted line).

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).

Consistent with previous studies, we compute a VSI-driven accretion stress α_acc ≈ 1–3 × 10⁻⁴ and a vertical turbulent mixing parameter α_mix ≈ 0.1–0.3 for our choice of disk parameters. To measure α 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 ${\overline{u}}_R=0$) and further adjusting the result by a factor of $\raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$3$}\right.$. Our method is described in detail in Appendix A, and our results for α_acc and α_mix as functions of R and z respectively are shown in Fig. 3.

Fig. 2 Radial (uR), vertical (uz), and azimuthal (uϕ) gas velocity components in our fiducial axisymmetric model, highlighting the vertically elongated sheets of motion generated by the VSI. The Keplerian backgroun uK has been subtracted from uϕ, and all components have been normalized to the local sound speed.

Fig. 3 Radial and vertical profiles of the turbulent accretion parameter αacc and the turbulent mixing parameter αmix in our fiducial axisymmetric model. We find a typical αacc ≈ 1.5 × 10−4 and αmix ≈ 0.15−0.3, in agreement with previous results. The faint curve on the upper panel refers to αacc without radial smoothing through Eq. (A.3).

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 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₀. The gas shows azimuthal structure, with spiral arms permeating the disk and two vortices orbiting at R ≈ l.l5 R₀· Nevertheless, the disk is still VSI-active, as indicated by the largely axisymmetric velocity component u_z with $u_z^{\text{mid}}\approx 0.2c_{\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[{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\frac{\rho }{\left(\nabla \times u\right)\cdot \widehat{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 ≈ 0.4R₀) 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.

We now measure turbulent stress by repeating the procedure in Sect. 3.1, this time computing α 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 α_acc compared to our 2D models. However, the disk shows a weaker vertical mixing parameter α_mix in the region R ~ R₀, which is likely more strongly affected by the vortex present there, as well as in the region R ≲ 0.7R₀, which is dominated by the spirals launched by the vortices that formed near the inner boundary.

By assuming that, approximately, the total stress is α^tot = α^VSI + α^spiral, we attempt to isolate the contribution of the VSI to the total turbulent stress by recomputing α after first averaging our snapshots along the ϕ 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 α 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 (δ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 (δu_r, δu_ϕ) 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}\text{for\hspace{0.17em}}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.

Fig. 4 Highlights of our 3D fiducial model. Left: surface density deviations in the disk compared to the azimuthally averaged surface density $\overline{\text{Σ}}$, 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 ϖ0, highlighting two vortices orbiting at R ≈ 1.15 R0

Fig. 5 Turbulent parameters αacc and α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, 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.4R0).

3.3 Plаnet-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 α^tot and α^VSI in order to highlight the effect of the planet on the VSI. We then collect our estimates on α_acc and α_mix with and without the contribution of nonaxisymmetric features for all 3D runs in Table 1.

Fig. 6 Surface density deviations and vertical velocity maps at the mid-plane for our disk with an embedded subthermal-mass planet (h = 0.1, Mp = 0.3 Mth) at R = R0 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.

3.3.1 Quasi-linear regime: h = 0.1 → M_p = 0.3 M_th

Figure 6 shows an example of the typical disk structure around a subthermal-mass planet with M_p = 0.3 M_th, taken after 100 planetary orbits in our model with h = 0.1, M_p = 100 M_⊕. 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 α_acc and vertical mixing parameter α_mix by first considering the full 3D structure of the disk (which includes the effects of spiral arms and results in a total stress α^tot), and then comparing our findings to an estimate of the VSI-driven stress α^VSI by first averaging all quantities along the ϕ 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{tot}}$ (dashed orange curve) by a factor of ≈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 α_mix on the right panel of Fig. 7, as the VSI is the main driver of vertical motion (and thus α_mix).

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 β = 0.01. An additional contribution to α_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 α_acc ~ 10⁻⁴ 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_p = 13 M_⊕, yielding 0.3 M_th, resulted in very similar behavior, further highlighting that M_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.

Fig. 7 Turbulent parameters αacc and 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, α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 → M_p = 2.4 M_th

Contrary to the disk model with h = 0.1 discussed in the previous paragraph, a planet with M_p = 100 M_⊕ amounts to 2.4M_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 (α_acc ≲ 5 × 10⁻⁵). 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_p = 1 M_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≳R₀, 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.

We further show that the VSI remains active in the inner disk (R ≲ R₀), 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 Ω_K, as well as the faster growth rate of the instability in the inner disk $\left(\propto {\Omega }_{\text{K}}^{-1}\right)$ when compared to the timescale over which spiral arms, which corotate with the planet $\left(\propto {\Omega }_{\text{p}}^{-1}\right)$, 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 ≳ R_p). We therefore conclude that the VSI is very unlikely to coexist with massive planets (M_p ≳ M_th), consistent with the qualitative findings of Stoll et al. (2017a) and the synthetic observations by Barraza-Alfaro (in prep.).

Fig. 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 Σ0 and ϖ0.

Fig. 9 Turbulence parameters αacc and αmix similar to Fig. 7, but for the superthermal-mass case with h = 0.05, Mp = 2.4Mth. The surface density snapshot of Fig. 8 is also shown on panel a, to highlight features present at different radii. The parameter α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 Rp ± 2.5 RHill

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₀, and place a blob of gas with total mass 1 M_J at the outer gap edge using a Gaussian distribution centered at R = 1.25 R₀, ϕ = π (see top panel of Fig. 10) $$\rho \left(R,\varphi \right)={\rho }^{\text{eq}}\times f_{\text{gap}}\left(R\right)\times f_{\text{vortex}}\left(R,\varphi \right),$$(10)

with ƒ_gap and ƒ_voгtex motivated by our embedded-planet models in Sect. 3.3: $$\begin{array}{l}{f}_{\text{gap}}\left(R\right)=1-0.95\exp \left[-\frac{{\left(R/R_0-1\right)}^2}{0.05}\right]\hfill \\ \begin{array}{c}{f}_{\text{vortex}}\left(R,\varphi \right)=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{\varphi -\pi }{0.35}\right)}^2\right].\end{array}\hfill \end{array}$$(11)

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 Ω(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 shows the results of our data reduction for these models. Consistent with our findings in Sects. 3.2 and 3.3, we find that α_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, α_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 α in any noticeable way.

Fig. 10 Surface density, vertical velocity and vortensity snapshots for our run with an embedded vortex after 400 orbits at R0. To eliminate the stripes of vortensity deviations (see Fig. 4) and highlight the ϖ minimum about the vortex, we normalize ϖ to its azimuthal median 〈ϖ〉.

Fig. 11 Turbulent αacc and αmix in our vortex runs, similar to Fig. 5. VSI activity is significantly lower near the radial vicinity of the vortex (R≈1.3R0).

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 ≳ 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 (in prep.) 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_z as a tracer for VSI activity, the instability should be weaker in the outer disk for M_p = 100 M_⊕, or 2.4M_th, in their models, but did not measure the vertical mixing parameter α_mix in their simulations and instead focused on accretion stress (α_acc).

Our models, for which we provide quantitative values of α_mix, support the findings of both studies and enrich them further by providing insight into the effect of a large-scale vortex on α_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 α_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 (in prep.), 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 ~ 0.1, q ≈ −0.5, such that α^VSI ~ 10⁻⁴) 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_p = 2.4 M_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_P = 1 M_J. Given the nearly-locally isothermal disk conditions (β = 0.01) and the lack of VSI activity in their vicinity resulting in a nearly inviscid disk with α ~ max(α_num, 10⁻⁷) ≪ 10⁻⁴, 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 α ≈ 5 × 10⁻⁵ (see top panel of Fig. 9), where it decays over ≈10³ local orbits at R ≈ 0.7 R₀. This lifetime is consistent with the estimates of Rometsch et al. (2021) for a disk with α = 10⁻⁴, β = 0.01 and h = 0.05, suggesting that the turbulence driven by the VSI is compatible with 2D laminar-α models. This last point complements the findings of Stoll et al. (2017a), who showed that α-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 (≳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_p = M_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 non-axisymmetric 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 (α_acc and α_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 (α^tot), which includes the contribution of nonaxisymmetric features. We then compared this to the stress calculated using azimuthally averaged datasets (α^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 α_acc and α_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_P ≳ M_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 $\left(\propto {\Omega }_{\text{K}}^{-1}\right)$ in contrast to the “disturbance timescale” by the planet’s spirals as they propagate radially inward $\left(\propto {\Omega }_{\text{p}}^{-1}\right)$, 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 β = 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 α_acc ~ 10⁻⁴ 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 (≈4 H 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): Barraza-Alfaro (in prep.).

Acknowledgements

A.Z. and R.P.N. 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. A.Z. and R.P.N. are supported by STFC grant ST/P000592/1, and R.P.N. 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).

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 〈q〉_y.

We begin by first computing an axisymmetric, time-averaged disk state for all quantities ρ, u and P using 150–500 snapshots depending on model, taken at intervals of 1 / P₀ . For a quantity q, we have $$\begin{array}{cc}{\langle q\rangle }_{\varphi t}=\frac{1}{N_{\text{files}}}{\displaystyle \sum _{i=i_{\text{start}}}^{i_{\text{end}}}{\langle q^i\rangle }_{\varphi },}& {\langle q^i\rangle }_{\varphi }=\frac{1}{2\pi }{\displaystyle \underset{0}{\overset{2\pi }{\int }}{q}^i\text{d}\varphi }\end{array},$$(A.1)

where i_start and i_end denote the indices of the first and last snapshots used, and N_files = i_end − i_start + 1 is the number of snapshots.

We then compute the Rϕ and zϕ components of the Reynolds stress tensor $\overline{\text{T}}\left(r\right)$ $$\begin{array}{cc}{T}_{x\varphi }^i\left(r\right)={\rho }^i\cdot \left[u_x^i-{\langle u_x\rangle }_{\varphi t}\right]\cdot \left[u_{\varphi }^i-{\langle u_{\varphi }\rangle }_{\varphi t}\right],& x\in \left\{R,z\right\}\end{array},$$(A.2)

and obtain two 3D data cubes 〈T_Rϕ〉_t and 〈T_zϕ〉_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: $$\begin{array}{cc}{T^{\prime }}_{x\varphi }\left(R,z\right)=\frac{1}{\text{Δ}R}{\displaystyle \underset{R-\text{Δ}R/2}{\overset{R+\text{Δ}R/2}{\int }}{\langle T_{x\varphi }\rangle }_{t\varphi }}\text{d}{R}^{\prime },& \text{\hspace{1em}\hspace{1em}Δ}R\end{array}=H\left(R\right).$$(A.3)

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 α_acc(R) we integrate ${T^{\prime }}_{R\varphi }$ vertically and finally obtain $$\begin{array}{cc}{\alpha }_{\text{acc}}\left(R\right)=-W_{R\varphi }\cdot {\left(R\frac{\text{dΩ}}{\text{d}R}\text{Σ}{c}_{\text{s}}H\right)}^{-1},& W_{R\varphi }\end{array}={\displaystyle \underset{z_{\min }}{\overset{z_{\max }}{\int }}{T^{\prime }}_{R\varphi }}\text{d}z,$$(A.4)

where Σ is computed using the time- and azimuthally-averaged 〈ρ〉_tϕ. We note that, since the mean azimuthal velocity and disk temperature do not change significantly due to the VSI $\left(\Delta u_{\varphi }/u_{{\varphi }_0}\sim {10}^{-2},\Delta T/T_0\sim {10}^{-6}\right)$, one can also instead use the initial, analytical profiles of Ω, c_s and H in these calculations. However, since the disk is constantly accreting (i.e., losing mass in our models), the mean density 〈ρ〉_tϕ is more appropriate than the initial profile ρ^eq (see Eq. (5)) when calculating Σ.

Finally, we compute α_mix by replacing the viscous stress tensor component in Eq. (10) of Stoll et al. (2017b) with ${T^{\prime }}_{z\varphi }$ and find $${\alpha }_{\text{mix}}=\frac{{T^{\prime }}_{z\varphi }}{{\langle P\rangle }_{\varphi t}}\cdot {\left(\frac{\left|q\right|h}{2}\frac{z}{H}\right)}^{-1}.$$(A.5)

We note that this only provides an estimate of α_mix, as it is not backed by the accretion model by Balbus & Papaloizou (1999). A more appropriate method of calculating α_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 ρⁱ and uⁱ in Eq. (A.2) with their azimuthal averages 〈ρⁱ〉_ϕ and 〈uⁱ〉_ϕ respectively. This means that azimuthal averaging in Eq. (A.3) is not necessary, as 〈T_xϕ〉_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 α~ 10⁻⁴. 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. B.1) 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.1H, 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 θ directions.

We find that a resolution of 16 cps captures a growth rate of ≈ 0.58/orbit, similar to that reported by Stoll & Kley (2014), with higher-resolution models showing very similar results (see Fig. B.1).

Fig. B.1 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. C.1, nonaxisymmetric features in the form of over- and under-densities form on the inner boundary edge at R = 0.4 R₀ 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{\left(\nabla \times u_{\text{mid}}\right)\cdot \widehat{z}}{{\rho }_{\text{mid}}}$$(C.1)

and plot it for various timestamps in Fig. C.2. We find that vortex-like features start forming as early as 10 P₀ into the simulation at both the inner boundary at R = 0.4 R₀ but also on the interface between the damping zone and the active domain at 0.5R₀. While these features grow in size in the absence of a planet (see panels a–e in Figs. C.1 & C.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).

Fig. C.1 Time evolution of deviations in the midplane density ρ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 (panelf). 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.

Fig. C.2 Time evolution of the deviations in the potential vorticity in the midplane according to Eq. (C.1) with respect to its azimuthal average. The figure is annotated similar to Fig. C.1. In particular, panel a shows that nonaxisymmetric features appear at the inner radial boundary edge at R = 0.4 R0 as well as at the edge of the wave-damping zone at 0.5 R0.

References

All Tables

All Figures

	Fig. 1 Timeline for our 3D models with embedded planets. Models in the blue, orange and green phases are discussed in Sects. 3.1–3.3 respectively. The planet finishes growing by t = t0 + 100 P0 (purple dotted line).
In the text

	Fig. 2 Radial (uR), vertical (uz), and azimuthal (uϕ) gas velocity components in our fiducial axisymmetric model, highlighting the vertically elongated sheets of motion generated by the VSI. The Keplerian backgroun uK has been subtracted from uϕ, and all components have been normalized to the local sound speed.
In the text

	Fig. 3 Radial and vertical profiles of the turbulent accretion parameter αacc and the turbulent mixing parameter αmix in our fiducial axisymmetric model. We find a typical αacc ≈ 1.5 × 10−4 and αmix ≈ 0.15−0.3, in agreement with previous results. The faint curve on the upper panel refers to αacc without radial smoothing through Eq. (A.3).
In the text

Fig. 4 Highlights of our 3D fiducial model. Left: surface density deviations in the disk compared to the azimuthally averaged surface density $\overline{\text{Σ}}$, 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 ϖ0, highlighting two vortices orbiting at R ≈ 1.15 R0

In the text

Fig. 5 Turbulent parameters αacc and α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, 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.4R0).

In the text

	Fig. 6 Surface density deviations and vertical velocity maps at the mid-plane for our disk with an embedded subthermal-mass planet (h = 0.1, Mp = 0.3 Mth) at R = R0 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.
In the text

	Fig. 7 Turbulent parameters αacc and 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, αacc increases but the overall mixing driven by the VSI is weaker by a factor of 2.
In the text

Fig. 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 Σ0 and ϖ0.

In the text

Fig. 9 Turbulence parameters αacc and αmix similar to Fig. 7, but for the superthermal-mass case with h = 0.05, Mp = 2.4Mth. The surface density snapshot of Fig. 8 is also shown on panel a, to highlight features present at different radii. The parameter α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 Rp ± 2.5 RHill

In the text

	Fig. 10 Surface density, vertical velocity and vortensity snapshots for our run with an embedded vortex after 400 orbits at R0. To eliminate the stripes of vortensity deviations (see Fig. 4) and highlight the ϖ minimum about the vortex, we normalize ϖ to its azimuthal median 〈ϖ〉.
In the text

	Fig. 11 Turbulent αacc and αmix in our vortex runs, similar to Fig. 5. VSI activity is significantly lower near the radial vicinity of the vortex (R≈1.3R0).
In the text

	Fig. B.1 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).
In the text

Fig. C.1 Time evolution of deviations in the midplane density ρ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 (panelf). 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.

In the text

	Fig. C.2 Time evolution of the deviations in the potential vorticity in the midplane according to Eq. (C.1) with respect to its azimuthal average. The figure is annotated similar to Fig. C.1. In particular, panel a shows that nonaxisymmetric features appear at the inner radial boundary edge at R = 0.4 R0 as well as at the edge of the wave-damping zone at 0.5 R0.
In the text