Vertical shear instability in two-moment radiation-hydrodynamical simulations of irradiated protoplanetary disks

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 dust-gas collisions and molecular line emission. Methods. We modeled the gas-dust mixture in a circumstellar disk around a T Tauri star by means of high-resolution axisymmetric radiation-hydrodynamical simulations including stellar irradiation with frequency-dependent 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 high-shear 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 VSI-unstable in surface layers up to tens of au for reasonably high gas emissivities. Conclusions. The significant transfer of kinetic energy to small-scale 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 millimeter-sized dust layers while appearing vertically extended in optical and near-infrared wavelengths.


Introduction
The vertical shear instability (VSI, Goldreich & Schubert 1967;Urpin & Brandenburg 1998) is one of several proposed hydrodynamical (HD) instabilities able to produce turbulence in weakly ionized zones of protoplanetary disks where magnetic instabilities are expected to be suppressed (see, e.g., Lesur et al. 2023).Even though numerical models predict VSI-induced turbulence to lead to negligible stellar accretion compared to current estimated rates (Hartmann et al. 1998;Manara et al. 2016), this phenomenon should still produce significant vertical stirring of dust grains (Stoll & Kley 2016;Flock et al. 2017), thereby increasing their collisional velocities and diffusivities and directly impacting their ability to coagulate and form planetesimals (Ormel & Cuzzi 2007;Johansen et al. 2014;Klahr & Schreiber 2020).Furthermore, the VSI has been shown in 3D simulations to lead to the formation of small-and large-scale anticyclonic vortices (Richard et al. 2016;Manger & Klahr 2018;Pfeil & Klahr 2021), which should locally concentrate dust particles, likely accelerating the formation of planetesimals (Johansen et al. 2007;Gerbig et al. 2020) and potentially explaining structures observed in disks (e.g., van der Marel et al. 2016;de Boer et al. 2021;Marr & Dong 2022).Determining where in a disk this instability can operate is thus a necessary step toward explaining current and future observations of protoplanetary disks, as well as improving current models of planet formation.So far, no direct observations have been made that confirm the occurrence of VSI in protoplanetary disks.As proposed by Barraza-Alfaro et al. (2021), large-scale velocity perturbations produced by the VSI may be detectable in CO kinematics observations with the current capabilities of the Atacama Large Millimeter/submillimeter Array (ALMA).On the other hand, current ALMA observations of razor-thin millimeter-dust distributions in disks up to ∼ 100 au (Pinte et al. 2016;Doi & Kataoka 2021;Villenave et al. 2022) seem to contradict a possible VSI activity close to the disk midplane at such distances, as the vertical stir-ring of grains would lead to thicker dust distributions than observed (Dullemond et al. 2022).On the theoretical side, predictions of VSI-stable and unstable regions have been made (Malygin et al. 2017;Pfeil & Klahr 2019;Fukuhara et al. 2021Fukuhara et al. , 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. 2023, Paper I henceforth).
Besides the question of the disk stability, little is known about the mechanisms that lead to the saturation of the VSI.In Latter & Papaloizou (2018), it was shown via a local Boussinesq analysis of nonlinear perturbations that, once the VSI modes reach high enough velocities, parasitic Kelvin-Helmholtz (KH) modes should be produced, limiting the maximum velocities reached by the VSI modes and eventually resulting in their saturation.Even though KH eddies are, in general, not resolved in hydrodynamical simulations, perturbations in between adjacent VSI-induced flows that could be caused by the Kelvin-Helmholtz instability (KHI) can typically be seen in both 2D and 3D global simulations (e.g., Nelson et al. 2013;Manger & Klahr 2018;Pfeil & Klahr 2021).On the other hand, an alternative saturation mechanism was recently introduced in Cui & Latter (2022), where it was proposed that VSI body modes can saturate by transferring energy to small scales via resonance with pairs of inertial waves produced by nonlinear interaction (parametric instability).This instability has not yet been seen in global VSI simulations, possibly due to the expected ≳ 30 grid points per wavelength required to resolve it, estimated in that work to correspond to ∼ 300 cells per scale height (these estimates, however, depend on the wavelength of the VSI modes).It is then clear that very high resolutions are required to get a better understanding of the nonlinear behavior and saturation mechanism of the VSI.Such resolution requirements can be excessive in 3D hydrodynamical simulations, but they can be achieved in axisymmetric simulations at the cost of neglecting 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;Barraza-Alfaro et al. 2021), and thus useful information on the nonlinear evolution of the instability can be extracted from high-resolution 2D simulations.However, saturation in real disks might also be linked to the formation of nonaxisymmetric structures, such as small-and large-scale anticyclonic vortices (Richard et al. 2016;Manger & Klahr 2018;Pfeil & Klahr 2021).
In this work we investigated the growth and evolution of secondary instabilities triggered by the VSI in the axisymmetric radiation-hydrodynamical (Rad-HD) simulations of protoplanetary disks in Paper I and analyzed their connection to the saturation of the VSI.We also extended the local stability analysis considered in that work by including variations of the local thermal relaxation timescale in regions where either the dust-gas collisional timescale or the radiative emission by gas molecules become important.We then applied this analysis to make predictions of different stability regions in our disk models considering varying degrees of dust depletion, and we linked our resulting stability maps to current and future observations of protoplanetary disks.
This article is organized as follows.In Section 2 we summarize the disk models and numerical methods employed in this work.In Section 3 we characterize the secondary instabilities triggered by the VSI in our simulations.In Section 4 we extend the stability analysis in Paper I and produce stability maps for our disk models.In Section 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 Section 6 we summarize our conclusions.Additional calculations are included in the appendices.

