© The Authors 2024

1. Introduction

The vertical shear instability (VSI, Goldreich & Schubert 1967; Urpin & Brandenburg 1998) is one of several proposed hydrodynamical (HD) instabilities able to produce turbulence in weakly ionized zones of protoplanetary disks where magnetic instabilities are expected to be suppressed (see, e.g., Lesur et al. 2023). Even though numerical models predict VSI-induced turbulence to lead to negligible stellar accretion compared to current estimated rates (Hartmann et al. 1998; Manara et al. 2016), this phenomenon should still produce significant vertical stirring of dust grains (Stoll & Kley 2016; Flock et al. 2017), thereby increasing their collisional velocities and diffusivities and directly impacting their ability to coagulate and form planetesimals (Ormel & Cuzzi 2007; Johansen et al. 2014; Klahr & Schreiber 2020). Furthermore, the VSI has been shown in 3D simulations to lead to the formation of small- and large-scale anticyclonic vortices (Richard et al. 2016; Manger & Klahr 2018; Pfeil & Klahr 2021), which should locally concentrate dust particles, likely accelerating the formation of planetesimals (Johansen et al. 2007; Gerbig et al. 2020) and potentially explaining structures observed in disks (e.g., van der Marel et al. 2016; de Boer et al. 2021; Marr & Dong 2022). Determining where in a disk this instability can operate is thus a necessary step toward explaining current and future observations of protoplanetary disks, as well as improving current models of planet formation.

So far, no direct observations have been made that confirm the occurrence of VSI in protoplanetary disks. As proposed by Barraza-Alfaro et al. (2021), large-scale velocity perturbations produced by the VSI may be detectable in CO kinematics observations with the current capabilities of the Atacama Large Millimeter/submillimeter Array (ALMA). On the other hand, current ALMA observations of razor-thin millimeter-dust distributions in disks up to ∼100 au (Pinte et al. 2016; Doi & Kataoka 2021; Villenave et al. 2022) seem to contradict a possible VSI activity close to the disk midplane at such distances, as the vertical stirring of grains would lead to thicker dust distributions than observed (Dullemond et al. 2022). On the theoretical side, predictions of VSI-stable and unstable regions have been made (Malygin et al. 2017; Pfeil & Klahr 2019; Fukuhara et al. 2021, 2023) by locally applying the global stability criterion by Lin & Youdin (2015) and using simplified prescriptions to obtain temperature distributions. One of the goals of this work is then to improve on these investigations by considering disk models with realistic temperature distributions obtained via radiative transfer, taking into account our previous analysis of local and global stability criteria (Melon Fuksman et al. 2024, Paper I henceforth).

Besides the question of the disk stability, little is known about the mechanisms that lead to the saturation of the VSI. In Latter & Papaloizou (2018), it was shown via a local Boussinesq analysis of nonlinear perturbations that, once the VSI modes reach high enough velocities, parasitic Kelvin–Helmholtz (KH) modes should be produced, limiting the maximum velocities reached by the VSI modes and eventually resulting in their saturation. Even though KH eddies are, in general, not resolved in hydrodynamical simulations, perturbations in between adjacent VSI-induced flows that could be caused by the Kelvin–Helmholtz instability (KHI) can typically be seen in both 2D and 3D global simulations (e.g., Nelson et al. 2013; Manger & Klahr 2018; Pfeil & Klahr 2021). On the other hand, an alternative saturation mechanism was recently introduced in Cui & Latter (2022), where it was proposed that VSI body modes can saturate by transferring energy to small scales via resonance with pairs of inertial waves produced by nonlinear interaction (parametric instability). This instability has not yet been seen in global VSI simulations, possibly due to the expected ≳30 grid points per wavelength required to resolve it, estimated in that work to correspond to ∼300 cells per scale height (these estimates, however, depend on the wavelength of the VSI modes). It is then clear that very high resolutions are required to get a better understanding of the nonlinear behavior and saturation mechanism of the VSI. Such resolution requirements can be excessive in 3D hydrodynamical simulations, but they can be achieved in axisymmetric simulations at the cost of neglecting non-axisymmetric phenomena. Given that the driving mechanism of the VSI is axisymmetric, so are approximately VSI modes in 3D (e.g., Flock et al. 2020; Pfeil & Klahr 2021; Barraza-Alfaro et al. 2021), and thus useful information on the nonlinear evolution of the instability can be extracted from high-resolution 2D simulations. However, saturation in real disks might also be linked to the formation of non-axisymmetric structures, such as small- and large-scale anticyclonic vortices (Richard et al. 2016; Manger & Klahr 2018; Pfeil & Klahr 2021).

In this work we investigated the growth and evolution of secondary instabilities triggered by the VSI in the axisymmetric radiation-hydrodynamical (Rad-HD) simulations of protoplanetary disks in Paper I and analyzed their connection to the saturation of the VSI. We also extended the local stability analysis considered in that work by including variations of the local thermal relaxation timescale in regions where either the dust-gas collisional timescale or the radiative emission by gas molecules become important. We then applied this analysis to make predictions of different stability regions in our disk models considering varying degrees of dust depletion, and we linked our resulting stability maps to current and future observations of protoplanetary disks.

This article is organized as follows. In Sect. 2 we summarize the disk models and numerical methods employed in this work. In Sect. 3 we characterize the secondary instabilities triggered by the VSI in our simulations. In Sect. 4 we extend the stability analysis in Paper I and produce stability maps for our disk models. In Sect. 5 we discuss caveats and further consequences of our results, focusing in particular on possible saturation mechanisms of the VSI and connecting our predicted stability regions with current observations of protoplanetary disks. Finally, in Sect. 6 we summarize our conclusions. Additional calculations are included in the appendices.

2. Disk models

In this work we consider the same two disk models as in Paper I and employ the same numerical scheme therein. We solve the Rad-HD equations by means of the Rad-HD module by Melon Fuksman et al. (2021) implemented in the open-source PLUTO code (version 4.4, Mignone et al. 2007). We use spherical coordinates (r, θ) with axial symmetry, solving as well for the azimuthal ϕ component of all vector fields. Some quantities in this work are instead expressed in cylindrical coordinates (R, z). We model the gravitational potential and heating produced by a T Tauri star of mass M_s = 0.5 M_⊙, effective temperature T_s = 4000 K, and radius R_s = 2.5 R_⊙. The heating rate of the disk gas-dust mixture due to stellar irradiation is computed via ray-tracing using frequency-dependent absorption opacities taken from (Krieger & Wolf 2020, 2022), assuming small dust grains dynamically well coupled to the gas with radii between 5 and 250 nm. The frequency-averaged values of these opacities are used for the computation of absorption and emission of reprocessed radiation.

Initial hydrostatic conditions are obtained by iterating the computation of temperature and density distributions via the radiative transfer method in the code and a numerical solution of hydrostatic equations, as detailed in Melon Fuksman & Klahr (2022). The gas column density is in every case kept as Σ(r) = 600 g cm⁻² (r/au)⁻¹. Once initial conditions are obtained in the domain (r, θ)∈[0.4, 100] au × [π/2 − 0.5, π/2 + 0.5], the Rad-HD equations are solved in the smaller domain (r, θ)∈[4, 7] au × [π/2 − 0.3, π/2 + 0.3].

All presented simulations are labeled in Table 1 together with their corresponding parameters. For the computation of absorption opacities for both stellar irradiation and radiative transfer of reprocessed infrared photons, we assume a constant dust-to-gas mass ratio f_dg of small (< 0.25 μm) grains. We consider two f_dg values, representing a nominal case with f_dg = 10⁻³ and a dust-depleted case with f_dg = 10⁻⁴. Assuming a total dust-to-gas mass ratio including all grain sizes of , these values correspond respectively to 10% and 1% of the total dust mass being in the small-grain population, which for the grain size distribution dn ∼ a^−3.5 da assumed for the opacity computation correspond to maximum grain sizes of 19 μm and 1.8 mm, respectively. The chosen disk region is either vertically thick or thin and has different stability regions depending on the dust content: while in the nominal case the entire domain is VSI-unstable, in the dust-depleted case the disk is stable at the middle layer located below the irradiation surfaces, defined as the regions of unity optical depth for photons emitted at the star. As detailed in Paper I, this results from the different space-dependent cooling timescales in the domain for varying dust content (see also Sect. 4).

Radiative transfer of reprocessed photons is modeled by evolving the frequency-integrated radiation energy density and flux employing the M1 closure by Levermore (1984). To increase the largest possible time step, the disparity between the radiation and HD timescales is reduced by employing an artificially small value of the speed of light to advect the radiation fields. In this work we take , which, as verified in Paper I via analytical and numerical testing, does not introduce unphysical phenomena in our simulations. Together with the employment of a reduced radial domain compared to the disk size, this allows us to achieve the high resolutions needed to investigate secondary instabilities. We employ a maximum resolution of N_r × N_θ = 1920 × 2048, which corresponds to approximately 200 cells per pressure scale height H (see Table 1) with an aspect ratio of 1 : 1. We limit numerical diffusion by employing third-order WENO spatial reconstruction (Yamaleev & Carpenter 2009; Mignone 2014) and HLLC-type Riemann solvers for both radiation (Melon Fuksman & Mignone 2019) and matter fields (Toro 2009). Dirichlet boundary conditions are applied by fixing all fields to the hydrostatic initial conditions at the ghost zones. Unless otherwise stated, all results presented in this work correspond to our highest-resolution runs.

3. Secondary instabilities

3.1. Overview

We begin by summarizing the main ideas presented in this section, which comprise our main results regarding secondary instabilities resulting from the evolution of the axisymmetric VSI. In its linear stage, the VSI transports gas in nearly vertical directions. Since the unstable directions of motion cross surfaces of constant specific angular momentum j_z = R²Ω, this leads to angular momentum redistribution. As shown in Paper I, the nonlinear development of the VSI results in the formation of nearly uniform-j_z bands, in which the initial vertical shear is significantly reduced. This can be seen in Figs. 1 and 2, showing velocity and j_z distributions after VSI saturation.

Fig. 1. Vertical velocity normalized by the local sound speed (top) and specific angular momentum jz = R2Ω (bottom) in the upper region of our highest-resolution simulations at t = 234 T5.5. White dotted lines indicate the z/r = 0.125 and 0.25 slices displayed in Fig. 2.

Fig. 2. Specific angular momentum at fixed z/r values in Fig. 1. The curves have have been vertically shifted proportionally to z/r for better visualization.

As shown in these figures and schematized in Fig. 3, the gas inside the constant-j_z bands rotates clockwise above the midplane (in a reference frame in which the orbital velocity enters the meridional plane) and counterclockwise below the midplane, with minimum speed ad the center of the bands and maximum speed at the edges. Due to the abrupt reversal of the vertical velocity at the interfaces between adjacent bands, significant meridional shear is created at such regions, eventually triggering the KHI. While this instability typically involves a balance of the destabilizing effect of shear and the stabilizing effect of buoyancy (e.g., in atmospheric physics), the latter is suppressed in our considered system by the fast cooling required for the VSI to operate. Instead, stabilization results from the radial angular momentum stratification, and a KH-stability criterion can be written in terms of a Richardson number in which the Brunt–Vaisälä frequency is replaced by the epicyclic frequency (Sect. 3.2). The VSI growth phase ends as soon as the KH eddies comprise a significant portion of the kinetic energy, supporting the hypothesis that the KHI is the main mechanism behind VSI saturation in disks (Sect. 3.3).

