Issue 
A&A
Volume 682, February 2024



Article Number  A140  
Number of page(s)  20  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/202346555  
Published online  16 February 2024 
Vertical shear instability in twomoment radiationhydrodynamical simulations of irradiated protoplanetary disks
II. Secondary instabilities and stability regions
Max Planck Institute for Astronomy, Königstuhl 17, 69117 Heidelberg, Germany
email: fuksman@mpia.de
Received:
31
March
2023
Accepted:
30
November
2023
Context. The vertical shear instability (VSI) is a hydrodynamical instability predicted to produce turbulence in magnetically inactive regions of protoplanetary disks. The regions in which this instability can occur and the physical phenomena leading to its saturation are a current matter of research.
Aims. We explore the secondary instabilities triggered by the nonlinear evolution of the VSI and their role in its saturation. We also expand on previous investigations on stability regions by considering temperature stratifications enforced by stellar irradiation and radiative cooling, and including the effects of dustgas collisions and molecular line emission.
Methods. We modeled the gasdust mixture in a circumstellar disk around a T Tauri star by means of highresolution axisymmetric radiationhydrodynamical simulations including stellar irradiation with frequencydependent opacities, considering different degrees of depletion of small dust grains.
Results. The flow pattern produced by the interplay of the axisymmetric VSI modes and the baroclinic torque forms bands of nearly uniform specific angular momentum. In the highshear regions in between these bands, the Kelvin–Helmholtz instability (KHI) is triggered. A third instability mechanism, consisting of an amplification of eddies by baroclinic torques, forms meridional vortices with Mach numbers up to ∼0.4. Our stability analysis suggests that protoplanetary disks can be VSIunstable in surface layers up to tens of au for reasonably high gas emissivities.
Conclusions. The significant transfer of kinetic energy to smallscale eddies produced by the KHI and possibly even the baroclinic acceleration of eddies limit the maximum energy of the VSI modes, likely leading to the saturation of the VSI. Depending on the gas molecular composition, the VSI can operate at surface layers even in regions where the midplane is stable. This picture is consistent with current observations of disks showing thin midplane millimetersized dust layers while appearing vertically extended in optical and nearinfrared wavelengths.
Key words: hydrodynamics / instabilities / radiative transfer / methods: numerical / protoplanetary disks
© The Authors 2024
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model.
Open access funding provided by Max Planck Society.
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 VSIinduced 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 largescale 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 BarrazaAlfaro et al. (2021), largescale 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 razorthin millimeterdust 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 VSIstable 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 VSIinduced 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 nonaxisymmetric 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; BarrazaAlfaro et al. 2021), and thus useful information on the nonlinear evolution of the instability can be extracted from highresolution 2D simulations. However, saturation in real disks might also be linked to the formation of nonaxisymmetric structures, such as small and largescale 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 radiationhydrodynamical (RadHD) 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 dustgas 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 RadHD equations by means of the RadHD module by Melon Fuksman et al. (2021) implemented in the opensource 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 gasdust mixture due to stellar irradiation is computed via raytracing using frequencydependent 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 frequencyaveraged 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^{−2} (r/au)^{−1}. Once initial conditions are obtained in the domain (r, θ)∈[0.4, 100] au × [π/2 − 0.5, π/2 + 0.5], the RadHD 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 dusttogas mass ratio f_{dg} of small (< 0.25 μm) grains. We consider two f_{dg} values, representing a nominal case with f_{dg} = 10^{−3} and a dustdepleted case with f_{dg} = 10^{−4}. Assuming a total dusttogas mass ratio including all grain sizes of , these values correspond respectively to 10% and 1% of the total dust mass being in the smallgrain 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 VSIunstable, in the dustdepleted 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 spacedependent cooling timescales in the domain for varying dust content (see also Sect. 4).
Varying parameters of the simulations analyzed in this work.
Radiative transfer of reprocessed photons is modeled by evolving the frequencyintegrated 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 thirdorder WENO spatial reconstruction (Yamaleev & Carpenter 2009; Mignone 2014) and HLLCtype 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 highestresolution 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^{2}Ω, this leads to angular momentum redistribution. As shown in Paper I, the nonlinear development of the VSI results in the formation of nearly uniformj_{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 j_{z} = R^{2}Ω (bottom) in the upper region of our highestresolution simulations at t = 234 T_{5.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 constantj_{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 KHstability 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 KHunstable zones in a VSIunstable axisymmetric disk. Gas rotates in the meridional plane inside of approximately constantj_{z} zones in the sense of rotation favored by the baroclinic torque, producing highshear regions that become KHunstable. 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 largescale structure of the vertical flows in the nonlinear stage of the VSI and the amplification of smallscale 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 constantj_{z} regions start forming, contributing to the clockwise rotation of the largescale 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 constantj_{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 counterrotating vortices in regions with sgnω_{ϕ} = sgnτ_{b}. These vortices can then be accelerated by τ_{b} in the constantj_{z} bands. Since sgnτ_{b} = sgnz, this leads to the formation of smallscale 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 constantj_{z} regions. In axisymmetric simulations, in which such regions can grow in size unimpeded by nonaxisymmetric 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. KelvinHelmholtz instability
We begin by analyzing the triggering of the KHI in between uniformj_{z} regions. As shown in Figs. 1, 3, and 4, the rotation pattern inside of the approximately uniformj_{z} bands causes downwardmoving gas parcels at z > 0 to experience significant shear on their largerr edges, while below the midplane shear is maximum on the smallerr edges. For f_{dg} = 10^{−3}, in which case the VSI flows cross the midplane, this produces an intercalated pattern of KHunstable regions, each of them occupying half of the vertical domain (see also Fig. 5). The same occurs for f_{dg} = 10^{−4} in the upper VSIunstable 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 counterclockwiserotating 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
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 ω^{2} 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, ω^{2} tends to for fast thermal relaxation compared to epicyclic motion (t_{cool} ≪ Ω^{−1}, 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 KHstability is Ri < 1/4 (Chandrasekhar 1961; Maeder et al. 2013). As soon as VSI modes are formed, this condition is only violated in the highshear 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 (∂_{R}v_{z})^{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 ∂_{R}v_{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 (∂_{R}v_{ϕ}). 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, longlived 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 counterrotating 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 highestresolution runs. As detailed in Appendix B, we achieve this by operating on the Fouriertransformed 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_{1}, k_{2}) 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^{−3} and 296 − 324 for f_{dg} = 10^{−4}. 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^{−5} 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 constantj_{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^{−3} and E(k)∼k^{−2.7 ± 0.1} at kH ∼ 50 − 200 for f_{dg} = 10^{−4}. 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 highk end of the VSIinduced 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 zerodiffusion 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 thirdorder RungeKutta scheme and WENO3 reconstruction. In that test, numerical diffusion leads to a steeper spectrum than the predicted E(k)∼k^{−3} (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 highend 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 highestresolution simulations (N_{θ} = 1024 and 2048), we observe a formation of small meridional vortices in the VSIactive regions that survive up to 1 − 10 orbital timescales and get accelerated up to hundreds of m s^{−1}, with Mach numbers up to 0.4 for f_{dg} = 10^{−3} and 0.2 for f_{dg} = 10^{−4} (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 VSIactive 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 highestresolution simulations at t = 234 T_{5.5}. The log(p/p_{0}) values are multiplied by 5 for f_{dg} = 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
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
where ω_{ϕ} = (∇ × v)_{ϕ}. In these expressions, we have defined the baroclinic and centrifugal terms, τ_{b} and τ_{c} respectively, as
where is the analog to quantifying the vertical squared specific angular momentum gradient. The first and second terms in the righthand 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Ω^{2} 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
In other words, ω_{ϕ} can vary either by transport throughout the domain, determined by the advection term
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 uniformj_{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 ω_{ϕ} ≈ −∂_{R}v_{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 constantj_{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^{1}, such regions form vortices with sgnω_{ϕ} = sgnz. Eddies created in this way can then be accelerated in the constantj_{z} zones by τ_{b}.
In the vorticity distribution for f_{dg} = 10^{−4} (Fig. 8), it can be seen that the amplified eddies are always adjacent to the zones with sgnω_{ϕ} = sgnz located at the smallerr side of the shear layers (it is harder to see this for f_{dg} = 10^{−3}, 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 loweramplitude counterrotating eddies survive.
In summary, the observed vortices with sgnω_{ϕ} = sgnτ_{b} accelerated in the uniformj_{z} bands are seeded at the band interfaces via τ_{a} by the nonlinear evolution of the KHI^{2}. 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 uniformj_{z} regions.
The prevailing baroclinic torque at the constantj_{z} bands may also explain why the VSIinduced largescale flows also respect a rotation pattern with sgnω_{ϕ} = sgnτ_{b}. In its linear stage, the VSI transports gas across constantj_{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 constantj_{z} bands. Upward and downwarddirected 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^{−4} − 10^{−2} Ω^{−1} for f_{dg} = 10^{−3} and 10 times larger for f_{dg} = 10^{−4} (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)^{2} 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_{0} inside of the band. This explains the positive pressure perturbations inside of the constantj_{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} ≪ Ω^{−1}), we can solve the resulting equation for p by writing with constant , which for small radial departures ΔR ≪ R_{0} results in , where Ω_{0} = Ω(R_{0}). 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 halfwidth 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 constantj_{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. Longtime 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 nonaxisymmetric 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^{−3}, 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 constantj_{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 eddydominated 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 RadHD 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 nonaxisymmetric simulations suggest that this transition is merely a consequence of the enforced axisymmetry. Thus, the longterm nonlinear evolution of the VSI can only be studied in highresolution 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 RadHD simulations in terms of the local criterion by Urpin (2003), which predicts instability as long as
where t_{rel} is the local thermal relaxation time, which in the RadHD 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 dustdepleted 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
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 micronsized, whereas in that work cm to metersized grains were assumed. Despite this difference, the functional form of the obtained timescale is identical as in that work^{3}:
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
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 dusttogas ratio , our assumed dusttogas ratios for small < 0.25 μm grains, f_{dg} = 10^{−3} and 10^{−4}, 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^{−3}, 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^{−4}.
Fig. 13. Thermal relaxation times normalized by the critical cooling time t_{crit} (top) and the orbital frequency (bottom) at r = 5.5 au for f_{dg} = 10^{−3}, assuming only radiative cooling (t_{rel} = t_{cool}), including the effect of dust collisions (), and including gas emission with defined in Eq. (13). Also shown are the corresponding values for t_{rel} = t_{crit}. The gas opacities reported in Freedman et al. (2008), Malygin et al. (2014, 2017) for solar abundances correspond to . 
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 freefree and boundfree 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
where t_{rel, thin} is the optically thin limit of this expression, given by
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 Planckaveraged gas absorption opacity, and . Gas opacities can in general span several orders of magnitude (∼10^{−6} − 1 cm^{2} g^{−1}) 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
where Θ is the Heaviside function and the constant parameterizes the uncertainty in the gas composition. For the computation of the combined dustgas 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 VSIunstable in our dustdepleted 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^{−2} for f_{dg} = 10^{−3} and 10^{−4}, 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 scaleindependent 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 zerogradient 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^{−3}, 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^{−4}, 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. VSIstability 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, t_{rel}/t_{crit} values are computed assuming only radiative cooling (t_{rel} = t_{cool}), including dust collisions (), and including as well molecular emission with (Eq. (13)). The top and bottom row correspond to the nominal and dustdepleted cases, respectively. Also shown in the leftmost panels is the location of the domain used in our RadHD 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 (righthand 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 t_{crit} as in the righthand 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 diffusiondominated 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 diffusiondominated between 20 and 100 au for f_{dg} = 10^{−3} 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 dustdepleted 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 lowerresolution 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 FloresRivera 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 KHunstable 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 KHunstable 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 smallscale 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 uniformj_{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 KHunstable 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 largescale 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 twodimensional. In particular the mechanisms of vortex stretching (Appendix A) and strain selfamplification^{4}, 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 (nonaxisymmetric) 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 largescale 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 uniformj_{z} bands, likely via inplane 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^{−3} 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 krange 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 constantj_{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^{−3} and ∼k^{−3 ± 0.2} for f_{dg} = 10^{−4}.
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 uniformj_{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 uniformj_{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 largescale 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 uniformj_{z} bands (Paper I and Papaloizou & Pringle 1984) to be inhibited by the growth of nonaxisymmetric modes, which should disallow the growth of largescale vortices at long times.
On top of this, it is even unclear whether smallscale 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 largescale VSIinduced flows maintain some degree of axisymmetry in some 3D simulations (e.g., Richard et al. 2016; Flock et al. 2020; BarrazaAlfaro 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 nonaxisymmetric 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 nonaxisymmetric 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 nonaxisymmetric 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^{−1} 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, nonaxisymmetric 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 uniformj_{z} zones (for sufficiently small ∂_{ϕ}v_{R}). This may be a consequence of the nonaxisymmetric transport of angular momentum, which should tend to counteract the formation of axisymmetric uniformj_{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 nonaxisymmetric structures are formed, ω_{z} ∝ ∂_{R}j_{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 nonaxisymmetric modes, as it may be the case near the midplane in Fig. 8 of Richard et al. (2016).
Despite the mentioned numerical challenges, the longterm 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 nonaxisymmetric modes.
5.5. Stability regions and observational consequences
The fact that our disk models are VSIstable 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 millimetersized 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 smalldust 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 edgeon disks at optical and nearinfrared 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 ^{12}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 BarrazaAlfaro 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 dusttogas mass ratios (Σ_{dust}/Σ_{gas} ≳ qH/R) due to the additional buoyancy caused by dust backreaction (Lin 2019). On the other hand, the strength of the VSI at the surface layers should depend on the largescale 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 selfconsistently 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 RadHD simulations presented in Paper I. We also extended the VSIstability analysis of Paper I including the effects of dustgas 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:

In its linear stage, the VSI produces adjacent upward and downwarddirected flows which are initially disconnected. The transport of angular momentum across constantj_{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 outwardmoving gas parcels or due to baroclinic acceleration.

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 VSIinduced flows survive despite the reduced angular momentum gradient inside of the bands.

While the meridional shear is also reduced inside of the formed uniformj_{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.

Large resolutions favor the KHI by increasing the shear in between uniformj_{z} bands and reducing the numerical diffusivity. Together with the increase of the number of KHunstable 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.

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 counterrotating eddies in adjacent zones with sgnω_{ϕ} = sgnτ_{b}. As a result, some of those vortices get accelerated by τ_{b} in the uniformj_{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^{−3} and ∼0.2 for f_{dg} = 10^{−4}, 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).

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

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 VSIunstable 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 edgeon disks at optical and nearinfrared 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 largescale azimuthal vortices. It is however unclear whether bands of constant specific angular momentum can form in real disks, as these are unstable to nonaxisymmetric 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 nonaxisymmetric modes requires rather expensive highresolution 3D simulations. These would make it possible to explore the longterm 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 uniformj_{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 selfconsistent 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).
Although the KHI is first triggered in the layer with sgnω_{ϕ} = −sgnz, Fig. 6 shows that the adjacent layer with the opposite vorticity should also be KHunstable (typically with a smaller growth rate), which could also facilitate the formation of eddies with sgnω_{ϕ} = sgnz.
Acknowledgments
We thank Marcelo BarrazaAlfaro, 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/161/2. We thank our collaboration partners on this project in Kiel, Sebastian Wolf and Anton Krieger, under contract WO 857/171/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 MaxPlanck 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
 Alexakis, A., & Biferale, L. 2018, Phys. Rep., 767, 1 [CrossRef] [Google Scholar]
 Barranco, J. A., Pei, S., & Marcus, P. S. 2018, ApJ, 869, 127 [NASA ADS] [CrossRef] [Google Scholar]
 BarrazaAlfaro, M., Flock, M., Marino, S., & Pérez, S. 2021, A&A, 653, A113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Booth, R. A., & Clarke, C. J. 2021, MNRAS, 502, 1569 [NASA ADS] [CrossRef] [Google Scholar]
 Bruneau, C. H., & Kellay, H. 2005, Phys. Rev. E, 71, 046305 [NASA ADS] [CrossRef] [Google Scholar]
 Burke, J. R., & Hollenbach, D. J. 1983, ApJ, 265, 223 [NASA ADS] [CrossRef] [Google Scholar]
 Butland, A. T. D., & Maddison, R. J. 1973, J. Nucl. Mater., 49, 45 [NASA ADS] [CrossRef] [Google Scholar]
 Carbone, M., & Bragg, A. D. 2020, J. Fluid Mech., 883, R2 [NASA ADS] [CrossRef] [Google Scholar]
 Chandrasekhar, S. 1961, Hydrodynamic and hydromagnetic stability (Oxford: Clarendon Press) [Google Scholar]
 Cui, C., & Bai, X.N. 2020, ApJ, 891, 30 [NASA ADS] [CrossRef] [Google Scholar]
 Cui, C., & Bai, X.N. 2022, MNRAS, 516, 4660 [NASA ADS] [CrossRef] [Google Scholar]
 Cui, C., & Latter, H. N. 2022, MNRAS, 512, 1639 [NASA ADS] [CrossRef] [Google Scholar]
 Cui, C., & Lin, M.K. 2021, MNRAS, 505, 2983 [NASA ADS] [CrossRef] [Google Scholar]
 de Boer, J., Ginski, C., Chauvin, G., et al. 2021, A&A, 649, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Delage, T. N., Okuzumi, S., Flock, M., Pinilla, P., & Dzyurkevich, N. 2022, A&A, 658, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Doi, K., & Kataoka, A. 2021, ApJ, 912, 164 [NASA ADS] [CrossRef] [Google Scholar]
 Dullemond, C. P., & Monnier, J. D. 2010, ARA&A, 48, 205 [NASA ADS] [CrossRef] [Google Scholar]
 Dullemond, C. P., Ziampras, A., Ostertag, D., & Dominik, C. 2022, A&A, 668, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Flaherty, K., Hughes, A. M., Simon, J. B., et al. 2020, ApJ, 895, 109 [NASA ADS] [CrossRef] [Google Scholar]
 Flock, M., Nelson, R. P., Turner, N. J., et al. 2017, ApJ, 850, 131 [Google Scholar]
 Flock, M., Turner, N. J., Nelson, R. P., et al. 2020, ApJ, 897, 155 [Google Scholar]
 FloresRivera, L., Flock, M., & Nakatani, R. 2020, A&A, 644, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Freedman, R. S., Marley, M. S., & Lodders, K. 2008, ApJS, 174, 504 [NASA ADS] [CrossRef] [Google Scholar]
 Fukuhara, Y., Okuzumi, S., & Ono, T. 2021, ApJ, 914, 132 [NASA ADS] [CrossRef] [Google Scholar]
 Fukuhara, Y., Okuzumi, S., & Ono, T. 2023, PASJ, 75, 233 [NASA ADS] [CrossRef] [Google Scholar]
 Gerbig, K., MurrayClay, R. A., Klahr, H., & Baehr, H. 2020, ApJ, 895, 91 [CrossRef] [Google Scholar]
 Glassgold, A. E., Najita, J., & Igea, J. 2004, ApJ, 615, 972 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P., & Schubert, G. 1967, ApJ, 150, 571 [Google Scholar]
 Hartmann, L., Calvet, N., Gullbring, E., & D’Alessio, P. 1998, ApJ, 495, 385 [Google Scholar]
 Helling, C., & Lucas, W. 2009, MNRAS, 398, 985 [NASA ADS] [CrossRef] [Google Scholar]
 Holton, J. R., & Hakim, G. J. 2013, An Introduction to Dynamic Meteorology (Cambridge: Academic Press) [Google Scholar]
 Hutchison, M. A., & Clarke, C. J. 2021, MNRAS, 501, 1127 [Google Scholar]
 James, H. A., & Kahn, F. D. 1970, A&A, 5, 232 [NASA ADS] [Google Scholar]
 Johansen, A., Oishi, J. S., Mac Low, M.M., et al. 2007, Nature, 448, 1022 [Google Scholar]
 Johansen, A., Blum, J., Tanaka, H., et al. 2014, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond,& T. Henning (Tucson: University of Arizona Press), 547 [Google Scholar]
 Johnson, P. L. 2021, J. Fluid Mech., 922, A3 [NASA ADS] [CrossRef] [Google Scholar]
 Kelley, D. H., & Ouellette, N. T. 2011, Phys. Fluids, 23, 115101 [NASA ADS] [CrossRef] [Google Scholar]
 Klahr, H. H., & Henning, T. 1997, Icarus, 128, 213 [NASA ADS] [CrossRef] [Google Scholar]
 Klahr, H., & Hubbard, A. 2014, ApJ, 788, 21 [Google Scholar]
 Klahr, H., & Schreiber, A. 2020, ApJ, 901, 54 [NASA ADS] [CrossRef] [Google Scholar]
 Klahr, H., Baehr, H., & Melon Fuksman, J. D. 2023, ApJ, submitted [arXiv:2305.08165] [Google Scholar]
 Knobloch, E., & Spruit, H. C. 1982, A&A, 113, 261 [NASA ADS] [Google Scholar]
 Kraichnan, R. H. 1967, Phys. Fluids, 10, 1417 [NASA ADS] [CrossRef] [Google Scholar]
 Krapp, L., Kratter, K. M., & Youdin, A. N. 2022, ApJ, 928, 156 [NASA ADS] [CrossRef] [Google Scholar]
 Krieger, A., & Wolf, S. 2020, A&A, 635, A148 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Krieger, A., & Wolf, S. 2022, A&A, 662, A99 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kundu, P. K., & Cohen, I. M. 2002, Fluid mechanics. Second edition (San Diego: Academic Press) [Google Scholar]
 Latter, H. N., & Kunz, M. W. 2022, MNRAS, 511, 1182 [CrossRef] [Google Scholar]
 Latter, H. N., & Papaloizou, J. 2018, MNRAS, 474, 3110 [NASA ADS] [CrossRef] [Google Scholar]
 Lesur, G., Flock, M., Ercolano, B., et al. 2023, ASP Conf. Ser., 534, 465 [NASA ADS] [Google Scholar]
 Levermore, C. D. 1984, J. Quant. Spectr. Rad. Transf., 31, 149 [NASA ADS] [CrossRef] [Google Scholar]
 Lin, M.K. 2019, MNRAS, 485, 5221 [CrossRef] [Google Scholar]
 Lin, M.K., & Youdin, A. N. 2015, ApJ, 811, 17 [NASA ADS] [CrossRef] [Google Scholar]
 Lyra, W. 2014, ApJ, 789, 77 [Google Scholar]
 Maeder, A., Meynet, G., Lagarde, N., & Charbonnel, C. 2013, A&A, 553, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Malygin, M. G., Kuiper, R., Klahr, H., Dullemond, C. P., & Henning, T. 2014, A&A, 568, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Malygin, M. G., Klahr, H., Semenov, D., Henning, T., & Dullemond, C. P. 2017, A&A, 605, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Manara, C. F., Rosotti, G., Testi, L., et al. 2016, A&A, 591, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Manger, N., & Klahr, H. 2018, MNRAS, 480, 2125 [NASA ADS] [CrossRef] [Google Scholar]
 Marr, M., & Dong, R. 2022, ApJ, 930, 80 [NASA ADS] [CrossRef] [Google Scholar]
 Melon Fuksman, J. D., & Klahr, H. 2022, ApJ, 936, 16 [NASA ADS] [CrossRef] [Google Scholar]
 Melon Fuksman, J. D., & Mignone, A. 2019, ApJS, 242, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Melon Fuksman, J. D., Klahr, H., Flock, M., & Mignone, A. 2021, ApJ, 906, 78 [NASA ADS] [CrossRef] [Google Scholar]
 Melon Fuksman, J. D., Flock, M., & Klahr, H. 2024, A&A, 682, A139 (Paper I) [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mignone, A. 2014, J. Comput. Phys., 270, 784 [Google Scholar]
 Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228 [Google Scholar]
 Muley, D., Melon Fuksman, J. D., & Klahr, H. 2023, A&A, 678, A162 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610 [Google Scholar]
 Ormel, C. W., & Cuzzi, J. N. 2007, A&A, 466, 413 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Papaloizou, J. C. B., & Pringle, J. E. 1984, MNRAS, 208, 721 [NASA ADS] [CrossRef] [Google Scholar]
 Pauly, T., & Garrod, R. T. 2016, ApJ, 817, 146 [Google Scholar]
 Pfeil, T., & Klahr, H. 2019, ApJ, 871, 150 [NASA ADS] [CrossRef] [Google Scholar]
 Pfeil, T., & Klahr, H. 2021, ApJ, 915, 130 [NASA ADS] [CrossRef] [Google Scholar]
 Pfeil, T., Birnstiel, T., & Klahr, H. 2023, ApJ, 959, 121 [NASA ADS] [CrossRef] [Google Scholar]
 Pinte, C., Dent, W. R. F., Ménard, F., et al. 2016, ApJ, 816, 25 [Google Scholar]
 Richard, S., Nelson, R. P., & Umurhan, O. M. 2016, MNRAS, 456, 3571 [NASA ADS] [CrossRef] [Google Scholar]
 Riols, A., & Lesur, G. 2018, A&A, 617, A117 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rüdiger, G., Arlt, R., & Shalybkov, D. 2002, A&A, 391, 781 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Seligman, D., & Laughlin, G. 2017, ApJ, 848, 54 [NASA ADS] [CrossRef] [Google Scholar]
 Sengupta, D., & Umurhan, O. M. 2023, ApJ, 942, 74 [NASA ADS] [CrossRef] [Google Scholar]
 Stoll, M. H. R., & Kley, W. 2016, A&A, 594, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Teague, R., Henning, T., Guilloteau, S., et al. 2018, ApJ, 864, 133 [NASA ADS] [CrossRef] [Google Scholar]
 Toro, E. F. 2009, Riemann Solvers and Numerical Methods for Fluid Dynamics (Berlin, Heidelberg: Springer) [Google Scholar]
 Urpin, V. 2003, A&A, 404, 397 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Urpin, V., & Brandenburg, A. 1998, MNRAS, 294, 399 [Google Scholar]
 van der Marel, N., Cazzoletti, P., Pinilla, P., & Garufi, A. 2016, ApJ, 832, 178 [Google Scholar]
 Villenave, M., Ménard, F., Dent, W. R. F., et al. 2020, A&A, 642, A164 [EDP Sciences] [Google Scholar]
 Villenave, M., Stapelfeldt, K. R., Duchêne, G., et al. 2022, ApJ, 930, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Yamaleev, N. K., & Carpenter, M. H. 2009, J. Comput. Phys., 228, 4248 [NASA ADS] [CrossRef] [Google Scholar]
 Zeller, R. C., & Pohl, R. O. 1971, Phys. Rev. B, 4, 2029 [NASA ADS] [CrossRef] [Google Scholar]
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
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
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 NavierStokes momentum equation, , and using Stokes’ theorem to convert the righthand side of this equation into a surface integral, we arrive to the desired expression,
where C = ∂S. The lefthand side of this equation can also be transformed using Stokes’ theorem, yielding the following equation for the vorticity (ω = ∇ × v):
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
The lefthand 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
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 Fouriertransformed 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 zextents, L, and rescale it by converting it into a uniform grid defined in a new domain {(x_{1}, x_{2})∈[0, L]×[0, L]}. We do so by defining each distribution V in the new grid as V(x_{1}(i),x_{2}(j)) = V(r(i),θ(N_{2} − j)), with i = 1, …, N_{1} and j = 1, …, N_{2}, where N_{1} × N_{2} 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_{1}, x_{2})∈[−L, L]×[−L, L]} with resolution 2N_{1} × 2N_{2}. We achieve this by assigning each considered distribution V to a single quadrant of the new grid and reflecting it across the x_{1} and x_{2} 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
where k is the wavenumber and ⟨ ⋅ ⟩ denotes average over S, in such a way that
With these definitions, it is straightforward to show that the volume average of the meridional kinetic energy, , can be written as
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
in such a way that E(k)Δk is the kinetic energy within a shell in kspace 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}Ω^{−1} at that location. We choose this upper layer region since it exhibits VSI growth for both considered dusttogas 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 rgrid with the uniform x1grid 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 largestdensity region.
Appendix C: Thermal relaxation timescale
C.1. Dustgas 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
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
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
where and c_{d} are the gas and dust specific heat capacities at constant volume and constant pressure, respectively. The radiative cooling term G^{0} in this expression is computed by replacing the term a_{R}T^{4} (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:
where , , and , with . In these equations we have dropped the δE_{r} term arising from the cG^{0} term, since the present analysis is only required in optically thin regions (more generally, for optically thin perturbations) where δE_{r} ≪ δ(a_{R}T^{4}) (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}, 0). In that case, the solution of Equation (C.4) is well approximated by δT_{g} = δT_{0} e^{−t/trel}, with
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
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^{−3} and 10^{−4}, respectively. This is a more reasonable assumption in the regions where collisional times become relevant, where we can expect centimetersized 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^{−3} and f_{dg} = 10^{−4}, 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:
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 dustdepleted 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 Planckaveraged gas absorption opacity (Kirchhoff’s law), and so we include a gasradiation interaction terms of the same form of G^{0}, specifically . This leads to the following system of equations:
which for small perturbations with respect to equilibrium becomes
where , with . In deriving Equation (C.9) we have again assumed δE_{r} ≪ δ(a_{R}T^{4}), 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,
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_{g}t_{r}/t_{r, g} ≪ 1 for our disk parameters, which leads to
The lowest accuracy of this approximation in our disk models, achieved for our highest considered value of cm^{2} g^{−1}, is 0.1% and 0.01% for f_{dg} = 10^{−3} and f_{dg} = 10^{−4}, respectively. For largely different t_{g}, t_{r, g}, and t_{cool, thin}, Equation (C.11) is equivalent to
which is the limit of the relaxation timescale considered in Malygin et al. (2017) for vanishing gasgas 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,
where and are respectively the Planck and Rosseland means of , while and are the frequencydependent 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
where
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
where t_{rel, thin} is given by Equation (C.11). Equation (C.16) also comprises the cases in which dustgas 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. Sizedependent 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
where ⟨T_{d}⟩=I[T_{d}(a)]/I[1] is the average dust temperature and
Taking , this leads to
where for b < 1/2. The radiative cooling term cG^{0} is now computed using the dust average temperature replacing with , where in the last term we dropped 𝒪(b^{2}) terms. Thus, linearizing these equations for small temperature perturbations, we obtain the system
where now , and . The resulting optically thin thermal relaxation time is
which is well approximated by
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 righthand side of Equation (C.22).
All Tables
All Figures
Fig. 1. Vertical velocity normalized by the local sound speed (top) and specific angular momentum j_{z} = R^{2}Ω (bottom) in the upper region of our highestresolution simulations at t = 234 T_{5.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 KHunstable zones in a VSIunstable axisymmetric disk. Gas rotates in the meridional plane inside of approximately constantj_{z} zones in the sense of rotation favored by the baroclinic torque, producing highshear regions that become KHunstable. 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 counterclockwiserotating 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 (∂_{R}v_{z})^{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 highestresolution simulations at t = 234 T_{5.5}. The log(p/p_{0}) values are multiplied by 5 for f_{dg} = 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 t_{crit} (top) and the orbital frequency (bottom) at r = 5.5 au for f_{dg} = 10^{−3}, assuming only radiative cooling (t_{rel} = t_{cool}), including the effect of dust collisions (), and including gas emission with defined in Eq. (13). Also shown are the corresponding values for t_{rel} = t_{crit}. 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 f_{dg} = 10^{−4}. 

In the text 
Fig. 15. VSIstability 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, t_{rel}/t_{crit} values are computed assuming only radiative cooling (t_{rel} = t_{cool}), including dust collisions (), and including as well molecular emission with (Eq. (13)). The top and bottom row correspond to the nominal and dustdepleted cases, respectively. Also shown in the leftmost panels is the location of the domain used in our RadHD 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 t_{crit} as in the righthand side of the global stability criterion (Eq. (8)). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.