Disk models
In this work we consider the same two disk models as in Paper I and employ the same numerical scheme therein.We solve the Rad-HD equations by means of the Rad-HD module by Melon Fuksman et al. (2021) implemented in the open-source PLUTO code (version 4.4, Mignone et al. 2007)).We use spherical coordinates (r, θ) with axial symmetry, solving as well for the azimuthal ϕ component of all vector fields.Some quantities in this work are instead expressed in cylindrical coordinates (R, z).We model the gravitational potential and heating produced by a T Tauri star of mass M s = 0.5M ⊙ , effective temperature T s = 4000 K, and radius R s = 2.5 R ⊙ .The heating rate of the disk gas-dust mixture due to stellar irradiation is computed via ray-tracing using frequency-dependent absorption opacities taken from (Krieger & Wolf 2020, 2022), assuming small dust grains dynamically well coupled to the gas with radii between 5 and 250 nm.The frequency-averaged values of these opacities are used for the computation of absorption and emission of reprocessed radiation.
Initial hydrostatic conditions are obtained by iterating the computation of temperature and density distributions via the radiative transfer method in the code and a numerical solution of hydrostatic equations, as detailed in Melon Fuksman & Klahr (2022).The gas column density is in every case kept as Σ(r) = 600 g cm −2 (r/au) −1 .Once initial conditions are obtained in the domain (r, θ) ∈ [0.4,100] au× [π/2 − 0.5, π/2 + 0.5], the Rad-HD equations are solved in the smaller domain (r, θ) All presented simulations are labeled in Table 1 together with their corresponding parameters.For the computation of absorption opacities for both stellar irradiation and radiative transfer of reprocessed infrared photons, we assume a constant dust-to-gas mass ratio f dg of small (< 0.25 µm) grains.We consider two f dg values, representing a nominal case with f dg = 10 −3 and a dust-depleted case with f dg = 10 −4 .Assuming a total dust-to-gas mass ratio including all grain sizes of ρ tot d /ρ = 10 −2 , these values correspond respectively to 10% and 1% of the total dust mass being in the small-grain population, which for the grain size distribution dn ∼ a −3.5 da assumed for the opacity computation correspond to maximum grain sizes of 19 µm and 1.8 mm, respectively.The chosen disk region is either vertically thick or thin and has different stability regions depending on the dust content: while in the nominal case the entire domain is VSI-unstable, in the dust-depleted case the disk is stable at the middle layer located below the irradiation Notes.Dust-to-gas mass ratio of small grains (f dg ), resolution (Nr × N θ ), approximate number of cells per scale height (H/∆r), and total simulation time (t f ) in units of the Keplerian orbital time at 5.5 au (T5.5).Also shown for reference are the normalized radiative cooling time Ωt cool and the ratio t cool t crit between the cooling time and the local critical value for instability in Equation (7) (Section 4.1), shown at the maximum height in the domain (z/r ≈ 0.3) and at the height z = ΓH 2 determining the vertically global stability criterion in Equation ( 8).
surfaces, defined as the regions of unity optical depth for photons emitted at the star.As detailed in Paper I, this results from the different space-dependent cooling timescales in the domain for varying dust content (see also Section 4).
Radiative transfer of reprocessed photons is modeled by evolving the frequency-integrated radiation energy density and flux employing the M1 closure by Levermore (1984).To increase the largest possible time step, the disparity between the radiation and HD timescales is reduced by employing an artificially small value ĉ < c of the speed of light to advect the radiation fields.In this work we take ĉ = 10 −4 c, which, as verified in Paper I via analytical and numerical testing, does not introduce unphysical phenomena in our simulations.Together with the employment of a reduced radial domain compared to the disk size, this allows us to achieve the high resolutions needed to investigate secondary instabilities.We employ a maximum resolution of N r × N θ = 1920 × 2048, which corresponds to approximately 200 cells per pressure scale height H (see Table 1) with an aspect ratio of 1 : 1.We limit numerical diffusion by employing third-order WENO spatial reconstruction (Yamaleev & Carpenter 2009;Mignone 2014) and HLLC-type Riemann solvers for both radiation (Melon Fuksman & Mignone 2019) and matter fields (Toro 2009).Dirichlet boundary conditions are applied by fixing all fields to the hydrostatic initial conditions at the ghost zones.Unless otherwise stated, all results presented in this work correspond to our highest-resolution runs.

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 uniform-j z bands, in which the initial vertical shear is significantly reduced.This can be seen in Figs. 1 and 2, showing velocity and j z distributions after VSI saturation.
As shown in these figures and schematized in Fig. 3, the gas inside the constant-j z bands rotates clockwise above the midplane (in a reference frame in which the orbital velocity enters the meridional plane) and counterclockwise below the midplane, with minimum speed ad the center of the bands and maximum speed at the edges.Due to the abrupt reversal of the vertical velocity at the interfaces between adjacent bands, significant meridional shear is created at such regions, eventually triggering the KHI.While this instability typically involves a balance of the destabilizing effect of shear and the stabilizing effect of buoyancy (e.g., in atmospheric physics), the latter is suppressed in our considered system by the fast cooling required for the VSI to operate.Instead, stabilization results from the radial angular momentum stratification, and a KH-stability criterion can be written in terms of a Richardson number in which the Brunt-Vaisälä frequency is replaced by the epicyclic frequency (Section 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 (Section 3.3).
Various aspects of both the large-scale structure of the vertical flows in the nonlinear stage of the VSI and the amplification of small-scale eddies in our simulations can be explained in terms of the generation of vorticity in rotating baroclinic flows.Locally, the time evolution of the azimuthal vorticity ω ϕ = (∇ × v) ϕ depends on a combination of three terms, denoted τ a , τ b , and τ c (Section 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 τ c ∝ ∂ z j 2 z as constant-j z regions start forming, contributing to the clockwise rotation of the large-scale flows at z > 0 and counterclockwise rotation at z < 0. Although the term τ a is generally dominant over the other two, the acceleration or deceleration of vortices advected with the flow only depends on the balance of τ b and τ c .Thus, in constant-j z regions, τ b is the dominant term contributing to vortex acceleration, which results in an acceleration of vortices with sgn(ω ϕ ) = sgn(τ b ) and, vice versa, a deceleration of vortices with sgn(ω ϕ ) = −sgn(τ b ).Even though the KHI mainly produces vortices rotating in the direction disfavored by the baroclinic torque, its nonlinear evolution also results, via τ a , in the formation of counter-rotating vortices in regions with sgn(ω ϕ ) = sgn(τ b ).These vortices can then be accelerated by τ b in the constant-j z bands.Since sgn(τ b ) = sgn(z), this leads to the formation of small-scale vortices rotating clockwise at z > 0 and counterclockwise at z < 0, which are accelerated up to significant fractions of the local sound speed (Sections 3.4 and 3.5).Further acceleration of these vortices is likely produced via advection from regions adjacent to the shear layers with sgn(ω ϕ ) = sgn(τ b ).The size of these vortices is limited by the width of the constant-j z regions.In axisymmetric simulations, in which such regions can grow in size unimpeded by non-axisymmetric instabilities, the baroclinically amplified vortices can eventually grow up to large sizes (∼ 1 au) dominating the gas dynamics, as detailed in Section 3.6.

Kelvin-Helmholtz instability
We begin by analyzing the triggering of the KHI in between uniform-j z regions.As shown in Figs. 1, 3, and 4, the rotation pattern inside of the approximately uniform-j z bands causes downward-moving gas parcels at z > 0 to experience significant shear on their larger-r edges, while below the midplane shear is maximum on the smaller-r edges.For f dg = 10 −3 , in which case the VSI flows cross the midplane, this produces an intercalated pattern of KH-unstable regions, each of them occupying half of the vertical domain (see also Fig. 5).The same occurs for f dg = 10 −4 in the upper VSI-unstable layers.The resulting eddies are clearly highlighted by ω ϕ , which also shows their sense of rotation.We show this quantity in Fig. 4 next to the normalized vertical velocity in a region close to the midplane for run dg3c4_2048.
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 z 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 2023).While in the adiabatic limit , where N R is the radial Brunt-Vaisälä frequency, ω 2 tends to κ 2 R for fast thermal relaxation compared to epicyclic motion (t cool ≪ Ω −1 , as shown in Section 4.1), leading to Eq. (1).In the context of rotating stars, an equivalent expression is obtained in Maeder et al. (2013) in the limits of instant cooling and uniform molecular weight.
A necessary condition for KH-stability is Ri < 1/4 (Chandrasekhar 1961;Maeder et al. 2013).As soon as VSI modes are formed, this condition is only violated in the high-shear regions in between vertical flows (shown in Fig. 5 in the nonlinear stage).At those locations, KH eddies are visible within a few orbits after this condition is met, as shown in Fig. 6.An alternative instability criterion based on a linear analysis of perturbations to the VSI modes was derived in Latter & Papaloizou (2018), requiring that kV /Ω ≳ 1.38 to trigger the KHI, where k is the wavenumber of the VSI modes and 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.
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.4. Velocity and vorticity distributions in a region close to the midplane in run dg3c4_2048 after 300 orbits.Red and blue regions in the vorticity map correspond to clockwise and counterclockwise rotation, respectively.KH eddies on the lower hemisphere rotate clockwise, while their nonlinear evolution produces counterclockwise-rotating eddies that are amplified by the baroclinic torque, visible in that region as dark blue spots.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), andFlock et al. (2020) (see Paper I), we can expect the KHI to produce azimuthal eddies in the same regions.Since j z increases with radius, such eddies should have prograde rotation with respect to the disk.Therefore, long-lived anticyclonic vortices such as those in Manger & Klahr (2018) and Pfeil & Klahr (2021), if formed via KHI, can only result from secondary eddies produced by the nonlinear evolution of this instability in adjacent regions with a negative vorticity perturbation, similarly to the formation of counterrotating eddies in the meridional plane (Section 3.4).

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 Fourier-transformed meridional momentum and velocity distributions in an upper disk region exhibiting VSI growth for both f dg values.Energy spectra between 6 and 37 orbital times are shown in Fig. 7 as a function of kH, where k = ||(k 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.
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 constant-j z bands.
The beginning of the saturation phase, occurring at ∼ 15 − 25 orbits depending on f dg , coincides with a deviation from the initial spectrum's shape.At that time, a 'hump' is formed at kH ∼ 50 − 200 due to the formation of KH eddies in that size range.The spectra converge to a new shape at saturation, with a short plateau following the main VSI peak up to kH ∼ 50.This plateau is followed by a shallower decrease with k than in the linear growth stage, with E(k) ∼ k −1.7±0.2 at kH ∼ 50 − 100 for f dg = 10 −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 high-k end of the VSI-induced structures and possibly a direct enstrophy cascade (a discussion on the energy transport over scales is given in Section 5.3).The precise limits of this range typically depend on the fluid's velocity, the density stratification, and the Richardson number, but for this analysis it is enough to consider that its upper limit is some number on the order of 1/d, where d is the typical width of the shear layers (Chandrasekhar 1961).In the region considered for the computation of the energy spectra, we typically have d ∼ 0.02−0.03H,corresponding to wavenumbers kH ∼ 200 − 300, coinciding with the region where we see the slope transition.In this regard, we must mention that these layers are typically resolved by ∼ 10 cells, and thus it is possible that numerical diffusion is enlarging their typical width with respect to its zero-diffusion value, artificially shifting the slope transition to lower wavenumbers.
On the other hand, numerical diffusion can also contribute to the steepening of the energy spectrum at high k.This occurs in the 2D homogeneous turbulence simulation made by Seligman & Laughlin (2017) with PLUTO, also using third-order Runge-Kutta scheme and WENO3 reconstruction.In that test, numerical diffusion leads to a steeper spectrum than the predicted E(k) ∼ k −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 6 orbits 9 orbits 11 orbits 16 orbits 18 orbits 23 orbits 25 orbits 30 orbits 32 orbits 37 orbits at kH ∼ 600 − 800 (about 3 cells per wavelength), which could indicate that the dissipation range is instead closer to this region.Yet, the relevance of numerical diffusion on the high-end of our VSI spectra can only be assessed via a proper comparison with 2D forced turbulence tests injecting energy at k ≲ 2π/d, as we plan to do in a future study.

Vorticity generation
In our highest-resolution simulations (N θ = 1024 and 2048), we observe a formation of small meridional vortices in the VSI-active regions that survive up to 1 − 10 orbital timescales and get accelerated up to hundreds of m s −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 VSI-active regions where the accelerated vortices form (Section 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).
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, v p = v R R + v z ẑ (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 κ 2 R quantifying the vertical squared specific angular momentum gradient.The first and second terms in the right-hand side of Equations ( 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 κ 2 z 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 (5) 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  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(ω ϕ ) = −sgn(z) 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(ω ϕ ) = sgn(z).This process can be seen in Fig. 6, where it is shown that the shear layer (with sgn(ω ϕ ) = −sgn(z)) is adjacent to a region of opposite ω ϕ -sign (see also Fig. 9).Once the KHI is triggered1 , such regions form vortices with sgn(ω ϕ ) = sgn(z).Eddies created in this way can then be accelerated in the constant-j z zones by τ b .
In the vorticity distribution for f dg = 10 −4 (Fig. 8), it can be seen that the amplified eddies are always adjacent to the zones with sgn(ω ϕ ) = sgn(z) located at the smaller-r 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 Fig. 9. Advection, baroclinic, and centrifugal terms (τa, τ b , and τc, respectively) determining the vorticity evolution (Equations ( 3) and ( 5)), computed at the same snapshot as Fig. 8.
In summary, the observed vortices with sgn(ω ϕ ) = sgn(τ b ) accelerated in the uniform-j z bands are seeded at the band interfaces via τ a by the nonlinear evolution of the KHI2 .The amplified vortices are rapidly destroyed when advected into the shear layers (possibly with some contribution of τ c ), and thus their size is limited by the width of the uniform-j z regions.
The prevailing baroclinic torque at the constant-j z bands may also explain why the VSI-induced large-scale flows also respect a rotation pattern with sgn(ω ϕ ) = sgn(τ b ).In its linear stage, the VSI transports gas across constant-j z surfaces, as determined by the unstable directions of motion predicted by linear theory (James & Kahn 1970;Knobloch & Spruit 1982;Klahr 2023).Since ∇j z decreases with height and increases with radius, this transport reduces the j z gradient between gas parcels moving away from the midplane and their adjacent outer parcels moving toward it, which starts forming constant-j z bands.Thus, the baroclinic acceleration may explain why the sat-urated 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.

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 (Equation 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 (Section 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.
An explanation for the fact that Ω v = O(Ω) 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 O(z/R) 2 terms, the sum of these accelerations is j 2 z /R 3 − GM s /R 2 , 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 constant-j z bands in Fig. 8. Assuming that relative temperature variations are much smaller than relative pressure variations, which can be done as long as the cooling time is much smaller than the dynamical time (t cool ≪ Ω −1 ), we can solve the resulting equation for p by writing ρ = p/c 2 s with constant c 2 s , which for small radial departures , 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 half-width of a VSI mode, this balance results in a relative pressure decrease ∆p/p ∼ . Thus, vortices can only be accel- erated as long as they can produce deeper pressure minima, but this process is limited by the pressure increase caused by the continuously enforced radial contraction of the constant-j z bands.We can then expect the saturation of the vortex amplification to be reached when the pressure perturbation caused by both processes is of a similar order of magnitude, which leads to ∆p/p ∼ , 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 (Equation (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.

Long-time evolution
In axisymmetric VSI simulations, bands of approximately constant j z can grow up to large widths (∼ 1 au in our simulations) for long enough times.This occurs due to the mixing of angular momentum enforced to occur in the meridional plane, given that these structures are unaffected by non-axisymmetric modes that should prevent them from growing to such extents in 3D.Since the size of the baroclinically amplified eddies is only limited by the width of the bands, these eddies are allowed to occupy larger regions as the bands grow.For f dg = 10 −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 Section 5.3).
Even if the resolution is not large enough to initially resolve the amplified eddies, as it occurs for N θ = 512, these are eventually resolved for sufficiently long times due to the increase of their typical size resulting from the growth of the bands.The transition into this regime occurs at earlier times for increasing resolution, as evidenced by the time evolution of the kinetic energy (Fig. 12).This shows that the mixing of angular momentum in between constant-j z regions is more efficient when smaller scales are resolved, and in particular when the shear in between these regions is larger, which may result from the angular momentum mixing between adjacent bands produced by the KH eddies, evidenced, for instance, in Fig. 1.Moreover, this also shows that the observed transition is not caused by numerical diffusion of angular momentum, as that would lead to the opposite trend.
A transition into an eddy-dominated regime is also seen in isothermal versions of our setup (not shown here) and in the isothermal runs by Klahr et al. (2023).This shows that, even though the increased optical depth of the large eddies with respect to the VSI modes and the subsequent increase of the local cooling timescales in the Rad-HD simulations could potentially trigger the COS, that instability is not what produces the transition into that regime.The required high resolutions and long times to obtain this transition (in our case ∼ 3000 and 800 orbits for 50 and 100 cells per scale height, respectively, with our employed numerical scheme) explain why this phenomenon has not been seen in other works.
The mentioned expected differences with nonaxisymmetric simulations suggest that this transition is merely a consequence of the enforced axisymmetry.Thus, the long-term nonlinear evolution of the VSI can only be studied in high-resolution 3D simulations (see Sections 5.3 and 5.4).

Influence of collisional timescales and gas emissivity
In Paper I, we analyzed the stability of our Rad-HD simulations in terms of the local criterion by Urpin (2003), which predicts instability as long as where t rel is the local thermal relaxation time, which in the Rad-HD simulations coincides with the radiative cooling time t cool , while N z is the vertical Brunt-Väisälä frequency (e.g., Rüdiger et al. 2002).The localized suppression of the VSI in the middle layer of our dust-depleted disk is then explained by the fact that this criterion is not met in that entire region.On the other hand, the vertically global stability criterion by Lin & Youdin (2015), predicting instability for where Γ is the heat capacity ratio, q = ∂ log T ∂ log R , and Ω K is the Keplerian rotation frequency, correctly predicts the stability properties of the middle layer but does not contemplate cases in which the disk surface layers, where the assumptions of vertically uniform temperature and cooling timescale break down, have different stability conditions than the midplane.
The stability analysis in Paper I relies on the assumption that gas and dust particles are in local thermal equilibrium.However, this is not the case at the surface layers, whose low densities make collisions between dust and gas particles infrequent enough that these species may have different temperatures.Considering this, we can still estimate t rel in those regions by means of a linear analysis of small perturbations of separate dust and gas temperatures.We follow this approach in Appendix C in the same way as Barranco et al. (2018), with the difference that our largest considered dust grains are micron-sized, whereas in that work cm-to meter-sized grains were assumed.Despite this difference, the functional form of the obtained timescale is identical as in that work3 : 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, ρ tot d is the total local dust density, a max and a min are respectively the maximum and minimum grain radii, and v g = 8k B T πµu is the mean thermal speed of gas particles.For an assumed total dust-to-gas ratio ρ tot d /ρ = 10 −2 , our assumed dustto-gas 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 Equation ( 9), are shown in Figures 13 and 14 (labeled as κ g P = 0), 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 .
This situation changes if cooling by molecular emission is not negligible.This effect should be most important at the surface layers, where the typical temperatures of up to ∼ 200 − 250 K in our domain may increase the amount of possible free-free and bound-free emission (and absorption) mechanisms (e.g., Freedman et al. 2008;Malygin et al. 2014) with respect to the cooler midplane.If optically thin radiative emission by gas molecules is considered, it can be shown (Appendix C) that the linear relaxation timescale is modified as where t rel,thin is the optically thin limit of this expression, given by where λ dg P and λ dg R 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 A&A proofs: manuscript no.ms , assuming only radiative cooling (t rel = t cool ), including the effect of dust collisions (κ g P 0 = 0), and including gas emission with κ g P 0 defined in Equation ( 13).Also shown are the corresponding values for t rel = tcrit.The gas opacities reported in Freedman et al. (2008), Malygin et al. (2014), and Malygin et al. (2017) for solar abundances correspond to κ g P 0 ∼ 10 −1 − 10 0 .
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 t r,g = cgT cκ g P [4a R T 4 +b g P (a R T 4 −Er)] , where c g = k B (Γ−1)µu is the gas specific heat capacity at constant volume, κ g P (T ) is the Planck-averaged gas absorption opacity, and b g P = d log κ g P d log T .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 (≲ 75K 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 κ g P 0 ∈ [10 −4 , 1] parameterizes the uncertainty in the gas composition.For the computation of the combined dust-gas Rosseland opacity, we only consider the contribution of dust absorption, assuming that our densities are low enough that absorption lines can be disregarded in this calculation (see, e.g., Freedman et al. 2008).
The obtained relaxation times at 5.5 au for different κ g P 0 values are shown in Figures 13 and 14.These figures show that the VSI can still operate at the surface layers provided gas opacities are large enough (κ g P 0 ≳ 10 −2 in our models).For κ g P 0 ≳ 10 −1 , the gas opacities are large enough that even the middle layer is VSI-unstable in our dust-depleted case.However, we can only expect the temperature and density distributions to be roughly unchanged for varying opacity if κ g P is smaller than both f dg κ d P and f dg κ d R , where κ d P and κ d R 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 κ g P 0 ≲ 10 −1 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 κ g P 0 ≲ 10 −1 .Thus, our conclusion regarding the possible instability of the upper layers remains the same.

Global disk stability
In this section we extend our stability analysis to the full domain of the hydrostatic disk models used as initial conditions (see Paper I), for which an estimation of the dominant VSI wavenumbers is required.To compute t cool , we assume the maximum dominant wavenumbers in our simulations to be approximately k ∼ 70/H (Fig. C.3 in Paper I) in all disk regions.For larger radii than in our simulation domain, this does not introduce any error, given that our obtained VSI modes are radially optically thin (Paper I) and therefore cool at the scale-independent optically thin rate  7)).From left to right, t rel /tcrit values are computed assuming only radiative cooling (t rel = t cool ), including dust collisions (κ g P 0 = 0), and including as well molecular emission with κ g P 0 ̸ = 0 (Equation ( 13)).The top and bottom row correspond to the nominal and dust-depleted cases, respectively.Also shown in the leftmost panels is the location of the domain used in our Rad-HD simulations (white dashed box).The gas opacities reported in Freedman et al. (2008), Malygin et al. (2014), andMalygin et al. (2017) for solar abundances correspond to κ g P 0 ∼ 10 −1 − 10 0 .given by Equation 16of that work, while for larger radii the optical depth along any wavelength only decreases with radius.For smaller radii, on the other hand, we can instead expect the dominant wavelengths in our simulations to become optically thick close to the midplane at some minimum radius, in which case VSI modes of such wavelengths cool via diffusion.If t rel > t crit at some location for our considered k value, that just means that only modes with larger wavenumbers can be unstable at that point.Considering that the strength of the instability decreases in our simulations for increasing average wavenumber (Paper I), it is also likely that this is also translated into either a partial or total suppression of the VSI.The resulting t rel /t crit values computed in this way are shown in Fig. 15 between 0.6 and 90 au (this is far enough from the radial boundaries to avoid boundary artifacts, as we used zero-gradient conditions at those boundaries to construct the hydrostatic models).Interpreting these distributions as approximate VSI stability maps, we obtain that the unstable surface regions extend up to a few tens of au for high enough gas opacity (κ g P 0 ≳ 10 −2 ).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 κ g P 0 ≳ 10 −1 , in which case an instability region is formed between 1 and 5 au due to the decrease of the optically thin cooling time with κ g P 0 .
For comparison with Malygin et al. (2017), we also computed stability maps based on a local evaluation of Equation 8 as done in that work, which we show in Figure 16.Both criteria coincide in that region for the following reasons: (I) the assumptions of the global criterion are met below the irradiation surface, and (II) the vertically global critical cooling time (right-hand side of Equation ( 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 κ g P 0 .
We defer a thorough comparison of local stability criteria with hydrodynamical simulations to future works.Differences with the stability maps at the midplane in Flock et al. (2017Flock et al. ( , 2020) ) and Pfeil & Klahr (2019, 2021), in particular the fact that an inner stable region and an outer stable region are obtained in the latter, are explained by the different assumed dust masses, opacity laws, and dominant wavelengths in those works and this one.The radiative cooling time is in general diffusion-dominated in inner regions and thus it increases with increasing opacity (e.g., Lin & Youdin 2015), which explains the inner stable region in Pfeil & Klahr (2021) and this work.Conversely, far enough from the star that perturbations are optically thin, t cool grows for increasing radius as the photon mean free path increases (see Paper I), eventually reaching the condition t cool > t crit , as shown in Fig. 15 (see also Malygin et al. (2017) and Dullemond et al. (2022)).The fact that the cooling times in Flock et al. (2017Flock et al. ( , 2020) ) are estimated to be diffusion-dominated 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 dust-depleted simulations.These regions appear as a consequence of the reduction of the thermal relaxation times above the irradiation surface induced by our temperature and density stratification, and thus they do not appear in stability maps of vertically isothermal disks (Malygin et al. 2017) nor in the models by (Pfeil & Klahr 2019), in which the temperature is computed via an effective 1D treatment of vertical radiative diffusion.Since gas radiative emission can only reduce the relaxation timescale for large enough temperatures (≳ 100 K is a rough estimate based on Malygin et al. (2017)), the radial extent of these unstable regions depends not only on the chemical composition of the gas but also on the dust opacity at the disk upper layers and the temperature of the central star, occupying larger radii for higher temperatures as evidenced by Equation 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).

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 Flores-Rivera et al. (2020).

Saturation mechanism
The reduction of the energy proportion in large scales once the KHI is triggered (Section 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 ( The distribution of the gas velocity in relation to the location of the KH-unstable zones provides a rough picture of this saturation mechanism: once the VSI modes reach sufficiently high velocities to trigger the KHI, the loss of energy to KH eddies in the unstable shear layers acts as an effective friction between the vertically moving flows (see Fig. 3), limiting the maximum kinetic energy that these can gain.Since in our simulations only a very small fraction of the energy is in the highest wavenumbers (Fig. 7), we expect numerical dissipation to play a minor role in the total energy dissipation as long as a direct energy cascade does not occur (see Section 5.3).In real 3D disks, viscous dissipation could become important if the KHI triggers a turbulent cascade leading to dissipation down to molecular viscous lengthscales.Instead, in our simulations, it is likely that energy is radiatively dissipated at the medium to large scales of the VSI flows and KH eddies.In particular, some dissipation may be produced by the baroclinic deceleration of the primary KH eddies, which grow with sgn(ω ϕ ) = −τ b .
Another possible saturation mechanism is the parametric instability considered in Cui & Latter (2022), whose longest wavelengths are expected to require ∼ 30 grid points per VSI wavelength to be resolved.In our simulations, the wavelength of the VSI modes is well resolved with λ/∆r ∼ 15 − 75, and in particular ∼ 75 in Fig. 4. In that figure we do indeed see small-scale perturbations inside of the VSI modes, but we do not see significant amounts of kinetic energy being transferred to them (see also Fig. 8).Rather, the wavelength of the VSI modes grows over time due to the mixing of angular momentum between approximately uniform-j z bands (Paper I and Section 3.6), and so even if the ripples inside of the VSI modes are produced via parametric instability, that process is not able to destroy such modes.On the other hand, it is unclear whether the observed ripples are instead sound waves emitted at the KH-unstable regions.Thus, we cannot conclude that a parametric instability occurs in our simulations due to the current lack of unequivocal signs of its occurrence in its nonlinear stage, but even if it does, its dynamical effect in our simulations appears to be marginal.

Energy transport between scales
The transfer of energy from VSI to KHI modes does not imply the occurrence of an energy flux from large to small scales as it occurs in 3D turbulent cascades.Instead, as first studied by Kraichnan (1967), 2D flows typically exhibit a split cascade transporting energy toward large scales and enstrophy toward small scales (some experimental verifications of such a cascade can be found in Bruneau & Kellay 2005;Kelley & Ouellette 2011).Thus, it is possible that a similar phenomenon occurs in our simulations.In this picture, energy is continuously injected in a broad range of wavelengths due to the geometry of the large-scale VSI flows (as we have verified) and transported toward small k starting from the largest k allowing for KHI growth (Section 3.3).However, it must be noted that our simulations are not strictly two-dimensional.In particular the mechanisms of vortex stretching (Appendix A) and strain selfamplification4 , often linked with the development of direct energy cascades in 3D (see, e.g., Carbone & Bragg 2020;Johnson 2021), are absent in 2D flows but not in axisymmetric flows.Moreover, as reviewed by Alexakis & Biferale (2018) (see also the discussion in Sengupta & Umurhan 2023), both direct and inverse energy cascades are possible in 3D (non-axisymmetric) turbulent flows when the effects of rotation, stratification, and even compressibility become important.
In our simulations, a hint possibly indicating transport of energy toward lower wavenumbers is given by the observation that vortices merge and grow in size, likely contributing to the formation of large-scale meridional vortices at long times (Section 3.6).Also resulting in the accumulation of energy in increasingly large scales is the growth over time of the uniform-j z bands, likely via in-plane angular momentum mixing between adjacent bands.If an inverse energy cascade occurs, then energy should be radiatively dissipated at the medium and large scales.In this scenario, it is possible that the slopes measured for 5 ≲ kH < 100 − 200 result from a combination of an inverse cascade and the geometry imposed by the VSI modes.In this regard, we must clarify that, despite the spectra decay close to k −5/3 for f dg = 10 −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 k-range dominated by KHI modes.At such scales, anisotropy can still be observed even for the velocity components perpendicular to the VSI flows, which track mostly KHI modes.Moreover, the spectra deviate from that functional form after several tens of orbits, as the constant-j z regions grow.This leads to a steepening over time, with E(k) eventually oscillating in that wavenumber range around ∼ k −2.4±0.2 for f dg = 10 −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 uniform-j z regions, nor the evaluation of the slopes of the energy spectra, are sufficient to estimate the influence of rotation, gravity, weak compressibility, and baroclinicity on the transport of energy and enstrophy between scales.Some information on this can likely be obtained by computing the energy and enstrophy fluxes resulting from each of these phenomena, as done in the multiple analyses in Alexakis & Biferale (2018) for incompressible flows.We defer such an analysis to a future work on the topic.

Caveats of axisymmetric models
Once the KHI is triggered in our simulations, amplified eddies never cease to form and dissipate, and thus they are able to continuously limit the maximum kinetic energy of the vertical flows even hundreds of orbits after saturation.Unfortunately, our simulations are insufficient to make predictions on the disk's state after thousands of orbits, given that the artificially enhanced growth of uniform-j z bands due to vertical angular momentum mixing in 2D eventually leads in some cases to a complete breakdown of the VSI flows and domination of large-scale azimuthal vortices, as shown in Section 3.6 (see also Klahr et al. 2023).If the growth of these vortices is enhanced via an inverse energy cascade, this process should likely be inhibited in 3D simulations, provided that a direct energy cascade occurs in that case.With enough azimuthal resolution, we expect the formation of wide uniform-j z bands (Paper I and Papaloizou & Pringle 1984) to be inhibited by the growth of non-axisymmetric modes, which should disallow the growth of large-scale vortices at long times.
On top of this, it is even unclear whether small-scale eddies can be baroclinically amplified and survive for several orbits if axisymmetry is not imposed.In principle, all of the ingredients needed for the formation of meridional vortices could still be present in that case, given that the VSI driving mechanism is axisymmetric.As a result, the large-scale VSI-induced flows maintain some degree of axisymmetry in some 3D simulations (e.g., Richard et al. 2016;Flock et al. 2020;Barraza-Alfaro et al. 2021) and are likely prone to KHI, as the velocity profiles in Manger & Klahr (2018) (Figure 10) seem to indicate.Axisymmetry, however, could be artificially enhanced by underresolving the ϕ-direction, and thus it is in general unclear to which extent non-axisymmetric modes can erase the initial axisymmetry of the VSI modes.Therefore, the question remains of whether bands of reduced ∇j z are formed in 3D, and even whether meridional vortices are destroyed by 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), andFlock et al. (2020) (see also Fig. A.1 in Paper I), but asserting the universality of this phenomenon requires further careful investigation.
Determining whether a baroclinic amplification of KH eddies can occur in non-axisymmetric disks has the additional difficulty that cell aspect ratios of ∼ 1 should be achieved to avoid purely axisymmetric effects, together with the high resolutions of at least ∼ 100 cells per scale height needed to resolve the vortices, which is computationally rather expensive.Interestingly, radial resolutions of H/∆R ∼ 100 − 150 are reached in the β-cooling 3D disk simulations by Richard et al. (2016) for an aspect ratio of H/R ∼ 0.05 similar to those in this work (see Paper I), albeit with a steeper temperature profile (T ∼ R −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, non-axisymmetric structures grow, forming azimuthal vortices in the regions of lowered ω z .This phenomenon can possibly be attributed to the Rossby Wave Instability, as suggested by the authors.The shear layers exhibit some ripples in the meridional plane (e.g., Fig. 13 in that work), which are possibly produced by the KHI.
For all parameters considered in Richard et al. (2016), the reduction of ω z in between the shear layers with respect to its Keplerian value appears lower than the value ∆ω z /Ω = −1/2 that would result from the formation of uniform-j z zones (for sufficiently small ∂ ϕ v R ).This may be a consequence of the non-axisymmetric transport of angular momentum, which should tend to counteract the formation of axisymmetric uniform-j z regions by the VSI (Paper I and Papaloizou & Pringle 1984).However, the fact that ∆ω z < 0 is still indicative that regions of lowered-∇j z are formed (before non-axisymmetric structures are formed, ω z ∝ ∂ 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 non-axisymmetric modes, as it may be the case near the midplane in Fig. 8 of Richard et al. (2016).
Despite the mentioned numerical challenges, the 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 non-axisymmetric modes.

Stability regions and observational consequences
The fact that our disk models are VSI-stable close to the midplane from a minimum radius is consistent with current dust continuum observations of protoplanetary disks according to the argument by some authors (Flock et al. 2017;Dullemond et al. 2022) that a suppression of the VSI, at least in the outer disk, is required to explain the thin distributions of millimeter-sized dust seen in ALMA dust continuum images of protoplanetary disks, for instance in HD 163296 at r ≈ 100 au (Doi & Kataoka 2021) and in Oph163131 up to ∼ 170 au (Villenave et al. 2022) (see also Pinte et al. 2016).In particular, a suppression of the VSI at large radii could explain why the dust scale height estimated in Doi & Kataoka (2021) is close to H at 68 au but much smaller at 100 au.The stable regions in the outer parts of our disk models grow in size as the density of small grains is reduced, which is why it has been proposed that these can naturally occur as a consequence of small-dust depletion due to coagulation (Fukuhara et al. 2021;Dullemond et al. 2022;Pfeil et al. 2023).To this picture, our local analysis adds the possibility that the VSI may still operate in surface layers even at distances from the star where it is suppressed at the midplane.
Considering the works by Stoll & Kley (2016) and Flock et al. (2017), we expect a VSI activity confined to regions above the irradiation surface to produce significant vertical stirring of small dust grains, and therefore to limit vertical settling in such regions.We conjecture that this phenomenon can contribute to the vertically extended appearance of edge-on disks at optical and near-infrared wavelengths (Villenave et al. 2020) even at radii where turbulence appears weak or nonexistent at the disk midplane (Villenave et al. 2022), as we expect the vertical stirring of dust grains to push the location of the τ = 1 surface for photon scattering at µm wavelengths away from the midplane.This scenario can be tested in numerical experiments by computing the distributions of different dust populations in a dust fluid approach (e.g., Krapp et al. 2022) and using these to produce synthetic dust continuum and scattered light images.It must be noted, however, that the distributions of small grains in the disk upper layers may also depend on nonideal MHD effects and magnetized winds (Riols & Lesur 2018;Booth & Clarke 2021;Hutchison & Clarke 2021).
Depending on whether the τ = 1 surfaces of 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 Barraza-Alfaro et al. (2021).In the same work, it is concluded that the nonthermal spectral broadening of molecular lines produced by the VSI should be negligible compared to current ALMA capabilities, which is consistent with current measured upper bounds (e.g., Teague et al. 2018;Flaherty et al. 2020), and which also applies to an instability localized at surface layers.

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

Conclusions
In this work we studied the growth and evolution of secondary instabilities parasitic to the VSI modes and their relation with the VSI saturation in the axisymmetric Rad-HD simulations presented in Paper I. We also extended the VSI-stability analysis of Paper I including the effects of dust-gas thermal equilibration and gas molecular emission, and applied it to construct stability maps in our disk models.Our main findings can be summarized as follows: 1.In its linear stage, the VSI produces adjacent upwardand downward-directed flows which are initially disconnected.The transport of angular momentum across constant-j z surfaces reduces the angular momentum gradient between flows moving away from the midplane and the immediate outward flows moving toward the midplane, causing the baroclinic torque (τ b ) to prevail over the centrifugal torque (τ c ) in those regions.Such adjacent flows eventually become connected in upper disk regions, either due to the angular momentum excess of outward-moving gas parcels or due to baroclinic acceleration.
2. Connected adjacent vertical flows exchange angular momentum in such a way that bands of approximately uniform j z are formed after a few orbits.The initial τ c is consequently suppressed inside of these bands, while the initial τ b is largely unaffected, leading to the production of vorticity in those regions by baroclinic torquing.This may explain why the saturated VSI-induced flows survive despite the reduced angular momentum gradient inside of the bands.
3. While the meridional shear is also reduced inside of the formed uniform-j z bands, it is maximum at the interface between them, with large enough values to overcome the stabilizing effect of the radial j z gradient.As a result, the KHI is triggered at those layers.This instability transfers significant kinetic energy from the VSI flows to KH eddies, acting as an effective friction between the vertical flows and limiting the maximum kinetic energy these can gain due to the VSI driving, possibly leading to the saturation of the VSI.
4. Large resolutions favor the KHI by increasing the shear in between uniform-j z bands and reducing the numerical diffusivity.Together with the increase of the number of KH-unstable interfaces with resolution, this may explain why both the average kinetic energy and the Reynolds stresses after saturation decrease with resolution without converging in our resolution range.We require about 100 cells per scale height to unequivocally identify the KH eddies in vorticity distributions, but the velocity perturbations in the VSI modes show that this instability is triggered in our lower resolution runs as well.Similar signs of developing KHI in other works suggest that this is a rather general phenomenon involved in the VSI saturation.
5. Due to the flow pattern produced by the VSI, the initially produced KH eddies rotate in the direction disfavored by τ b , namely, with sgn(ω ϕ ) = −sgn(τ b ).However, the nonlinear evolution of the KHI at the shear layers also generates counter-rotating eddies in adjacent zones with sgn(ω ϕ ) = sgn(τ b ).As a result, some of those vortices get accelerated by τ b in the uniform-j z zones, possibly with some contribution of vorticity advection from regions near the shear layers.Vortices accelerated in this way reach Mach numbers up to ∼ 0.4 for f dg = 10 −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).
6. 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.
7. An extension of this stability analysis to regions in our disk models left out of the simulation domain shows that, for sufficient depletion of small grains via dust coagulation and sufficiently high molecular emissivity, protoplanetary disks can be VSI-unstable in surface layers up to several tens of au while remaining stable at the midplane.This picture is consistent with current observations of vertically extended edge-on disks at optical and near-infrared wavelengths (e.g., Villenave et al. 2020) appearing as thin layers in ALMA dust continuum images (e.g., Villenave et al. 2022), evidencing a low level of turbulent dust stirring at the midplane.The stability of the surface layers is subject to the gas molecular composition and their local temperature.
If a baroclinic amplification of meridional eddies does indeed occur in protoplanetary disks, it may affect dust coagulation and fragmentation by increasing the collisional velocities of small dust grains, affecting their diffusivity, and potentially concentrating them (Klahr & Henning 1997).On top of this, this process could play a role in the birth of large-scale azimuthal vortices.It is however unclear whether bands of constant specific angular momentum can form in real disks, as these are unstable to non-axisymmetric modes which cannot be captured by our axisymmetric simulations.However, reduced-∇j z bands seem to be present in a number of 3D simulations (Richard et al. 2016;Manger & Klahr 2018;Flock et al. 2020).Since the formation of such regions is a requirement for the formation of amplified eddies, determining whether this process can occur and, in that case, how it may interact with non-axisymmetric modes requires rather expensive high-resolution 3D simulations.These would make it possible to explore the long-term evolution of the VSI and its induced secondary instabilities without running into artifacts produced in axisymmetric simulations, such as the seemingly unrestricted growth of large uniform-j z regions.
Further work is also needed to make better predictions on the stability of disk surface layers.In particular, improvements to this work can be made by modeling the distribution and temperature of small dust grains, which determine the collisional timescale and therefore affect the strength of the VSI.Conversely, the balance between vertical stirring and settling of dust grains depends on the level of turbulence, and therefore a self-consistent model is needed.Such models should also consider the likely VSI suppression in strongly ionized regions, and, conversely, the expected prevalence of the VSI in regions where nonideal MHD effects become important (e.g., Latter & Kunz 2022;Cui & Bai 2022).
Appendix C: Thermal relaxation timescale

C.1. Dust-gas collisional timescales
In general, the relaxation timescale t rel for gas temperature perturbations differs from t cool if the timescale of energy transfer between dust and gas particles is long enough that the dust and gas temperatures, T d and T g respectively, can be different during the relaxation.To consider this effect, we obtain expressions for t rel following an analogous procedure to the analysis in Barranco et al. (2018).As in Paper I (Appendix B), we neglect advection by taking v = 0. We also neglect gas opacities and assume the heating rate of spherical dust grains due to collisions with gas particles to be given, as in Burke & Hollenbach (1983), by where a is the grain radius, n g = ρ/(µu) and v g = 8k B T πµu are the number density and the mean thermal speed of gas particles, respectively, while A 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 I[f (a)] = amax amin da f (a)n(a), 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 ρ tot d (i.e., not the density of small grains ρ d used for the computation of opacities) can be computed as I 4πa 3 3 ρ gr , where ρ gr is the bulk density of the dust grains, and thus I[πa 2 ] = 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/t rel , with || /(t g t r ) −1 , (C.5) where t −1 || = t −1 g + t −1 d + t −1 r .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 t rel ≈ t g (1 + t r /t d ) = t g + t cool,thin , (C.6) where t cool,thin is the optically thin limit of the radiative cooling time (Equations 16 and 17 in Paper I).This approximation was based on the disk parameters in that work, where grain sizes in the range [1, 100] cm were considered, whereas in the present work we assume the much smaller a max = 19 µm and 1.8 mm for f dg = 10 −3 and 10 −4 , respectively.This is a more reasonable assumption in the regions where collisional times become relevant, where we can expect centimeter-sized grains to be removed by vertical settling, as shown for example in the model by Flock et al. (2020).Yet, with our employed size distribution, the relations t r , t d ≪ t g are still satisfied, the latter of which we have verified taking the values of the specific heats capacities of graphite and amorphous silicates from Butland & Maddison (1973) and Zeller & Pohl (1971), respectively.Therefore, as we have verified in our disk models, the approximation in Eq. (C.6) is still valid with 0.1% and 0.01% accuracy for f dg = 10 −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 ρ tot d /ρ = 10 −2 , 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 t g ∝ √ a max a min /ρ tot d .Since ρ tot d ∝ √ a max − √ a min and a min ≪ a max , this results in t g ∝ √ a max a min /( √ a max − √ a min ) ≈ √ a min , and thus t g is rather insensitive to the precise value of a max .

C.2. Inclusion of gas radiative emission
The same analysis can be applied if the emission of radiation by gas molecules is taken into account.We assume for this that the gas emits thermally at a rate determined by the Planck-averaged gas absorption opacity κ g P (T ) (Kirchhoff's law), and so we include a gas-radiation interaction terms of the same form of G 0 , specifically G 0 g = κ g P (T )ρ(E r − a R T 4 ).This leads to the following system of equations: (C.8)

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

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

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

Fig. 6 .
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.

Fig. 7 .
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 Sections 3.3 and 5.3 are shown for comparison.
equation relating the vertical shear to the disk's baroclinicity; see Equation 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 (Equation (2)).As soon as uniform-j z bands start to form, the term τ c ∝ ∂ z j 2 z 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 ) = sgn(z), the amplified vortices in our simulations rotate clockwise for z > 0 and counterclockwise for z < 0.
Upwardand downward-directed flows are eventually connected at the upper disk layers, either due to the angular momentum excess of the parcels moving away from the midplane or due to the rotation produced by the baroclinic torque once |τ b | > |τ c |.This rotation further enhances the angular momentum exchange between the connected regions, until bands of uniform j z are formed, inside of which |τ b | ≫ |τ c |.

Fig. 10 .
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 Equation (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.

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

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

Fig. 16 .
Fig. 16.Same as Figure 15 but locally computing tcrit as in the right-hand side of the global stability criterion (Equation (8)).
Fig. A.2 in Paper I), as the shear between adjacent bands increases and the KH eddies are better resolved.Moreover, larger resolutions result in smaller average VSI wavelengths (Fig. C.3 in Paper I), and therefore also in a larger density of KH-unstable interfaces.All of these effects may contribute to the obtained reduction of the saturated average kinetic energy and Reynolds stresses with increasing resolution shown in Paper I.
max a min ) −1/2 .The resulting equations for the evolution of the temperature of gas and dust particles arec g ρ ∂ t T g = −Λ dg c d ρ tot d ∂ t T d = Λ dg + c G 0 − ∇ • F Irr , (C.3) where c g = k B(Γ−1)µu 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 (Equation5in Paper I) with a R T 4 d .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:∂ t δT g = −(δT g − δT d )/t g ∂ t δT d = −(δT d − δT g )/t d − δT d /t r , RT 4 +b P (a R T 4 −Er)] , with b P = d log κ d P d log T .

Table 1 .
Varying parameters of the simulations analyzed in this work.