Fig. 3. Schematic view of the KH-unstable zones in a VSI-unstable axisymmetric disk. Gas rotates in the meridional plane inside of approximately constant-jz zones in the sense of rotation favored by the baroclinic torque, producing high-shear regions that become KH-unstable. The transfer of kinetic energy to KH eddies acts as an effective friction between the vertical flows and contributes to the saturation of the VSI.

Various aspects of both the large-scale structure of the vertical flows in the nonlinear stage of the VSI and the amplification of small-scale eddies in our simulations can be explained in terms of the generation of vorticity in rotating baroclinic flows. Locally, the time evolution of the azimuthal vorticity ω_ϕ = (∇×v)_ϕ depends on a combination of three terms, denoted τ_a, τ_b, and τ_c (Sect. 3.4). The terms τ_b and τ_c account for baroclinic and centrifugal torques, respectively, while τ_a accounts for vorticity advection. Once the VSI starts growing, τ_b starts overcoming as constant-j_z regions start forming, contributing to the clockwise rotation of the large-scale flows at z > 0 and counterclockwise rotation at z < 0. Although the term τ_a is generally dominant over the other two, the acceleration or deceleration of vortices advected with the flow only depends on the balance of τ_b and τ_c. Thus, in constant-j_z regions, τ_b is the dominant term contributing to vortex acceleration, which results in an acceleration of vortices with sgnω_ϕ = sgnτ_b and, vice versa, a deceleration of vortices with sgnω_ϕ = −sgnτ_b. Even though the KHI mainly produces vortices rotating in the direction disfavored by the baroclinic torque, its nonlinear evolution also results, via τ_a, in the formation of counter-rotating vortices in regions with sgnω_ϕ = sgnτ_b. These vortices can then be accelerated by τ_b in the constant-j_z bands. Since sgnτ_b = sgnz, this leads to the formation of small-scale vortices rotating clockwise at z > 0 and counterclockwise at z < 0, which are accelerated up to significant fractions of the local sound speed (Sects. 3.4 and 3.5). Further acceleration of these vortices is likely produced via advection from regions adjacent to the shear layers with sgnω_ϕ = sgnτ_b. The size of these vortices is limited by the width of the constant-j_z regions. In axisymmetric simulations, in which such regions can grow in size unimpeded by non-axisymmetric instabilities, the baroclinically amplified vortices can eventually grow up to large sizes (∼1 au) dominating the gas dynamics, as detailed in Sect. 3.6.

3.2. Kelvin-Helmholtz instability

We begin by analyzing the triggering of the KHI in between uniform-j_z regions. As shown in Figs. 1, 3, and 4, the rotation pattern inside of the approximately uniform-j_z bands causes downward-moving gas parcels at z > 0 to experience significant shear on their larger-r edges, while below the midplane shear is maximum on the smaller-r edges. For f_dg = 10⁻³, in which case the VSI flows cross the midplane, this produces an intercalated pattern of KH-unstable regions, each of them occupying half of the vertical domain (see also Fig. 5). The same occurs for f_dg = 10⁻⁴ in the upper VSI-unstable layers. The resulting eddies are clearly highlighted by ω_ϕ, which also shows their sense of rotation. We show this quantity in Fig. 4 next to the normalized vertical velocity in a region close to the midplane for run dg3c4_2048.

Fig. 4. Velocity and vorticity distributions in a region close to the midplane in run dg3c4_2048 after 300 orbits. Red and blue regions in the vorticity map correspond to clockwise and counterclockwise rotation, respectively. KH eddies on the lower hemisphere rotate clockwise, while their nonlinear evolution produces counterclockwise-rotating eddies that are amplified by the baroclinic torque, visible in that region as dark blue spots.

To operate, the KHI requires that the kinetic energy available to displace a fluid element across a shear layer overcomes the counteracting work of restoring forces. This balance is generally quantified in terms of the Richardson number (e.g., Chandrasekhar 1961), defined as the ratio between these energies. While typically (e.g. in shear flows perpendicular to gravity in the atmosphere or the ocean) buoyancy provides the restoring force, this effect is in our case suppressed by fast cooling. Conversely, the increase of j_z with radius contributes to the stabilization of radial displacements. In light of these considerations, we define the radial Richardson number as

(1)

where is the local epicyclic frequency. To obtain this expression, we have simply replaced the squared restoring frequency in the numerator, typically given by the buoyancy frequency (e.g., Chandrasekhar 1961), with the squared frequency ω² of radially oscillating parcels in stratified disks, given in Goldreich & Schubert (1967), Urpin (2003), Klahr (in prep.). While in the adiabatic limit , where N_R is the radial Brunt–Vaisälä frequency, ω² tends to for fast thermal relaxation compared to epicyclic motion (t_cool ≪ Ω⁻¹, as shown in Sect. 4.1), leading to Eq. (1). In the context of rotating stars, an equivalent expression is obtained in Maeder et al. (2013) in the limits of instant cooling and uniform molecular weight.

A necessary condition for KH-stability is Ri < 1/4 (Chandrasekhar 1961; Maeder et al. 2013). As soon as VSI modes are formed, this condition is only violated in the high-shear regions in between vertical flows (shown in Fig. 5 in the nonlinear stage). At those locations, KH eddies are visible within a few orbits after this condition is met, as shown in Fig. 6. An alternative instability criterion based on a linear analysis of perturbations to the VSI modes was derived in Latter & Papaloizou (2018), requiring that kV/Ω ≳ 1.38 to trigger the KHI, where k is the wavenumber of the VSI modes and V their maximum velocity amplitude. This criterion is also consistent with our results, as we only see signs of KHI growth in regions where kV/Ω is at least of order unity (typically at least ∼2), which we have verified using the wavelength distributions estimated in Paper I.

Fig. 5. Same as Fig. 4, showing instead (∂Rvz)2 and the squared epicyclic frequency , both normalized by the squared rotation angular velocity.

As the VSI starts forming bands of approximately uniform j_z, where, consequently, κ_R ≪ Ω (Fig. 5), the meridional shear in the nonlinear state causes the condition Ri < 1/4 to be satisfied in the entire regions where the VSI operates (top right panel in Fig. 6). Despite this, we do not observe a similar spontaneous formation of eddies inside of such bands as at their interfaces. This can possibly be explained by a lower KHI growth rate in those regions, which we expect to be proportional to the shear ∂_Rv_z. For instance, the growth rate for a uniform fluid with a discontinuous velocity transition Δv is k|Δv|/2 (e.g., Chandrasekhar 1961), where k is the component of the wavenumber parallel to the discontinuity interface.

Fig. 6. Onset of the KHI in run dg3c4_2048 in a region at the disk upper layers, where the instability first occurs. Top: snapshots at different times (in units of the orbital period at 5.5 au) showing the deviation of the radial Richardson number (Eq. (1)) from its critical value of 1/4 required to trigger the KHI. Black regions correspond to Ri > 1/4. Bottom: same snapshot showing the azimuthal vorticity, highlighting the formation of eddies. The rightmost panels show the same quantities > 200 orbits after VSI saturation.

The VSI not only increases the shear in the meridional plane. As evidenced by the κ_R distribution in Fig. 5, jumps in the specific angular momentum in between vertical flows produce regions of large azimuthal shear (∂_Rv_ϕ). If the same happens in 3D, as it appears to be the case in the simulations by Richard et al. (2016), Manger & Klahr (2018), and Flock et al. (2020; see Paper I), we can expect the KHI to produce azimuthal eddies in the same regions. Since j_z increases with radius, such eddies should have prograde rotation with respect to the disk. Therefore, long-lived anticyclonic vortices such as those in Manger & Klahr (2018) and Pfeil & Klahr (2021), if formed via KHI, can only result from secondary eddies produced by the nonlinear evolution of this instability in adjacent regions with a negative vorticity perturbation, similarly to the formation of counter-rotating eddies in the meridional plane (Sect. 3.4).

3.3. Energy spectra

Further information linking the VSI saturation process with the onset of the KHI can be obtained by inspecting the evolution of the kinetic energy spectrum in our highest-resolution runs. As detailed in Appendix B, we achieve this by operating on the Fourier-transformed meridional momentum and velocity distributions in an upper disk region exhibiting VSI growth for both f_dg values. Energy spectra between 6 and 37 orbital times are shown in Fig. 7 as a function of kH, where k = ||(k₁, k₂)|| is the wavenumber norm and H is the average local pressure scale height in the selected region. This value is larger than at the midplane, and therefore so is the resolution relative to H, with H/Δr ∼ 263 − 288 for f_dg = 10⁻³ and 296 − 324 for f_dg = 10⁻⁴. The curves have been averaged every 3 snapshots (∼2.3 orbits) to reduce noise. To highlight the energy redistribution over time, we also show in that figure the same spectra normalized by their integrated values.

Fig. 7. Energy spectra. Top: Values in code units at different times computed as described in Appendix B. Bottom: Normalized values computed at the same times. Some relevant slopes mentioned in Sects. 3.3 and 5.3 are shown for comparison.

As the VSI modes grow in amplitude, energy increases over time for all k. At the linear stage, the VSI primary modes produce an energy peak around kH ∼ 20 − 30, which, after sharply decreasing at kH ∼ 30, decays approximately as k⁻⁵ up to kH ≳ 100. During that phase, energy at the largest wavenumbers exists due to the functional form of the developed VSI flows. In particular, as the VSI grows, the resulting flows develop sharp variations in the vicinity of the interfaces between the forming constant-j_z bands.

The beginning of the saturation phase, occurring at ∼15 − 25 orbits depending on f_dg, coincides with a deviation from the initial spectrum’s shape. At that time, a “hump” is formed at kH ∼ 50 − 200 due to the formation of KH eddies in that size range. The spectra converge to a new shape at saturation, with a short plateau following the main VSI peak up to kH ∼ 50. This plateau is followed by a shallower decrease with k than in the linear growth stage, with E(k)∼k^{−1.7 ± 0.2} at kH ∼ 50 − 100 for f_dg = 10⁻³ and E(k)∼k^{−2.7 ± 0.1} at kH ∼ 50 − 200 for f_dg = 10⁻⁴. For higher wavenumbers, the spectra decay as k^−ς, with ς ≳ 7.

The described transition in the shape of the spectra over time has the consequence that the proportion of the total kinetic energy contained in large scales is smaller in the saturated state than in the linear growth phase. More precisely, taking kH = 35 as a limit between small and large scales (this value is in between the main peak produced by the VSI modes and the region dominated by KH eddies), the energy in large scales decreases between ∼95% at the growth stage and 50%−60% after saturation. This shows that saturation occurs as soon as KH eddies become energetically important, suggesting that some kinetic energy injected by the VSI driving mechanism is diverted into creating KH modes, eventually halting the growth of the VSI.

The obtained steepening of the spectra at kH > 100 − 200 can have several possible explanations. First of all, the width of the shear layers triggering the KHI imposes a wavenumber range in which KH linear modes can grow (Chandrasekhar 1961). Thus, we should expect the energy spectrum to decrease faster above the maximum wavenumber allowing for KHI growth, with a slope determined by a combination of the energy spectrum of the high-k end of the VSI-induced structures and possibly a direct enstrophy cascade (a discussion on the energy transport over scales is given in Sect. 5.3). The precise limits of this range typically depend on the fluid’s velocity, the density stratification, and the Richardson number, but for this analysis it is enough to consider that its upper limit is some number on the order of 1/d, where d is the typical width of the shear layers (Chandrasekhar 1961). In the region considered for the computation of the energy spectra, we typically have d ∼ 0.02 − 0.03H, corresponding to wavenumbers kH ∼ 200 − 300, coinciding with the region where we see the slope transition. In this regard, we must mention that these layers are typically resolved by ∼10 cells, and thus it is possible that numerical diffusion is enlarging their typical width with respect to its zero-diffusion value, artificially shifting the slope transition to lower wavenumbers.

On the other hand, numerical diffusion can also contribute to the steepening of the energy spectrum at high k. This occurs in the 2D homogeneous turbulence simulation made by Seligman & Laughlin (2017) with PLUTO, also using third-order Runge-Kutta scheme and WENO3 reconstruction. In that test, numerical diffusion leads to a steeper spectrum than the predicted E(k)∼k⁻³ (for a direct enstrophy cascade) up to scales 10 times larger than the grid’s resolution. However, unlike in our simulations, no energy is continuously injected into the system, and thus the spectrum does not reach a steady state over time. Moreover, we note that our spectra steepen exhibit further steepening at kH ∼ 600 − 800 (about 3 cells per wavelength), which could indicate that the dissipation range is instead closer to this region. Yet, the relevance of numerical diffusion on the high-end of our VSI spectra can only be assessed via a proper comparison with 2D forced turbulence tests injecting energy at k ≲ 2π/d, as we plan to do in a future study.

3.4. Vorticity generation

In our highest-resolution simulations (N_θ = 1024 and 2048), we observe a formation of small meridional vortices in the VSI-active regions that survive up to 1 − 10 orbital timescales and get accelerated up to hundreds of m s⁻¹, with Mach numbers up to 0.4 for f_dg = 10⁻³ and 0.2 for f_dg = 10⁻⁴ (see Figs. 8 and 10). These vortices, which are typically responsible for the maximum velocities reached in the domain (Fig. 5 in Paper I), rotate clockwise at z > 0 and counterclockwise at z < 0. Similar meridional vortices also form in the 2D isothermal and β-cooling simulations by Klahr et al. (2023). As argued by these authors, the observed acceleration is not produced by the convective overstability (COS, Klahr & Hubbard 2014; Lyra 2014), given that the cooling timescale is much smaller than Ω in the VSI-active regions where the accelerated vortices form (Sect. 4.1). Further evidence for this is the fact that these vortices appear even in isothermal simulations, that is, in absence of buoyancy, as we also verified in isothermal versions of our setup at the same resolution (not shown here). On the other hand, the driving mechanism cannot be the VSI, since the same phenomenon occurs if the gravitational potential is modified in such a way that the disk is still baroclinic but has no vertical shear (Klahr et al. 2023). On the contrary, the vortices thrive in bands of approximately uniform specific angular momentum, where they are instead accelerated by the meridional baroclinic torque. This is a novel instability mechanism in the context of protoplanetary disks, analogous to the sea breeze phenomenon in atmospheric physics (see, e.g., Holton & Hakim 2013).

Fig. 8. Meridional velocity norm (top), vorticity (middle), and pressure perturbation (bottom) distributions in the upper region of our highest-resolution simulations at t = 234 T5.5. The log(p/p0) values are multiplied by 5 for fdg = 10−4 for better visualization.

We can understand how the driving mechanism works by examining the time evolution of the integrated vorticity inside of closed streamlines in the (R, z)-plane advected with the fluid. We begin by considering a closed streamline C in the meridional plane defined in such a way that it is at all locations tangent to the (R, z)-projected velocity, (C is not strictly a streamline in the sense that its tangent vector is not v but v_p). In Appendix A, it is shown that the circulation around C, namely, the clockwise line integral of v along C, follows the evolution equation

(2)

where S is the surface enclosed by C and the time derivative takes into account the evolution of C as the fluid evolves. Equivalently, using Stokes’ theorem, the circulation can be replaced as a surface integral as

(3)

where ω_ϕ = (∇ × v)_ϕ. In these expressions, we have defined the baroclinic and centrifugal terms, τ_b and τ_c respectively, as

(4)

where is the analog to quantifying the vertical squared specific angular momentum gradient. The first and second terms in the right-hand side of Eqs. (2) and (3) quantify respectively the work per unit density along close curves of the centrifugal force and the expansion work. An illustrative example, taken from Kundu & Cohen (2002), explains why the second term leads to the generation of vorticity: a mass element in a baroclinic fluid experiences a nonzero net torque due to the misalignment of the center of gravity with respect to the line that passes through the center of buoyancy in the direction of the net pressure force. In a similar way, a nonzero produces a nonzero net torque on meridional vortices due to the mismatch of the centrifugal acceleration RΩ² when the gas moves outward and inward in R.

A similar expression quantifying the vorticity evolution at fixed locations can be obtained by taking the curl of the velocity evolution equation (see Appendix A), namely

(5)

In other words, ω_ϕ can vary either by transport throughout the domain, determined by the advection term

(6)

or by local generation, determined by the torque terms τ_b and τ_c. If τ_b = τ_c = 0, Eq. (5) becomes a conservation equation for the total integrated ω_ϕ.

In the initial hydrostatic state of our simulations, we have τ_a = 0 and τ_c + τ_b = 0, and thus no vorticity is initially generated (the latter expression is the thermal wind equation relating the vertical shear to the disk’s baroclinicity; see Eq. (5) in Paper I). As shown in Fig. 9, this situation changes as the VSI develops, until eventually |τ_a|≫|τ_b|,|τ_c| in most of the domain. However, τ_a can only transport existing vorticity without creating nor destroying it. Thus, as long as vortices are advected but not destroyed via nonlinear interaction with the flow, their acceleration or deceleration depends solely on the balance of τ_b and τ_c in their enclosed surfaces (Eq. (2)). As soon as uniform-j_z bands start to form, the term significantly decreases in such regions until it becomes negligible compared to τ_b. Therefore, the acceleration of closed streamlines transported into those regions is entirely determined by the baroclinic torque. This means that vortices in these bands with sgnω_ϕ = sgnτ_b are accelerated and, vice versa, vortices with sgnω_ϕ = −sgnτ_b are decelerated. Since sgnτ_b = sgnz, the amplified vortices in our simulations rotate clockwise for z > 0 and counterclockwise for z < 0.

Fig. 9. Advection, baroclinic, and centrifugal terms (τa, τb, and τc, respectively) determining the vorticity evolution (Eqs. (3) and (5)), computed at the same snapshot as Fig. 8.

The KHI is first triggered in shear layers with ω_ϕ of the opposite sign as τ_b (see, e.g., Figs. 4 and 6) due to the spatial distribution of the VSI flows (Fig. 3). This can be easily seen by taking v_z ≫ v_R, leading to ω_ϕ ≈ −∂_Rv_z, and resulting in sgnω_ϕ = −sgnz at the shear layers. Therefore, the initially produced KH eddies can only be decelerated by τ_b if they are transported into the constant-j_z zones. However, the KHI also produces via advection vortices of opposite ω_ϕ as the initially produced KH eddies, some of which can be seen in Figs. 4, 6, and 8. Since the term τ_a can advect vorticity but not create it, such vortices can only be created by accumulating flows which already have sgnω_ϕ = sgnz. This process can be seen in Fig. 6, where it is shown that the shear layer (with sgnω_ϕ = −sgnz) is adjacent to a region of opposite ω_ϕ-sign (see also Fig. 9). Once the KHI is triggered¹, such regions form vortices with sgnω_ϕ = sgnz. Eddies created in this way can then be accelerated in the constant-j_z zones by τ_b.

In the vorticity distribution for f_dg = 10⁻⁴ (Fig. 8), it can be seen that the amplified eddies are always adjacent to the zones with sgnω_ϕ = sgnz located at the smaller-r side of the shear layers (it is harder to see this for f_dg = 10⁻³, in which case the flow structure is less orderly). It is then possible that these vortices are being constantly fed gas with their same vorticity sign, possibly contributing to their acceleration. This seems to be a minor effect in the irradiated upper layers, where only the vortices with sgnω_ϕ = sgnz prevail, while vortices with sgnω_ϕ = −sgnz are inhibited by τ_b. Close to the midplane, instead, the baroclinic torquing is weaker (Fig. 9), and eddies rotate with sgnω_ϕ = −sgnz as determined by the KHI driving, while only a few lower-amplitude counter-rotating eddies survive.

In summary, the observed vortices with sgnω_ϕ = sgnτ_b accelerated in the uniform-j_z bands are seeded at the band interfaces via τ_a by the nonlinear evolution of the KHI². The amplified vortices are rapidly destroyed when advected into the shear layers (possibly with some contribution of τ_c), and thus their size is limited by the width of the uniform-j_z regions.

The prevailing baroclinic torque at the constant-j_z bands may also explain why the VSI-induced large-scale flows also respect a rotation pattern with sgnω_ϕ = sgnτ_b. In its linear stage, the VSI transports gas across constant-j_z surfaces, as determined by the unstable directions of motion predicted by linear theory (James & Kahn 1970; Knobloch & Spruit 1982; Klahr, in prep.). Since ∇j_z decreases with height and increases with radius, this transport reduces the j_z gradient between gas parcels moving away from the midplane and their adjacent outer parcels moving toward it, which starts forming constant-j_z bands. Upward- and downward-directed flows are eventually connected at the upper disk layers, either due to the angular momentum excess of the parcels moving away from the midplane or due to the rotation produced by the baroclinic torque once |τ_b|> |τ_c|. This rotation further enhances the angular momentum exchange between the connected regions, until bands of uniform j_z are formed, inside of which |τ_b|≫|τ_c|. Thus, the baroclinic acceleration may explain why the saturated VSI modes survive despite the lowered angular momentum gradient: even though the VSI driving is initially of centrifugal origin, it is the baroclinic torque what drives this rotation pattern once the centrifugal acceleration gradient vanishes.

3.5. A closer look at meridional vortices

A gas parcel embedded in a vortex travels through regions of varying temperature, and thus experiences successive cooling and heating cycles. If the thermal relaxation time was longer than the eddy overturn time, the vortex temperature would become approximately uniform, which would in turn get rid of the local baroclinicity, since τ_b ∝ ∇ρ × ∇T (Eq. (5) in Paper I). Therefore, in order to work, this process requires shorter cooling times than the rotation period of the vortices. This condition is satisfied in the entire domain, as the vortex rotation angular frequency, Ω_v, is typically on the order of the orbital frequency, while the cooling time above the irradiation surface, where vortices are formed, is 10⁻⁴ − 10⁻² Ω⁻¹ for f_dg = 10⁻³ and 10 times larger for f_dg = 10⁻⁴ (Sect. 4.1). On the other hand, despite these short cooling times, we can see some heating and cooling throughout the vortices, as seen in the temperature perturbations in Fig. 10. This occurs due to the expansion and contraction, evidenced by the sign of the work term p∇⋅v in that figure, caused by thermal relaxation as gas is transported to and from regions of increasing entropy, respectively (the direction of the entropy gradient is indicated on the same figure). The integral of this term in a region occupied by an amplified vortex is positive, indicating the conversion of internal to kinetic energy expected from the acceleration mechanism.

Fig. 10. Vortex highlighted in the upper left panel of Fig. 8. Clockwise from the top left: meridional velocity norm, expansion work term (code units), temperature perturbation, and sum of the baroclinic and centrifugal torque terms in Eq. (3) normalized by Ω2. The direction of the meridional velocity and the entropy gradient are indicated with arrows in the left and right top panels, respectively.

An explanation for the fact that Ω_v = 𝒪(Ω) can be given by noting that bands of approximately constant angular momentum are radially compressed by the imbalance between gravity and centrifugal acceleration (see, e.g., Papaloizou & Pringle 1984). Neglecting 𝒪(z/R)² terms, the sum of these accelerations is , which is negative for large radii and positive for small radii. If one such band is radially kept in place by the pressure of its adjacent regions, this quantity must be be positive at the inner edge and negative at the outer edge of the band. In that case, neglecting meridional velocity fluctuations, hydrostatic balance in the radial direction determines that , which means that the pressure must reach a maximum at some radius R = R₀ inside of the band. This explains the positive pressure perturbations inside of the constant-j_z bands in Fig. 8. Assuming that relative temperature variations are much smaller than relative pressure variations, which can be done as long as the cooling time is much smaller than the dynamical time (t_cool ≪ Ω⁻¹), we can solve the resulting equation for p by writing with constant , which for small radial departures ΔR ≪ R₀ results in , where Ω₀ = Ω(R₀). On the other hand, the vortices produce pressure minima that compensate for the outward centrifugal force due to their own rotation, as it can be seen in Fig. 8. For a vortex of radius ∼ΔR, taking ΔR as the half-width of a VSI mode, this balance results in a relative pressure decrease . Thus, vortices can only be accelerated as long as they can produce deeper pressure minima, but this process is limited by the pressure increase caused by the continuously enforced radial contraction of the constant-j_z bands. We can then expect the saturation of the vortex amplification to be reached when the pressure perturbation caused by both processes is of a similar order of magnitude, which leads to , and so Ω_v ∼ Ω.

In their saturated state, the amplified eddies alter the fluid state in such a way that τ_c is no longer zero but of the same order of magnitude as τ_b, which is also altered with respect to the initial state. Each of these terms individually have opposite signs in different regions of the saturated amplified vortices, and so does their sum, which controls the vortex’s acceleration (Eq. (3)). Despite the generally complex form of this quantity (see Fig. 10), we typically see that the negative regions approximately balance the positive regions when integrated inside of closed streamlines. This behavior is expected, as such a balance is required in order to halt the acceleration process.

3.6. Long-time evolution

In axisymmetric VSI simulations, bands of approximately constant j_z can grow up to large widths (∼1 au in our simulations) for long enough times. This occurs due to the mixing of angular momentum enforced to occur in the meridional plane, given that these structures are unaffected by non-axisymmetric modes that should prevent them from growing to such extents in 3D. Since the size of the baroclinically amplified eddies is only limited by the width of the bands, these eddies are allowed to occupy larger regions as the bands grow. For f_dg = 10⁻³, in which case the VSI is active in the entire domain, the j_z gradient is eventually largely reduced in most of the domain (except at the few remaining interfaces between bands) after thousands of orbits. In that state, most of the kinetic energy is contained in the amplified vortices, as shown in Fig. 11. Although we do see vortex mergers, potentially indicating an inverse energy cascade resulting from the enforced axisymmetry, it is unclear whether it is this phenomenon or the baroclinic driving that is the main mechanism behind the formation of the large eddies (see Sect. 5.3).

Fig. 11. Kinetic energy density due to meridional motion (; left) and specific angular momentum (right) after 1360 orbits in run dg3c4_1024.

Even if the resolution is not large enough to initially resolve the amplified eddies, as it occurs for N_θ = 512, these are eventually resolved for sufficiently long times due to the increase of their typical size resulting from the growth of the bands. The transition into this regime occurs at earlier times for increasing resolution, as evidenced by the time evolution of the kinetic energy (Fig. 12). This shows that the mixing of angular momentum in between constant-j_z regions is more efficient when smaller scales are resolved, and in particular when the shear in between these regions is larger, which may result from the angular momentum mixing between adjacent bands produced by the KH eddies, evidenced, for instance, in Fig. 1. Moreover, this also shows that the observed transition is not caused by numerical diffusion of angular momentum, as that would lead to the opposite trend.

Fig. 12. Time evolution of the average meridional kinetic energy in runs dg3c4_512 and dg3c4_1024.

A transition into an eddy-dominated regime is also seen in isothermal versions of our setup (not shown here) and in the isothermal runs by Klahr et al. (2023). This shows that, even though the increased optical depth of the large eddies with respect to the VSI modes and the subsequent increase of the local cooling timescales in the Rad-HD simulations could potentially trigger the COS, that instability is not what produces the transition into that regime. The required high resolutions and long times to obtain this transition (in our case ∼3000 and 800 orbits for 50 and 100 cells per scale height, respectively, with our employed numerical scheme) explain why this phenomenon has not been seen in other works.

The mentioned expected differences with non-axisymmetric simulations suggest that this transition is merely a consequence of the enforced axisymmetry. Thus, the long-term nonlinear evolution of the VSI can only be studied in high-resolution 3D simulations (see Sects. 5.3 and 5.4).

4. Stability regions

4.1. Influence of collisional timescales and gas emissivity

In Paper I, we analyzed the stability of our Rad-HD simulations in terms of the local criterion by Urpin (2003), which predicts instability as long as

(7)

where t_rel is the local thermal relaxation time, which in the Rad-HD simulations coincides with the radiative cooling time t_cool, while N_z is the vertical Brunt–Väisälä frequency (e.g., Rüdiger et al. 2002). The localized suppression of the VSI in the middle layer of our dust-depleted disk is then explained by the fact that this criterion is not met in that entire region. On the other hand, the vertically global stability criterion by Lin & Youdin (2015), predicting instability for

(8)

where Γ is the heat capacity ratio, , and Ω_K is the Keplerian rotation frequency, correctly predicts the stability properties of the middle layer but does not contemplate cases in which the disk surface layers, where the assumptions of vertically uniform temperature and cooling timescale break down, have different stability conditions than the midplane.

The stability analysis in Paper I relies on the assumption that gas and dust particles are in local thermal equilibrium. However, this is not the case at the surface layers, whose low densities make collisions between dust and gas particles infrequent enough that these species may have different temperatures. Considering this, we can still estimate t_rel in those regions by means of a linear analysis of small perturbations of separate dust and gas temperatures. We follow this approach in Appendix C in the same way as Barranco et al. (2018), with the difference that our largest considered dust grains are micron-sized, whereas in that work cm- to meter-sized grains were assumed. Despite this difference, the functional form of the obtained timescale is identical as in that work³:

(9)

where the radiative cooling timescale, t_cool, is computed as in Paper I, while the collisional characteristic timescale for gas thermal relaxation, t_g, is obtained as

(10)

where ρ_gr is the bulk density of the dust grains, is the total local dust density, a_max and a_min are respectively the maximum and minimum grain radii, and is the mean thermal speed of gas particles. For an assumed total dust-to-gas ratio , our assumed dust-to-gas ratios for small < 0.25 μm grains, f_dg = 10⁻³ and 10⁻⁴, correspond respectively to a_max = 19 μm and 1.8 mm. Even if dust settling imposes a lower cutoff to the maximum dust size than a_max, the value of t_g is unchanged, as shown in Appendix C. We also show in the same appendix that our obtained expression for t_rel is unchanged if the dust temperature depends slowly on the grain size.

The thermal relaxation times in our simulations, computed using Eq. (9), are shown in Figs. 13 and 14 (labeled as ), where it can be seen that the increase of the thermal relaxation time at the surface layers can inhibit the growth of VSI modes, as it occurs, for example, in the hydrodynamical simulations by Pfeil & Klahr (2021). For f_dg = 10⁻³, the VSI can still operate up to a few scale heights above the midplane, while the instability should be nearly completely suppressed for f_dg = 10⁻⁴.

Fig. 13. Thermal relaxation times normalized by the critical cooling time tcrit (top) and the orbital frequency (bottom) at r = 5.5 au for fdg = 10−3, assuming only radiative cooling (trel = tcool), including the effect of dust collisions (), and including gas emission with defined in Eq. (13). Also shown are the corresponding values for trel = tcrit. The gas opacities reported in Freedman et al. (2008), Malygin et al. (2014, 2017) for solar abundances correspond to .

Fig. 14. Same as Fig. 13 but for fdg = 10−4.

This situation changes if cooling by molecular emission is not negligible. This effect should be most important at the surface layers, where the typical temperatures of up to ∼200 − 250 K in our domain may increase the amount of possible free-free and bound-free emission (and absorption) mechanisms (e.g., Freedman et al. 2008; Malygin et al. 2014) with respect to the cooler midplane. If optically thin radiative emission by gas molecules is considered, it can be shown (Appendix C) that the linear relaxation timescale is modified as

(11)

where t_rel, thin is the optically thin limit of this expression, given by

(12)

where and are respectively the Planck and Rosseland mean free paths due to both gas and dust photon absorption (see Appendix C), while t_cool, thin is the optically thin limit of t_cool (Eq. (17) in Paper I) and t_r, g is the optically thin cooling time due to gas emission. Assuming that the gas emission rate is determined by Kirchhoff’s law, we compute the latter timescale as , where is the gas specific heat capacity at constant volume, is the Planck-averaged gas absorption opacity, and . Gas opacities can in general span several orders of magnitude (∼10⁻⁶ − 1 cm² g⁻¹) depending on chemical composition and molecular abundances (Freedman et al. 2008; Helling & Lucas 2009; Dullemond & Monnier 2010; Malygin et al. 2014). If we assume, as we have done throughout this work, that the temperatures of the gas and the reprocessed radiation are similar, then we should also expect the gas opacity to drop at low temperatures (≲75 K in Freedman et al. 2008; Malygin et al. 2014, assuming solar composition) due to the lack of absorption lines at low frequencies. Taking these considerations into account, we adopt a simplified gas opacity law of the form

(13)

where Θ is the Heaviside function and the constant parameterizes the uncertainty in the gas composition. For the computation of the combined dust-gas Rosseland opacity, we only consider the contribution of dust absorption, assuming that our densities are low enough that absorption lines can be disregarded in this calculation (see, e.g., Freedman et al. 2008).

The obtained relaxation times at 5.5 au for different values are shown in Figs. 13 and 14. These figures show that the VSI can still operate at the surface layers provided gas opacities are large enough ( in our models). For , the gas opacities are large enough that even the middle layer is VSI-unstable in our dust-depleted case. However, we can only expect the temperature and density distributions to be roughly unchanged for varying opacity if is smaller than both and , where and are respectively the Planck and Rosseland means of the dust absorption opacity. Using these criteria, we estimate the computed t_rel/t_crit values in our disk models to be reliable at the middle layers for and 10⁻² for f_dg = 10⁻³ and 10⁻⁴, respectively, while at the upper layers the temperatures the temperature and density distributions should be unchanged for . Thus, our conclusion regarding the possible instability of the upper layers remains the same.

4.2. Global disk stability

In this section we extend our stability analysis to the full domain of the hydrostatic disk models used as initial conditions (see Paper I), for which an estimation of the dominant VSI wavenumbers is required. To compute t_cool, we assume the maximum dominant wavenumbers in our simulations to be approximately k ∼ 70/H (Fig. C.3 in Paper I) in all disk regions. For larger radii than in our simulation domain, this does not introduce any error, given that our obtained VSI modes are radially optically thin (Paper I) and therefore cool at the scale-independent optically thin rate given by Eq. (16) of that work, while for larger radii the optical depth along any wavelength only decreases with radius. For smaller radii, on the other hand, we can instead expect the dominant wavelengths in our simulations to become optically thick close to the midplane at some minimum radius, in which case VSI modes of such wavelengths cool via diffusion. If t_rel > t_crit at some location for our considered k value, that just means that only modes with larger wavenumbers can be unstable at that point. Considering that the strength of the instability decreases in our simulations for increasing average wavenumber (Paper I), it is also likely that this is also translated into either a partial or total suppression of the VSI.

The resulting t_rel/t_crit values computed in this way are shown in Fig. 15 between 0.6 and 90 au (this is far enough from the radial boundaries to avoid boundary artifacts, as we used zero-gradient conditions at those boundaries to construct the hydrostatic models). Interpreting these distributions as approximate VSI stability maps, we obtain that the unstable surface regions extend up to a few tens of au for high enough gas opacity (). For f_dg = 10⁻³, the midplane is unstable between ∼2 au (VSI modes cool via diffusion for smaller radii), and ∼100 au, after which densities are low enough that the optically thin cooling timescale is above t_crit. For f_dg = 10⁻⁴, the entire midplane is stable unless , in which case an instability region is formed between 1 and 5 au due to the decrease of the optically thin cooling time with .

Fig. 15. VSI-stability maps showing the ratio of the relaxation time to the critical cooling time defined in the local stability criterion (Eq. (7)). From left to right, trel/tcrit values are computed assuming only radiative cooling (trel = tcool), including dust collisions (), and including as well molecular emission with (Eq. (13)). The top and bottom row correspond to the nominal and dust-depleted cases, respectively. Also shown in the leftmost panels is the location of the domain used in our Rad-HD simulations (white dashed box). The gas opacities reported in Freedman et al. (2008), Malygin et al. (2014, 2017) for solar abundances correspond to .

For comparison with Malygin et al. (2017), we also computed stability maps based on a local evaluation of Eq. (8) as done in that work, which we show in Fig. 16. Both criteria coincide in that region for the following reasons: (I) the assumptions of the global criterion are met below the irradiation surface, and (II) the vertically global critical cooling time (right-hand side of Eq. (8)) can be obtained by evaluating t_crit at z = ΓH/2, which is much smaller than the vertical extent of the middle layer. We also see a general agreement on the predicted unstable regions at the surface layers, with some differences for small . We defer a thorough comparison of local stability criteria with hydrodynamical simulations to future works.

Fig. 16. Same as Fig. 15 but locally computing tcrit as in the right-hand side of the global stability criterion (Eq. (8)).

Differences with the stability maps at the midplane in Flock et al. (2017, 2020) and Pfeil & Klahr (2019, 2021), in particular the fact that an inner stable region and an outer stable region are obtained in the latter, are explained by the different assumed dust masses, opacity laws, and dominant wavelengths in those works and this one. The radiative cooling time is in general diffusion-dominated in inner regions and thus it increases with increasing opacity (e.g., Lin & Youdin 2015), which explains the inner stable region in Pfeil & Klahr (2021) and this work. Conversely, far enough from the star that perturbations are optically thin, t_cool grows for increasing radius as the photon mean free path increases (see Paper I), eventually reaching the condition t_cool > t_crit, as shown in Fig. 15 (see also Malygin et al. 2017; Dullemond et al. 2022). The fact that the cooling times in Flock et al. (2017, 2020) are estimated to be diffusion-dominated between 20 and 100 au for f_dg = 10⁻³ despite the similar parameters in that model and ours is explained by the approximation made in those works that the dominant VSI wavelength is ∼H.

The most striking difference with other stability maps is the presence of unstable regions at the surface layers for reasonably high gas opacities, where we should expect a localized VSI activity as seen in our dust-depleted simulations. These regions appear as a consequence of the reduction of the thermal relaxation times above the irradiation surface induced by our temperature and density stratification, and thus they do not appear in stability maps of vertically isothermal disks (Malygin et al. 2017) nor in the models by Pfeil & Klahr (2019), in which the temperature is computed via an effective 1D treatment of vertical radiative diffusion. Since gas radiative emission can only reduce the relaxation timescale for large enough temperatures (≳100 K is a rough estimate based on Malygin et al. 2017), the radial extent of these unstable regions depends not only on the chemical composition of the gas but also on the dust opacity at the disk upper layers and the temperature of the central star, occupying larger radii for higher temperatures as evidenced by Eq. (4) in Paper I. Future studies can test these results in hydrodynamical simulations by employing methods decoupling the dust and gas temperatures, such as that in Muley et al. (2023).

5. Discussion

5.1. KHI in other works

Although clearly identifying the KH eddies in our vorticity distributions requires resolutions of at least N_θ ≥ 1024, we also see similar velocity perturbations to the VSI modes as in Figs. 1 and 4 for lower resolutions (Fig. A.2 in Paper I). This indicates that the KHI is also triggered in our lower-resolution runs (N_θ ≥ 256, H/Δr ≥ 25) even though eddies are not well resolved. Similar signs of developing KHI can be seen in other works on the VSI showing velocity distributions resulting from both 2D and 3D simulations (e.g., Nelson et al. 2013; Manger & Klahr 2018; Pfeil & Klahr 2021), which indicates that this is a rather general phenomenon. This also suggests that the KHI is the instability mechanism behind the meridional vortices in the axisymmetric isothermal simulations by Flores-Rivera et al. (2020).

5.2. Saturation mechanism

The reduction of the energy proportion in large scales once the KHI is triggered (Sect. 3.3) supports the hypothesis, proposed by Latter & Papaloizou (2018), that the growth of the KHI plays an important role in the saturation of the VSI. In Paper I, it can be seen that the perturbations produced by the KHI to the VSI modes become more prominent for increasing resolution (Fig. A.2 in Paper I), as the shear between adjacent bands increases and the KH eddies are better resolved. Moreover, larger resolutions result in smaller average VSI wavelengths (Fig. C.3 in Paper I), and therefore also in a larger density of KH-unstable interfaces. All of these effects may contribute to the obtained reduction of the saturated average kinetic energy and Reynolds stresses with increasing resolution shown in Paper I.

The distribution of the gas velocity in relation to the location of the KH-unstable zones provides a rough picture of this saturation mechanism: once the VSI modes reach sufficiently high velocities to trigger the KHI, the loss of energy to KH eddies in the unstable shear layers acts as an effective friction between the vertically moving flows (see Fig. 3), limiting the maximum kinetic energy that these can gain. Since in our simulations only a very small fraction of the energy is in the highest wavenumbers (Fig. 7), we expect numerical dissipation to play a minor role in the total energy dissipation as long as a direct energy cascade does not occur (see Sect. 5.3). In real 3D disks, viscous dissipation could become important if the KHI triggers a turbulent cascade leading to dissipation down to molecular viscous lengthscales. Instead, in our simulations, it is likely that energy is radiatively dissipated at the medium to large scales of the VSI flows and KH eddies. In particular, some dissipation may be produced by the baroclinic deceleration of the primary KH eddies, which grow with sgnω_ϕ = −τ_b.

Another possible saturation mechanism is the parametric instability considered in Cui & Latter (2022), whose longest wavelengths are expected to require ∼30 grid points per VSI wavelength to be resolved. In our simulations, the wavelength of the VSI modes is well resolved with λ/Δr ∼ 15 − 75, and in particular ∼75 in Fig. 4. In that figure we do indeed see small-scale perturbations inside of the VSI modes, but we do not see significant amounts of kinetic energy being transferred to them (see also Fig. 8). Rather, the wavelength of the VSI modes grows over time due to the mixing of angular momentum between approximately uniform-j_z bands (Paper I and Sect. 3.6), and so even if the ripples inside of the VSI modes are produced via parametric instability, that process is not able to destroy such modes. On the other hand, it is unclear whether the observed ripples are instead sound waves emitted at the KH-unstable regions. Thus, we cannot conclude that a parametric instability occurs in our simulations due to the current lack of unequivocal signs of its occurrence in its nonlinear stage, but even if it does, its dynamical effect in our simulations appears to be marginal.

5.3. Energy transport between scales

The transfer of energy from VSI to KHI modes does not imply the occurrence of an energy flux from large to small scales as it occurs in 3D turbulent cascades. Instead, as first studied by Kraichnan (1967), 2D flows typically exhibit a split cascade transporting energy toward large scales and enstrophy toward small scales (some experimental verifications of such a cascade can be found in Bruneau & Kellay 2005; Kelley & Ouellette 2011). Thus, it is possible that a similar phenomenon occurs in our simulations. In this picture, energy is continuously injected in a broad range of wavelengths due to the geometry of the large-scale VSI flows (as we have verified) and transported toward small k starting from the largest k allowing for KHI growth (Sect. 3.3). However, it must be noted that our simulations are not strictly two-dimensional. In particular the mechanisms of vortex stretching (Appendix A) and strain self-amplification⁴, often linked with the development of direct energy cascades in 3D (see, e.g., Carbone & Bragg 2020; Johnson 2021), are absent in 2D flows but not in axisymmetric flows. Moreover, as reviewed by Alexakis & Biferale (2018; see also the discussion in Sengupta & Umurhan 2023), both direct and inverse energy cascades are possible in 3D (non-axisymmetric) turbulent flows when the effects of rotation, stratification, and even compressibility become important.

In our simulations, a hint possibly indicating transport of energy toward lower wavenumbers is given by the observation that vortices merge and grow in size, likely contributing to the formation of large-scale meridional vortices at long times (Sect. 3.6). Also resulting in the accumulation of energy in increasingly large scales is the growth over time of the uniform-j_z bands, likely via in-plane angular momentum mixing between adjacent bands. If an inverse energy cascade occurs, then energy should be radiatively dissipated at the medium and large scales. In this scenario, it is possible that the slopes measured for 5 ≲ kH < 100 − 200 result from a combination of an inverse cascade and the geometry imposed by the VSI modes. In this regard, we must clarify that, despite the spectra decay close to k^−5/3 for f_dg = 10⁻³ as one would expect in fully developed homogeneous and isotropic turbulence (even in 2D), the flows in our simulations depart from these conditions, even in the k-range dominated by KHI modes. At such scales, anisotropy can still be observed even for the velocity components perpendicular to the VSI flows, which track mostly KHI modes. Moreover, the spectra deviate from that functional form after several tens of orbits, as the constant-j_z regions grow. This leads to a steepening over time, with E(k) eventually oscillating in that wavenumber range around ∼k^{−2.4 ± 0.2} for f_dg = 10⁻³ and ∼k^{−3 ± 0.2} for f_dg = 10⁻⁴.

It is unclear whether the accumulation of energy at increasingly large scales is necessarily the result of an inverse cascade. Moreover, neither the observation of sporadic vortex mergers and the growth of uniform-j_z regions, nor the evaluation of the slopes of the energy spectra, are sufficient to estimate the influence of rotation, gravity, weak compressibility, and baroclinicity on the transport of energy and enstrophy between scales. Some information on this can likely be obtained by computing the energy and enstrophy fluxes resulting from each of these phenomena, as done in the multiple analyses in Alexakis & Biferale (2018) for incompressible flows. We defer such an analysis to a future work on the topic.

5.4. Caveats of axisymmetric models

Once the KHI is triggered in our simulations, amplified eddies never cease to form and dissipate, and thus they are able to continuously limit the maximum kinetic energy of the vertical flows even hundreds of orbits after saturation. Unfortunately, our simulations are insufficient to make predictions on the disk’s state after thousands of orbits, given that the artificially enhanced growth of uniform-j_z bands due to vertical angular momentum mixing in 2D eventually leads in some cases to a complete breakdown of the VSI flows and domination of large-scale azimuthal vortices, as shown in Sect. 3.6 (see also Klahr et al. 2023). If the growth of these vortices is enhanced via an inverse energy cascade, this process should likely be inhibited in 3D simulations, provided that a direct energy cascade occurs in that case. With enough azimuthal resolution, we expect the formation of wide uniform-j_z bands (Paper I and Papaloizou & Pringle 1984) to be inhibited by the growth of non-axisymmetric modes, which should disallow the growth of large-scale vortices at long times.

On top of this, it is even unclear whether small-scale eddies can be baroclinically amplified and survive for several orbits if axisymmetry is not imposed. In principle, all of the ingredients needed for the formation of meridional vortices could still be present in that case, given that the VSI driving mechanism is axisymmetric. As a result, the large-scale VSI-induced flows maintain some degree of axisymmetry in some 3D simulations (e.g., Richard et al. 2016; Flock et al. 2020; Barraza-Alfaro et al. 2021) and are likely prone to KHI, as the velocity profiles in Manger & Klahr (2018; Fig. 10) seem to indicate. Axisymmetry, however, could be artificially enhanced by underresolving the ϕ-direction, and thus it is in general unclear to which extent non-axisymmetric modes can erase the initial axisymmetry of the VSI modes. Therefore, the question remains of whether bands of reduced ∇j_z are formed in 3D, and even whether meridional vortices are destroyed by non-axisymmetric modes. As discussed in Paper I and later in this section, such regions appear to form in the 3D HD simulations by Richard et al. (2016), Manger & Klahr (2018), and Flock et al. (2020; see also Fig. A.1 in Paper I), but asserting the universality of this phenomenon requires further careful investigation.

Determining whether a baroclinic amplification of KH eddies can occur in non-axisymmetric disks has the additional difficulty that cell aspect ratios of ∼1 should be achieved to avoid purely axisymmetric effects, together with the high resolutions of at least ∼100 cells per scale height needed to resolve the vortices, which is computationally rather expensive. Interestingly, radial resolutions of H/ΔR ∼ 100 − 150 are reached in the β-cooling 3D disk simulations by Richard et al. (2016) for an aspect ratio of H/R ∼ 0.05 similar to those in this work (see Paper I), albeit with a steeper temperature profile (T ∼ R⁻¹ vs. T ∼ R^−1/2 in this work) and with lower azimuthal and meridional resolutions than in the radial direction (H/(RΔϕ)≈H/(RΔθ)≈20). In those simulations, the VSI initially grows axisymmetrically, forming bands of reduced vertical vorticity ω_z with respect to the initial Keplerian state. These bands are separated by shear layers in which ω_z abruptly increases. When the VSI flows reach sufficient amplitude, non-axisymmetric structures grow, forming azimuthal vortices in the regions of lowered ω_z. This phenomenon can possibly be attributed to the Rossby Wave Instability, as suggested by the authors. The shear layers exhibit some ripples in the meridional plane (e.g., Fig. 13 in that work), which are possibly produced by the KHI.

For all parameters considered in Richard et al. (2016), the reduction of ω_z in between the shear layers with respect to its Keplerian value appears lower than the value Δω_z/Ω = −1/2 that would result from the formation of uniform-j_z zones (for sufficiently small ∂_ϕv_R). This may be a consequence of the non-axisymmetric transport of angular momentum, which should tend to counteract the formation of axisymmetric uniform-j_z regions by the VSI (Paper I and Papaloizou & Pringle 1984). However, the fact that Δω_z < 0 is still indicative that regions of lowered-∇j_z are formed (before non-axisymmetric structures are formed, ω_z ∝ ∂_Rj_z). While we cannot determine from these results whether such regions still maintain some degree of axisymmetry above the midplane, a reduction of ∇j_z would mean that τ_b should overcome τ_c in those locations. It is then possible that a baroclinic amplification of eddies occurs in those simulations, provided this is not restricted by the employed vertical resolution. However, it is challenging to identify this phenomenon in the data presented in that work (Fig. 13) without resorting to 3D vortex detection methods. It is also unclear whether coherent regions of reduced-∇j_z can in some cases be destroyed by non-axisymmetric modes, as it may be the case near the midplane in Fig. 8 of Richard et al. (2016).

Despite the mentioned numerical challenges, the long-term evolution of the VSI and the eventual baroclinic amplification of meridional eddies are problems worth studying in the future. Firstly, besides potentially intervening in the VSI saturation, amplified eddies may have an impact on dust evolution by enhancing the stirring of small dust grains, increasing their average collisional velocities, and even concentrating them in cases where Ω_v does not exceed ∼Ω (Klahr & Henning 1997), for instance during their acceleration stage. Furthermore, understanding the formation and evolution of meridional vortices may cast some light on the birth of azimuthal anticyclonic vortices produced in 3D VSI simulations (Richard et al. 2016; Manger & Klahr 2018; Flock et al. 2020; Pfeil & Klahr 2021) and potentially observed in images of protoplanetary disks (e.g., van der Marel et al. 2016; de Boer et al. 2021; Marr & Dong 2022), which can boost planet formation by acting as dust traps. If the KHI has a role in the initial formation of such vortices (Latter & Papaloizou 2018), it is likely that these are affected by baroclinic torques in their early stages. Moreover, it is also possible that some vortices form as a result of the breakup of reduced-∇j_z bands by non-axisymmetric modes.

5.5. Stability regions and observational consequences

The fact that our disk models are VSI-stable close to the midplane from a minimum radius is consistent with current dust continuum observations of protoplanetary disks according to the argument by some authors (Flock et al. 2017; Dullemond et al. 2022) that a suppression of the VSI, at least in the outer disk, is required to explain the thin distributions of millimeter-sized dust seen in ALMA dust continuum images of protoplanetary disks, for instance in HD 163296 at r ≈ 100 au (Doi & Kataoka 2021) and in Oph163131 up to ∼170 au (Villenave et al. 2022) (see also Pinte et al. 2016). In particular, a suppression of the VSI at large radii could explain why the dust scale height estimated in Doi & Kataoka (2021) is close to H at 68 au but much smaller at 100 au. The stable regions in the outer parts of our disk models grow in size as the density of small grains is reduced, which is why it has been proposed that these can naturally occur as a consequence of small-dust depletion due to coagulation (Fukuhara et al. 2021; Dullemond et al. 2022; Pfeil et al. 2023). To this picture, our local analysis adds the possibility that the VSI may still operate in surface layers even at distances from the star where it is suppressed at the midplane.

Considering the works by Stoll & Kley (2016) and Flock et al. (2017), we expect a VSI activity confined to regions above the irradiation surface to produce significant vertical stirring of small dust grains, and therefore to limit vertical settling in such regions. We conjecture that this phenomenon can contribute to the vertically extended appearance of edge-on disks at optical and near-infrared wavelengths (Villenave et al. 2020) even at radii where turbulence appears weak or nonexistent at the disk midplane (Villenave et al. 2022), as we expect the vertical stirring of dust grains to push the location of the τ = 1 surface for photon scattering at μm wavelengths away from the midplane. This scenario can be tested in numerical experiments by computing the distributions of different dust populations in a dust fluid approach (e.g., Krapp et al. 2022) and using these to produce synthetic dust continuum and scattered light images. It must be noted, however, that the distributions of small grains in the disk upper layers may also depend on nonideal MHD effects and magnetized winds (Riols & Lesur 2018; Booth & Clarke 2021; Hutchison & Clarke 2021).

Depending on whether the τ = 1 surfaces of ¹²CO lines are located inside the unstable surface layers or not (if these occur), it might be possible to detect signs of velocity structures produced by the VSI in such layers in ALMA CO kinematics observations, as proposed by Barraza-Alfaro et al. (2021). In the same work, it is concluded that the nonthermal spectral broadening of molecular lines produced by the VSI should be negligible compared to current ALMA capabilities, which is consistent with current measured upper bounds (e.g., Teague et al. 2018; Flaherty et al. 2020), and which also applies to an instability localized at surface layers.

5.6. Additional physics affecting the disk stability

Some physical effects have been left out of our stability analysis. On the one hand, the disk middle layer can be stabilized for large enough vertically integrated dust-to-gas mass ratios (Σ_dust/Σ_gas ≳ qH/R) due to the additional buoyancy caused by dust back-reaction (Lin 2019). On the other hand, the strength of the VSI at the surface layers should depend on the large-scale magnetic field and on the degree of ionization in such locations, which are in general poorly constrained (Lesur et al. 2023). In particular, recent studies on the interplay of the VSI and magnetic phenomena indicate that strong magnetic fields well coupled to the gas in highly ionized regions (i.e., in the ideal MHD limit) should suppress the VSI by providing additional buoyancy via magnetic tension (Latter & Papaloizou 2018; Cui & Lin 2021), whereas, in poorly ionized regions, nonideal MHD effects such as the Hall effect and ambipolar and Ohmic diffusion can diminish this effect and favor the growth of the VSI (Cui & Bai 2020, 2022; Cui & Lin 2021; Latter & Kunz 2022). Most of these models assume locally isothermal disks and therefore favor VSI growth. Future nonideal MHD studies including finite thermal relaxation times or radiative transfer, ideally self-consistently estimating the local ionization degree (e.g., Delage et al. 2022), are therefore required in order to better constrain the expected degree of turbulence at the surface layers of protoplanetary disks.

6. Conclusions

In this work we studied the growth and evolution of secondary instabilities parasitic to the VSI modes and their relation with the VSI saturation in the axisymmetric Rad-HD simulations presented in Paper I. We also extended the VSI-stability analysis of Paper I including the effects of dust-gas thermal equilibration and gas molecular emission, and applied it to construct stability maps in our disk models. Our main findings can be summarized as follows:

1. In its linear stage, the VSI produces adjacent upward- and downward-directed flows which are initially disconnected. The transport of angular momentum across constant-j_z surfaces reduces the angular momentum gradient between flows moving away from the midplane and the immediate outward flows moving toward the midplane, causing the baroclinic torque (τ_b) to prevail over the centrifugal torque (τ_c) in those regions. Such adjacent flows eventually become connected in upper disk regions, either due to the angular momentum excess of outward-moving gas parcels or due to baroclinic acceleration.

2. Connected adjacent vertical flows exchange angular momentum in such a way that bands of approximately uniform j_z are formed after a few orbits. The initial τ_c is consequently suppressed inside of these bands, while the initial τ_b is largely unaffected, leading to the production of vorticity in those regions by baroclinic torquing. This may explain why the saturated VSI-induced flows survive despite the reduced angular momentum gradient inside of the bands.

3. While the meridional shear is also reduced inside of the formed uniform-j_z bands, it is maximum at the interface between them, with large enough values to overcome the stabilizing effect of the radial j_z gradient. As a result, the KHI is triggered at those layers. This instability transfers significant kinetic energy from the VSI flows to KH eddies, acting as an effective friction between the vertical flows and limiting the maximum kinetic energy these can gain due to the VSI driving, possibly leading to the saturation of the VSI.

4. Large resolutions favor the KHI by increasing the shear in between uniform-j_z bands and reducing the numerical diffusivity. Together with the increase of the number of KH-unstable interfaces with resolution, this may explain why both the average kinetic energy and the Reynolds stresses after saturation decrease with resolution without converging in our resolution range. We require about 100 cells per scale height to unequivocally identify the KH eddies in vorticity distributions, but the velocity perturbations in the VSI modes show that this instability is triggered in our lower resolution runs as well. Similar signs of developing KHI in other works suggest that this is a rather general phenomenon involved in the VSI saturation.

5. Due to the flow pattern produced by the VSI, the initially produced KH eddies rotate in the direction disfavored by τ_b, namely, with sgnω_ϕ = −sgnτ_b. However, the nonlinear evolution of the KHI at the shear layers also generates counter-rotating eddies in adjacent zones with sgnω_ϕ = sgnτ_b. As a result, some of those vortices get accelerated by τ_b in the uniform-j_z zones, possibly with some contribution of vorticity advection from regions near the shear layers. Vortices accelerated in this way reach Mach numbers up to ∼0.4 for f_dg = 10⁻³ and ∼0.2 for f_dg = 10⁻⁴, corresponding in each case to the maximum velocities in the simulation domain. Conversely, the growth of vortices with sgnω_ϕ = −sgnτ_b is disfavored by τ_b. Besides Klahr et al. (2023), this mechanism has not been seen in other works due to the high resolutions required to resolve the amplified vortices (≳100 cells per scale height with our integration methods).

6. Our stability analysis shows that the infrequent dust-gas collisions at the disk surface layers can locally suppress the VSI, but these regions can still be unstable for reasonably high gas emissivity.

7. An extension of this stability analysis to regions in our disk models left out of the simulation domain shows that, for sufficient depletion of small grains via dust coagulation and sufficiently high molecular emissivity, protoplanetary disks can be VSI-unstable in surface layers up to several tens of au while remaining stable at the midplane. This picture is consistent with current observations of vertically extended edge-on disks at optical and near-infrared wavelengths (e.g., Villenave et al. 2020) appearing as thin layers in ALMA dust continuum images (e.g., Villenave et al. 2022), evidencing a low level of turbulent dust stirring at the midplane. The stability of the surface layers is subject to the gas molecular composition and their local temperature.

If a baroclinic amplification of meridional eddies does indeed occur in protoplanetary disks, it may affect dust coagulation and fragmentation by increasing the collisional velocities of small dust grains, affecting their diffusivity, and potentially concentrating them (Klahr & Henning 1997). On top of this, this process could play a role in the birth of large-scale azimuthal vortices. It is however unclear whether bands of constant specific angular momentum can form in real disks, as these are unstable to non-axisymmetric modes which cannot be captured by our axisymmetric simulations. However, reduced-∇j_z bands seem to be present in a number of 3D simulations (Richard et al. 2016; Manger & Klahr 2018; Flock et al. 2020). Since the formation of such regions is a requirement for the formation of amplified eddies, determining whether this process can occur and, in that case, how it may interact with non-axisymmetric modes requires rather expensive high-resolution 3D simulations. These would make it possible to explore the long-term evolution of the VSI and its induced secondary instabilities without running into artifacts produced in axisymmetric simulations, such as the seemingly unrestricted growth of large uniform-j_z regions.

Further work is also needed to make better predictions on the stability of disk surface layers. In particular, improvements to this work can be made by modeling the distribution and temperature of small dust grains, which determine the collisional timescale and therefore affect the strength of the VSI. Conversely, the balance between vertical stirring and settling of dust grains depends on the level of turbulence, and therefore a self-consistent model is needed. Such models should also consider the likely VSI suppression in strongly ionized regions, and, conversely, the expected prevalence of the VSI in regions where nonideal MHD effects become important (e.g., Latter & Kunz 2022; Cui & Bai 2022).

Acknowledgments

We thank Marcelo Barraza-Alfaro, Remo Burn, Riccardo Franceschi, Dhruv Muley, and Bhargav Vaidya for useful comments and discussions that helped improve this manuscript. The research of J.D.M.F. and H.K. is supported by the German Science Foundation (DFG) under the priority program SPP 1992: “Exoplanet Diversity” under contract KL 1469/16-1/2. We thank our collaboration partners on this project in Kiel, Sebastian Wolf and Anton Krieger, under contract WO 857/17-1/2, for providing us with the employed tabulated opacity coefficients. M.F. acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No. 757957). All numerical simulations were run on the ISAAC and VERA clusters of the MPIA and the COBRA cluster of the Max Planck Society, all of these hosted at the Max-Planck Computing and Data Facility in Garching (Germany). We thank the anonymous referee for many insightful comments that helped us improve the quality and presentation of this work.

References

Appendix A: Vorticity generation

In this appendix we derive a series of expressions determining the generation and transport of vorticity in axisymmetric disks. We start by considering a closed streamline C in the (R, z)-plane, defined such that its tangent vector is at all locations parallel to , this is, the projection of v onto the (R, z)-plane, where are the cylindrical unit vectors. By defining the circulation of v along this curve as the line integral

(A.1)

and taking into account the variation of C with time, specifically, , we can follow analogous steps as in the common derivation of Kelvin’s circulation theorem (e.g., Kundu & Cohen 2002) to obtain

(A.2)

where we have used for the second step the expression of ∇v_p in cylindrical coordinates for axisymmetric flows and the fact that . Replacing the velocity terms using the Navier-Stokes momentum equation, , and using Stokes’ theorem to convert the right-hand side of this equation into a surface integral, we arrive to the desired expression,

(A.3)

where C = ∂S. The left-hand side of this equation can also be transformed using Stokes’ theorem, yielding the following equation for the vorticity (ω = ∇ × v):

(A.4)

Finally, taking (i.e., choosing clockwise circulation for the computation of Γ_C), Eqs. (A.3) and (A.4) become Eqs. (2) and (3).

An equation quantifying the local evolution of the ϕ component of the vorticity can also be obtained by taking the curl of the momentum equation, leading to

(A.5)

The left-hand side of this expression can be simplified by making use of the identity , after which the resulting term ∇ × (v⋅∇v) = ∇ × (ω×v) can be converted into ∇ ⋅ (ωv)−ω ⋅ ∇v, resulting in

(A.6)

This expression is valid in general for 3D flows, as we have not yet assumed axisymmetry. Here it is noteworthy that, unlike in strictly 2D flows, in axisymmetric flows the vortex stretching term ω ⋅ ∇v is not necessarily zero provided v_ϕ ≠ 0. Equation (5) is finally obtained as the ϕ component of this equation after zeroing all ϕ derivatives.

Appendix B: Spectral analysis

In Section 3.3, we show energy spectra obtained by operating on Fourier-transformed velocity and momentum distributions. Since computing discrete Fourier transforms requires employing an evenly spaced grid, we take a small portion of the (r, θ) domain of approximately equal r- and z-extents, L, and rescale it by converting it into a uniform grid defined in a new domain {(x₁, x₂)∈[0, L]×[0, L]}. We do so by defining each distribution V in the new grid as V(x₁(i),x₂(j)) = V(r(i),θ(N₂ − j)), with i = 1, …, N₁ and j = 1, …, N₂, where N₁ × N₂ is the grid resolution of the chosen domain region. To enforce the periodicity of the transformed distributions, we then perform a fourfold replication of our data on an extended grid defined in the domain S = {(x₁, x₂)∈[−L, L]×[−L, L]} with resolution 2N₁ × 2N₂. We achieve this by assigning each considered distribution V to a single quadrant of the new grid and reflecting it across the x₁ and x₂ axes.

To compute energy spectra, we follow a procedure similar to that in Alexakis & Biferale (2018) with some modifications that allow us to consider nonuniform ρ distributions (a similar approach is followed, for example, in Sengupta & Umurhan 2023). We define the Fourier transform of a periodic function F defined on S as

(B.1)

where k is the wavenumber and ⟨ ⋅ ⟩ denotes average over S, in such a way that

(B.2)

With these definitions, it is straightforward to show that the volume average of the meridional kinetic energy, , can be written as

(B.3)

where ^* denotes complex conjugation, while m_p = ρv_p denotes the momentum density in the (R, z)-plane (see Appendix A). Defining a wavenumber norm step Δk, we can then use this expression to approximate the kinetic energy spectrum as

(B.4)

in such a way that E(k)Δk is the kinetic energy within a shell in k-space of radius k and width Δk.

We take the v_p and m_p distributions from the region r ∈ [5.25, 5.75] au and θ ∈ [θ_min + δ, θ_min + δ + Δθ], where Δθ is such that L = 0.5 au, which is approximately equal to the local pressure scale height H = c_sΩ⁻¹ at that location. We choose this upper layer region since it exhibits VSI growth for both considered dust-to-gas ratios, and we choose δ = 0.02 to consider a region separate from the θ-boundary. The domain is large enough that it contains between 5 and 6 radial VSI wavelengths, and small enough that the maximum relative stretching of cells when replacing the logarithmic r-grid with the uniform x1-grid is below 5%. Since in this domain ρ decreases with z down to a factor of 1/20, the obtained spectra are mostly representative of the kinetic energy distribution in the largest-density region.

Appendix C: Thermal relaxation timescale

C.1. Dust-gas collisional timescales

In general, the relaxation timescale t_rel for gas temperature perturbations differs from t_cool if the timescale of energy transfer between dust and gas particles is long enough that the dust and gas temperatures, T_d and T_g respectively, can be different during the relaxation. To consider this effect, we obtain expressions for t_rel following an analogous procedure to the analysis in Barranco et al. (2018). As in Paper I (Appendix B), we neglect advection by taking v = 0. We also neglect gas opacities and assume the heating rate of spherical dust grains due to collisions with gas particles to be given, as in Burke & Hollenbach (1983), by

(C.1)

where a is the grain radius, n_g = ρ/(μu) and are the number density and the mean thermal speed of gas particles, respectively, while 𝒜_H is the thermal accommodation coefficient for atomic and molecular hydrogen interacting with graphite and silicate grains, which we henceforth assume to be ∼0.5 (Burke & Hollenbach 1983).

Let us first assume that all dust grains are at the same temperature. Using the grain size distribution n(a)∼a^−3.5 assumed in our employed opacity model, we can compute the total energy exchange rate per unit volume between gas and dust as

(C.2)

with , where a_min and a_max are respectively the minimum and maximum grain effective radii. The size distribution is normalized in such a way that the total dust density (i.e., not the density of small grains ρ_d used for the computation of opacities) can be computed as , where ρ_gr is the bulk density of the dust grains, and thus . The resulting equations for the evolution of the temperature of gas and dust particles are

(C.3)

where and c_d are the gas and dust specific heat capacities at constant volume and constant pressure, respectively. The radiative cooling term G⁰ in this expression is computed by replacing the term a_RT⁴ (Equation 5 in Paper I) with . A linearization of these equations for small departures δT_g and δT_d with respect to the equilibrium temperature T = T_g = T_d leads to the following system:

(C.4)

where , , and , with . In these equations we have dropped the δE_r term arising from the cG⁰ term, since the present analysis is only required in optically thin regions (more generally, for optically thin perturbations) where δE_r ≪ δ(a_RT⁴) (see Appendix B in Paper I). The gas thermal relaxation time can then be obtained in terms of the evolution of δT_g given an initial condition of the form (δT_g, δT_d) = (δT₀, 0). In that case, the solution of Equation (C.4) is well approximated by δT_g = δT₀ e^−t/t_rel, with

(C.5)

where . As shown in Barranco et al. (2018), an approximation to this expression can be obtained under the condition that t_g ≫ t_d, t_r, in which case t_rel is well approximated by

(C.6)

where t_cool, thin is the optically thin limit of the radiative cooling time (Equations 16 and 17 in Paper I). This approximation was based on the disk parameters in that work, where grain sizes in the range [1, 100] cm were considered, whereas in the present work we assume the much smaller a_max = 19 μm and 1.8 mm for f_dg = 10⁻³ and 10⁻⁴, respectively. This is a more reasonable assumption in the regions where collisional times become relevant, where we can expect centimeter-sized grains to be removed by vertical settling, as shown for example in the model by Flock et al. (2020). Yet, with our employed size distribution, the relations t_r, t_d ≪ t_g are still satisfied, the latter of which we have verified taking the values of the specific heats capacities of graphite and amorphous silicates from Butland & Maddison (1973) and Zeller & Pohl (1971), respectively. Therefore, as we have verified in our disk models, the approximation in Eq. (C.6) is still valid with 0.1% and 0.01% accuracy for f_dg = 10⁻³ and f_dg = 10⁻⁴, respectively. Moreover, since in optically thick regions t_g ≪ t_cool, we can replace Equation (C.6) with the following expression, valid in all radiative transport regimes:

(C.7)

where t_cool is computed as in Paper I (Equation 16). The resulting collisional timescale, t_g, is independent on the gas density, as both the heat capacity ρc_g and the energy exchange rate Λ_dg in Equation (C.3) are proportional to ρ.

Away from the midplane, dust settling might impose a lower cutoff to the dust distribution than the assumed a_max satisfying , especially in the dust-depleted case, in which a_max = 1.8 mm. However, this should not affect the value of t_g, which depends on a_max as . Since and a_min ≪ a_max, this results in , and thus t_g is rather insensitive to the precise value of a_max.

C.2. Inclusion of gas radiative emission

The same analysis can be applied if the emission of radiation by gas molecules is taken into account. We assume for this that the gas emits thermally at a rate determined by the Planck-averaged gas absorption opacity (Kirchhoff’s law), and so we include a gas-radiation interaction terms of the same form of G⁰, specifically . This leads to the following system of equations:

(C.8)

which for small perturbations with respect to equilibrium becomes

(C.9)

where , with . In deriving Equation (C.9) we have again assumed δE_r ≪ δ(a_RT⁴), which is valid for optically thin perturbations. A solution of this system for initial T_g perturbations leads to the optically thin thermal relaxation timescale,

(C.10)

where now . This expression can be further approximated by noting that the factor proportional to is much smaller than 1 for low enough , namely, for low enough , and that t_gt_r/t_r, g ≪ 1 for our disk parameters, which leads to

(C.11)

The lowest accuracy of this approximation in our disk models, achieved for our highest considered value of cm² g⁻¹, is 0.1% and 0.01% for f_dg = 10⁻³ and f_dg = 10⁻⁴, respectively. For largely different t_g, t_r, g, and t_cool, thin, Equation (C.11) is equivalent to

(C.12)

which is the limit of the relaxation timescale considered in Malygin et al. (2017) for vanishing gas-gas collisional times (albeit with some differences in the way of computing t_g).

In dense regions where t_g ≪ t_cool, thin, we can safely assume T_g = T_d and denote these temperatures by T. In that case, t_rel tends to the radiative cooling time that can be obtained by repeating the linear analysis in Paper I (Appendix B) including gas opacities. This can be achieved by simply replacing the Planck and Rosseland photon mean free paths due to dust absorption by the corresponding combined mean free paths for dust and gas absorption, namely,

(C.13)

where and are respectively the Planck and Rosseland means of , while and are the frequency-dependent dust and gas absorption opacities per dust and gas mass, respectively (scattering is neglected in this work, but otherwise should be replaced by the total opacity coefficient ). In this way, the same system of equations as in Paper I (Appendix B) is obtained, and the resulting cooling time is

(C.14)

where

(C.15)

is the optically thin limit of this expression valid for , with . With this result, we can write a general expression for the relaxation timescale valid in all radiative transport and collisional regimes, given by

(C.16)

where t_rel, thin is given by Equation (C.11). Equation (C.16) also comprises the cases in which dust-gas collisional times and gas opacities are not included. In the optically thick limit, in which case t_g ≪ t_cool, this expression tends to the optically thick limit of t_cool. For optically thin perturbations for which t_g ≪ t_cool, it tends to t_cool, thin. Finally, in regions where collisional times become important, the gas perturbations are always optically thin, and thus this expression tends to the correct limit t_rel, thin.

C.3. Size-dependent dust temperature

We lastly explore how the thermal relaxation timescale is changed if the dust temperature has a dependence on size T_d(a)∝a^−b arising from the frequency dependence of the grain emissivity (Glassgold et al. 2004, where b ∈ [1/5, 1/3], and Pauly & Garrod 2016, where b = 1/6), in which case the collisional energy exchange rate per unit volume becomes

(C.17)

where ⟨T_d⟩=I[T_d(a)]/I[1] is the average dust temperature and

(C.18)

Taking , this leads to

(C.19)

where for b < 1/2. The radiative cooling term cG⁰ is now computed using the dust average temperature replacing with , where in the last term we dropped 𝒪(b²) terms. Thus, linearizing these equations for small temperature perturbations, we obtain the system

(C.20)

where now , and . The resulting optically thin thermal relaxation time is

(C.21)

which is well approximated by

(C.22)

given that both and are ≪1. Since the size dependence of T_d results from the balance between absorption of starlight and infrared emission, values of B ≠ 1 can only be assumed above the irradiation surface, where t_cool, thin < t_g. Therefore, Equation (C.22) shows that a slow dependence of the dust temperature with the grain size does not significantly alter the thermal relaxation time. The same conclusion is reached if gas radiative emission is included, in which case the optically thin relaxation time is given by Equation (C.11) replacing t_g + t_cool, thin with the right-hand side of Equation (C.22).

All Tables

All Figures

	Fig. 1. Vertical velocity normalized by the local sound speed (top) and specific angular momentum jz = R2Ω (bottom) in the upper region of our highest-resolution simulations at t = 234 T5.5. White dotted lines indicate the z/r = 0.125 and 0.25 slices displayed in Fig. 2.
In the text

	Fig. 2. Specific angular momentum at fixed z/r values in Fig. 1. The curves have have been vertically shifted proportionally to z/r for better visualization.
In the text

Fig. 3. Schematic view of the KH-unstable zones in a VSI-unstable axisymmetric disk. Gas rotates in the meridional plane inside of approximately constant-jz zones in the sense of rotation favored by the baroclinic torque, producing high-shear regions that become KH-unstable. The transfer of kinetic energy to KH eddies acts as an effective friction between the vertical flows and contributes to the saturation of the VSI.

In the text

Fig. 4. Velocity and vorticity distributions in a region close to the midplane in run dg3c4_2048 after 300 orbits. Red and blue regions in the vorticity map correspond to clockwise and counterclockwise rotation, respectively. KH eddies on the lower hemisphere rotate clockwise, while their nonlinear evolution produces counterclockwise-rotating eddies that are amplified by the baroclinic torque, visible in that region as dark blue spots.

In the text

	Fig. 5. Same as Fig. 4, showing instead (∂Rvz)2 and the squared epicyclic frequency , both normalized by the squared rotation angular velocity.
In the text

Fig. 6. Onset of the KHI in run dg3c4_2048 in a region at the disk upper layers, where the instability first occurs. Top: snapshots at different times (in units of the orbital period at 5.5 au) showing the deviation of the radial Richardson number (Eq. (1)) from its critical value of 1/4 required to trigger the KHI. Black regions correspond to Ri > 1/4. Bottom: same snapshot showing the azimuthal vorticity, highlighting the formation of eddies. The rightmost panels show the same quantities > 200 orbits after VSI saturation.

In the text

	Fig. 7. Energy spectra. Top: Values in code units at different times computed as described in Appendix B. Bottom: Normalized values computed at the same times. Some relevant slopes mentioned in Sects. 3.3 and 5.3 are shown for comparison.
In the text

	Fig. 8. Meridional velocity norm (top), vorticity (middle), and pressure perturbation (bottom) distributions in the upper region of our highest-resolution simulations at t = 234 T5.5. The log(p/p0) values are multiplied by 5 for fdg = 10−4 for better visualization.
In the text

	Fig. 9. Advection, baroclinic, and centrifugal terms (τa, τb, and τc, respectively) determining the vorticity evolution (Eqs. (3) and (5)), computed at the same snapshot as Fig. 8.
In the text

Fig. 10. Vortex highlighted in the upper left panel of Fig. 8. Clockwise from the top left: meridional velocity norm, expansion work term (code units), temperature perturbation, and sum of the baroclinic and centrifugal torque terms in Eq. (3) normalized by Ω2. The direction of the meridional velocity and the entropy gradient are indicated with arrows in the left and right top panels, respectively.

In the text

	Fig. 11. Kinetic energy density due to meridional motion (; left) and specific angular momentum (right) after 1360 orbits in run dg3c4_1024.
In the text

	Fig. 12. Time evolution of the average meridional kinetic energy in runs dg3c4_512 and dg3c4_1024.
In the text

Fig. 13. Thermal relaxation times normalized by the critical cooling time tcrit (top) and the orbital frequency (bottom) at r = 5.5 au for fdg = 10−3, assuming only radiative cooling (trel = tcool), including the effect of dust collisions (), and including gas emission with defined in Eq. (13). Also shown are the corresponding values for trel = tcrit. The gas opacities reported in Freedman et al. (2008), Malygin et al. (2014, 2017) for solar abundances correspond to .

In the text

	Fig. 14. Same as Fig. 13 but for fdg = 10−4.
In the text

Fig. 15. VSI-stability maps showing the ratio of the relaxation time to the critical cooling time defined in the local stability criterion (Eq. (7)). From left to right, trel/tcrit values are computed assuming only radiative cooling (trel = tcool), including dust collisions (), and including as well molecular emission with (Eq. (13)). The top and bottom row correspond to the nominal and dust-depleted cases, respectively. Also shown in the leftmost panels is the location of the domain used in our Rad-HD simulations (white dashed box). The gas opacities reported in Freedman et al. (2008), Malygin et al. (2014, 2017) for solar abundances correspond to .

In the text

	Fig. 16. Same as Fig. 15 but locally computing tcrit as in the right-hand side of the global stability criterion (Eq. (8)).
In the text