Vertical shear instability in two-moment radiation-hydrodynamical simulations of irradiated protoplanetary disks I. Angular momentum transport and turbulent heating

Context. Hydrodynamical instabilities are likely the main source of turbulence in weakly ionized regions of protoplanetary disks. Among these, the vertical shear instability (VSI) stands out as a rather robust mechanism due to its few requirements to operate, namely a baroclinic stratiﬁcation, which is enforced by the balance of stellar heating and radiative cooling, and short thermal relaxation timescales. Aims. Our goal is to characterize the transport of angular momentum and the turbulent heating produced by the nonlinear evolution of the VSI in axisymmetric models of disks around T Tauri stars, considering varying degrees of depletion of small dust grains resulting from dust coagulation. We also explore the local applicability of both local and global VSI-stability criteria. Methods. We modeled the gas-dust mixture in our disk models by means of high-resolution axisymmetric radiation-hydrodynamical simulations including stellar irradiation with frequency-dependent opacities. This is the ﬁrst study of this instability to rely on two-moment radiative transfer methods. Not only are these able to handle transport in both the optically thin and thick limits, but also they can be integrated via implicit-explicit methods, thus avoiding the employment of expensive global matrix solvers. This is done at the cost of artiﬁcially reducing the speed of light, which, as we veriﬁed for this work, does not introduce unphysical phenomena. Results. Given su ﬃ cient depletion of small grains (with a dust-to-gas mass ratio of 10% of our nominal value of 10 − 3 for < 0 . 25 µ m grains), the VSI can operate in surface disk layers while being inactive close to the midplane, resulting in a suppression of the VSI body modes. The VSI reduces the initial vertical shear in bands of approximately uniform speciﬁc angular momentum, whose formation is likely favored by the enforced axisymmetry. Similarities with Reynolds stresses and angular momentum distributions in 3D simulations suggest that the VSI-induced angular momentum mixing in the radial direction may be predominantly axisymmetric. The stability regions in our models are well explained by local stability criteria, while the employment of global criteria is still justiﬁable up to a few scale heights above the midplane, at least as long as VSI modes are radially optically thin. Turbulent heating produces only marginal temperature increases of at most 0 . 1% and 0 . 01% in the nominal and dust-depleted models, respectively, peaking at a few (approximately three) scale heights above the midplane. We conclude that it is unlikely that the VSI can, in general, lead to any signiﬁcant temperature increase since that would either require it to e ﬃ ciently operate in largely optically thick disk regions or to produce larger levels of turbulence than predicted by models of passive irradiated disks.


Introduction
Understanding the diverse mechanisms that produce turbulence in protoplanetary disks is a necessary step toward improving our understanding of planet formation.Turbulence directly impacts the transport, settling, coagulation, and fragmentation of dust grains by controlling their diffusivity and collisional velocities (Ormel & Cuzzi 2007;Johansen et al. 2014;Klahr & Schreiber 2020), thus playing an important role in planetesimal formation.It can also contribute to stellar accretion by producing outward transport of angular momentum, although it is unclear whether the low level of transport produced by different hydrodynamical (HD) and magnetohydrodynamical (MHD) instabilities (e.g., Lesur et al. 2023) can explain current observations of accretion rates (Hartmann et al. 1998;Manara et al. 2016), or whether instead accretion is mainly driven by magnetized winds (e.g., Bai & Stone 2013;Béthune et al. 2017;Gressel et al. 2020).On the other hand, large-scale structures produced by different instabilities, such as zonal flows and vortices, are able to concentrate dust grains, possibly triggering the streaming instability and leading to the formation of planetesimals (Johansen et al. 2007;Gerbig et al. 2020;Schäfer et al. 2020).In poorly ionized "dead" zones of protoplanetary disks, where the magnetorotational instability (MRI, e.g., Balbus & Hawley 1991) is expected to be suppressed by nonideal phenomena such as Ohmic and ambipolar diffusion and, depending on the orientation of large-scale magnetic fields, the Hall effect (Bai 2015;Simon et al. 2018), purely HD instabilities have been predicted to take over as the main mechanisms driving turbulence (e.g., Lyra & Umurhan 2019;Pfeil & Klahr 2019).Among these, the vertical shear instability (VSI) stands out as a rather robust phenomenon dominating the production of turbulence in numerous HD (Nelson et al. 2013;Manger & Klahr 2018;Flock et al. 2017b) and even nonideal MHD simulations (e.g., Cui & Bai 2022).
Local and vertically global stability analyses (e.g., Urpin & Brandenburg 1998;Urpin 2003;Lin & Youdin 2015) predict that the VSI can be triggered as long as two conditions are met.On the one hand, the disk's pressure and density stratification must be baroclinic (∇p×∇ρ ̸ = 0, where p and ρ are the gas pressure and density, respectively), which in vertical and radial hydrostatic equilibrium naturally leads to an azimuthal velocity stratification varying with distance from the midplane, that is, with vertical shear.Baroclinic stratifications naturally result from temperature gradients either in the radial or in the vertical direction, which are typically enforced by the equilibrium between stellar or turbulent heating and radiative cooling.On the other hand, buoyancy forces, which can restore gas perturbations into equilibrium, must be suppressed.Under these conditions, a gas parcel moving in a direction between its corresponding constant specific angular momentum surface and the disk rotation axis is accelerated away from is original position due to its angular momentum excess (e.g., Knobloch & Spruit 1982).Globally, this process can destabilize radially traveling inertial waves (Nelson et al. 2013;Barker & Latter 2015), leading to the production of vertically elongated vertical flows reaching significant fractions of the sound speed (e.g., Flock et al. 2017b).
If the additional buoyancy produced by magnetic tension (e.g., Latter & Papaloizou 2018) and dust backreaction (Lin 2019) are negligible, the stabilizing effect of buoyancy can be suppressed for short enough thermal relaxation timescales.A local constraint on the thermal relaxation timescale leading to instability was heuristically derived in Urpin (2003), following the work by Townsend (1958), in the spirit of a Richardson number (e.g., Chandrasekhar 1961), quantifying the ratio between the work of restoring forces and the energy that can be extracted from the vertical shear.Local stability analyses are however insufficient to describe both the linear and nonlinear global evolution of the instability, which produces flows spanning different local conditions over several scale heights in the vertical direction.To overcome this limitation, a global stability criterion was derived in Lin & Youdin (2015) for disks with vertically uniform temperatures and cooling rates.Even though global HD simulations support this criterion in both the linear and nonlinear regimes (e.g., Nelson et al. 2013;Manger et al. 2021), this analysis cannot account for deviations from the cited assumptions.In such cases, local stability criteria may still be useful to determine the stability of specific regions, but little progress has been made in this direction (see, however, Stoll & Kley 2016;Pfeil & Klahr 2021;Fukuhara et al. 2023).
If dust and gas species are kept in thermal equilibrium by collisions, the thermal relaxation timescale is determined by the radiative cooling rate.Computing this quantity is in general not an easy task, as it depends on the optical depth along the typical wavelengths of the perturbations produced by the VSI (see, e.g., Lin & Youdin 2015).These are cur-rently unknown due to a lack of theoretical predictions, together with the fact that in HD computations the minimum wavelengths are limited by grid diffusion.Even if characteristic VSI modes are radially optically thin, in which case the cooling rate becomes scale-independent, the temperature perturbations produced by the VSI flows may span several characteristic sizes with different optical depths exceeding the typical VSI wavelengths.On top of this, if larger-scale structures such as vortices are produced by secondary instabilities or independent mechanisms (Marcus et al. 2015;Barranco et al. 2018;Manger & Klahr 2018;Pfeil & Klahr 2021), their optical depths may largely exceed that of the VSI modes, and hence they may cool at a different rate.Thus, the issue of computing accurate cooling timescales in HD simulations is only avoided if these naturally result from an effective solution of the radiative transfer equation, as it is done in radiation hydrodynamics (Rad-HD) codes.Moreover, unlike for isothermal and β-cooling prescriptions, such codes can be used to quantify the turbulent heating produced in VSI-unstable disk regions, which is particularly relevant for planet formation models relying on α-viscous heating prescriptions (e.g., Burn et al. 2022;Voelkel et al. 2022).
In this work we investigated the nonlinear evolution of the VSI in axisymmetric Rad-HD simulations of protoplanetary disks including heating by stellar irradiation with realistic tabulated opacities, focusing on the VSI growth and saturation phases, the produced transport of angular momentum and heating, and the applicability of local stability criteria.Both the baroclinic stratification and the thermal relaxation timescales required to trigger the VSI in our simulations arise from the balance of irradiation heating and radiative cooling.We studied a small disk region between 4 and 7 au around a T Tauri star and consider different degrees of depletion of small dust grains due to coagulation and settling, making the disk either vertically thick or thin and resulting in the formation different stability regions.
We neglected departures from thermal equilibrium between dust and gas species occurring in low-density regions, where the thermal relaxation timescale depends on a combination of the radiative cooling timescale and the collisional timescales of dust and gas species (see, e.g., Malygin et al. 2017;Barranco et al. 2018).This effect is considered in an extension of our stability analysis presented in a second part of this work (Melon Fuksman et al. 2023, Paper II henceforth).
We employed resolutions of up to ∼ 200 cells per scale height, which are unprecedented in Rad-HD VSI simulations, aiming to investigate the role of secondary instabilities on the VSI saturation in Paper II.The resulting computational overhead is limited by artificially reducing the speed of light, which in this context can be done as long as thermal relaxation timescales are unaffected by this operation, as we verified both analytically and numerically in this work.Also advantageous in terms of computational efficiency is the good parallel scalability of our employed numerical code (Melon Fuksman et al. 2021), which is the first Rad-HD method applied to the study of disk instabilities not to rely on globally implicit solvers.
This article is organized as follows.In Section 2 we outline our employed numerical scheme and describe our disk models.In Section 3 we characterize the velocity and vertical shear distributions in the different growth and saturation phases of the VSI, while in Section 4 we analyze the resulting redistribution of angular momentum in relation with the small and large-scale flows in our simulations.In Section 5 we explore the thermal evolution of our simulated disks, and in Section 6 we analyze our obtained stability regions in terms of local criteria.Finally, in Section 7 we discuss our results, and in Section 8 we summarize our conclusions.Additional figures and calculations are included in the appendices.

Radiation hydrodynamics
We modeled the gas-dust mixture in a circumstellar disk and its heating and cooling processes by means of the Rad-HD module (Melon Fuksman et al. 2021) implemented in the open-source PLUTO code (version 4.4, Mignone et al. 2007).In the chosen configuration for this work, the following system of equations is solved: where ρ, p, and v are the gas density, pressure and velocity, while E r , F r , and P r are respectively the radiation energy, flux, and pressure tensor.We normalize the radiation flux by the speed of light c in such a way that both E r and F r are measured in energy density units.The gas energy density is defined as E = ρϵ + 1 2 ρv 2 imposing an ideal equation of state to the internal energy density, namely ρϵ = p Γ−1 , where we set the heat capacity ratio as Γ = 1.41 (Decampli et al. 1978).This system is completely defined by imposing the M1 closure (Levermore 1984) to define the radiation pressure tensor as a function of F r and E r as , and δ ij is the Kronecker delta.
We assume axisymmetry and solved these equations in spherical coordinates (r, θ), where θ is the polar angle, solving also for the azimuthal ϕ-component of all vector fields.We include a gravitational potential modeling a star of mass M s located at the center of coordinates, that is, Φ = −M s G/r.We include stellar heating by computing after each HD step the stellar irradiation flux produced by the central star as where T s is the star temperature, R s is the star radius, B ν (T ) = (2hν 3 /c 2 )/(e hν/k B T − 1) is the Planck radiative intensity at a temperature T and frequency ν, and h and k B are the Planck and Boltzmann constants, respectively.The gas temperature, assumed to be the same for dust grains, is computed according to the ideal gas law T = µu kB p ρ , where µ = 2.35 and u are respectively the assumed gas mean molecular weight and the atomic mass unit.The frequencydependent optical depth τ ν (r, θ) is computed via ray tracing after each density update.We include for this calculation the optical depth of the part of the disk left out of the simulation domain.When constructing initial conditions (Section 2.2), we estimate this value following the prescription in Melon Fuksman & Klahr (2022).In the Rad-HD simulations, we retrieve it from the initial conditions.
Optical depths are computed considering only opacities due to dust absorption, which at optical and near-infrared wavelengths is dominated by small (≲ 1 µm) grains.This does not mean that there is not significant mass in larger grain sizes, but that larger sizes do not contribute to the total opacity, as specified in the next section.We neglect dust settling of small grains and assumed a perfect coupling between the dust and gas distributions with a constant dust-to-gas mass ratio f dg of small grains, thus computing the dust density of small grains as ρ d = f dg ρ.We employ for this calculation the same absorption opacities as in Krieger & Wolf (2020, 2022), which are obtained for spherical grains composed of 62.5% silicate and 37.5% graphite with a density of 2.5 g cm −3 and a size distribution dn ∼ a −3.5 da with radii in the range a ∈ [5,250] nm.Opacity values are logarithmically sampled in the frequency range [ν min , ν max ] = [1.5 × 10 11 , 6 × 10 15 ] Hz.We also use these opacities to compute the frequency-integrated radiation-matter interaction terms, which in the gas comoving frame have the expressions where tildes indicate comoving quantities, while a r = 4σ SB /c is the radiation constant, σ SB is the Stefan-Boltzmann constant, and κ d P and χ d R are, respectively, the Planck-and Rosseland-averaged absorption and extinction coefficients.The laboratory frame coefficients (G 0 , G) are obtained via a Lorentz transformation of ( G0 , G), which introduces terms of orders v/c and v 2 /c 2 as detailed in Melon Fuksman et al. (2021).In the simulations presented in this work these terms are negligible, as is the contribution of G to the gas momentum, which has no role in the resulting gas dynamics.Throughout this work we neglect scattering, namely, we take χ d R = κ d R , which is the Rosselandaveraged absorption opacity, although the derivations in Appendix B consider a general χ d R .Graphs of the frequencydependent and frequency-averaged opacities used in this work are shown in Melon Fuksman & Klahr (2022).
Our resolution method relies on the reduced speed of light approximation (RSLA), in which the speed of light in the radiation block is reduced to an artificially low value ĉ < c in order to reduce the scale disparity between gas and radiation characteristic speeds.This avoids the employment of the prohibitively small time steps that would be given by the Courant-Friedrichs-Lewy condition using the real speed of light, thus reducing the overall computational cost.The obtained solutions of the Rad-HD equations are in general unchanged for large enough ĉ values, as we later verify in our setup.).Dashed, white lines indicate the τP (Ts) = 1 surface for stellar photons computed using Planck-averaged opacities, solid white lines show the τP = 1 surfaces for vertically traveling infrared photons, and green dash-dotted lines correspond to the z = ±H surfaces.

Initial conditions
We considered disks initially in hydrostatic equilibrium, with temperature distributions resulting from the balance of heating due to stellar irradiation and cooling by infrared emission.We obtained our initial configurations by means of the same iterative procedure in Melon Fuksman & Klahr (2022), in which we alternated the computation of (ρ, v ϕ ) solving hydrostatic equilibrium equations and of (T, E r , F r ) by means of our treatment of stellar and reprocessed radiation.The latter step was carried out by evolving the radiation fields in each iteration for a total time t iter of 0.1 orbits at 1 au, which is enough to prevent artefacts induced by self-shadowing (Melon Fuksman & Klahr 2022) and reach convergence.
We constructed models of protoplanetary disks around a T Tauri star of mass M s = 0.5 M ⊙ , effective temperature T s = 4000 K, and radius R s = 2.5 R ⊙ , which leads to an effective luminosity of ∼ 1.4 L ⊙ .We assumed a gas column density of Σ(r) = 600 g cm −2 (r/au) −1 which is kept constant throughout iterations.The hydrostatic and thermal equilibrium computations were carried out in a domain (r, θ) ∈ [0.4,100] au × [π/2 − 0.5, π/2 + 0.5] discretized with a resolution of N r × N θ = 240 × 400.The total gas mass in that region is thus ∼ 0.04M ⊙ .If a vertically integrated dust-to-gas mass ratio of approximately 0.01 is assumed, the total dust mass is approximately on the upper end of the distribution of estimated dust masses as a function of stellar masses for disks observed with ALMA (Manara et al. 2018).
The only physical parameter that varies among our simulations is the assumed dust-to-gas ratio of small grains, f dg , whereas all other physical quantities are kept constant.We considered two values of f dg : a nominal case with f dg = 10 −3 and a dust-depleted case with f dg = 10 −4 .If we assume a total dust-to-gas mass ratio of 10 −2 including all grain sizes, these values correspond respectively to maximum grain sizes of 19 µm and 1.8 mm, in which cases only 10% and 1% of the total dust mass is in the small-size pop-ulation responsible for the heating and cooling processes in the disk.The latter case can represent the late evolution of the disk after a few Myr (Birnstiel et al. 2012) or simply a disk where significant coagulation and settling took place.As shown by the synthetic images in Dullemond et al. (2022), that degree of dust depletion can still be consistent with the vertically extended appearance of scattered light images of protoplanetary disks (see also Pfeil et al. 2023).

Hydrostatic distributions
In our Rad-HD simulations, we employed a higherresolution grid than the one used to construct initial conditions.We defined this new grid in the domain of study (r, θ) ∈ [4, 7] au × [π/2 − 0.3, π/2 + 0.3], in which we retrieved all initial distributions from the hydrostatic solution via bilinear interpolation.Inner regions of protoplanetary disks (≲ 10 au) are rather unexplored in numerical simulations (see, however, Stoll & Kley 2016;Pfeil & Klahr 2021) as they are currently hard to resolve in mm observations (e.g., Andrews et al. 2018).However, our chosen disk region serves as a laboratory to study the thermal evolution produced by the VSI (Section 5) in regions of different vertical optical depth, as this domain becomes either vertically optically thick or thin depending on the dust content (Fig. 1 and this section).On the other hand, its reduced radial extent also allows us to resolve the secondary instabilities studied in Paper II, where we also investigate the stability of outer disk regions via local criteria.
The resulting initial temperature distributions for both f dg values are shown in Fig. 1, while vertical slices at r = 5.5 au are shown in Fig. 2. Close to the midplane, the disk is approximately vertically isothermal up to a few scale heights above the midplane (∼ 2.5 and 3.5 for f dg = 10 −3 and 10 −4 , respectively), with a temperature well described by a power law T ∝ R −q , where q = 0.49 and 0.48 for f dg = 10 −3 and 10 −4 , respectively, where R is the cylindrical radius.This corresponds to a disk aspect ratio H/R = h 0 R 4 au ) where h 0 = 0.054 and 0.052 for f dg = 10 −3 and 10 −4 , respectively, while H is the pressure scale height measured at the midplane.The flaring index, 1−q 2 , is in every case positive and close to 1/4.
At higher distances from the midplane, the temperature distributions transition from the midplane values of 65 − 90 K to higher values up to ∼ 250 K.This transition occurs near the location of the τ ν = 1 surfaces for photons emitted at the star with frequencies corresponding to maximum stellar luminosity, where the irradiation heating rate is maximal.As shown in Fig. 1, this region is approximately located at the irradiation surface, which is defined through the condition τ P (T s ) = 1, where τ P (T s ) is the same optical depth computed with the Planck-averaged opacity at the star temperature.Throughout this work, we shall use the terms "upper" or "surface layers" and "middle layer" to refer to the disk regions above and below that surface, respectively.In the same figure it is shown that this temperature transition occurs a few scale heights above the region where the disk becomes optically thin for vertically trav-eling reprocessed photons, defined as the vertical τ P = 1 surfaces computed with κ d P (T ) (Fig. 1).In our models, the total vertical optical depth computed in this way from z = 0 to infinity is τ P ≈ 10 − 20 for f dg = 10 −3 and 0.7 − 2 for f dg = 10 −4 , in which case most of the domain is optically thin for vertically traveling photons.
Above the irradiation surface, the disk is optically thin to radially traveling photons, and thus the temperature results from the balance of irradiation heating and optically thin cooling.Specifically, in terms of Equation (1), which is roughly proportional to r −1/2 if opacities vary slowly with temperature and τ P (T s ) ≪ 1.This explains the obtained midplane temperature power-law index q ∼ 0.5 (steeper profiles would require the inclusion of an additional heat source, for instance, due to viscosity).Under this surface, the temperature is determined by the transport of reprocessed radiation.As described in Melon Fuksman & Klahr (2022), the fact that we are using a fluid approach to describe this process causes a slight deviation from isothermality right below the temperature transition (z/r ∼ 0.1 − 0.2 in Fig. 2), which barely alters the resulting shear distribution with respect to an isothermal stratification (see Fig. 2).

Vertical shear distribution
When the disk is in hydrostatic equilibrium, an expression for the vertical shear can be obtained by writing the curl of the momentum equation (second line in Equation ( 1)) in cylindrical coordinates (R, ϕ, z), which leads to the following relations: where , with an analogous expression for the epicyclic frequency z , quantifies the vertical variation of the specific angular momentum j z = R 2 Ω, while Ω = v ϕ /R is the orbital angular frequency and T = p/ρ.This equation, which is a form of the so-called thermal wind equation in geophysical scenarios, shows that a vertically varying rotation velocity occurs if the disk is baroclinic, namely, such that ∇p × ∇ρ ̸ = 0.As evidenced by the last line, baroclinic stratifications naturally occur in vertically isothermal disks as long as ∂ R T ̸ = 0.For this reason, below the irradiation surface, κ 2 z approximates with good accuracy its value for a vertically isothermal temperature distribution, as shown in Fig. 2. On the contrary, above the irradiation surface, the term proportional to ∂ z T becomes important, and the obtained vertical shear results from the balance of both terms on the right-hand side of Eq. ( 5), which have opposite signs.In particular, ∂ z T ∂ R log ρ is positive above the midplane and negative below, and the opposite is true for −∂ R T ∂ z log ρ.In every case, the term proportional to the radial temperature gradient dominates, and its effect is diminished by the opposite contribution of the vertical temperature gradient.At the region of maximal shear, located at the temperature transition region, both |∂ z T ∂ R log ρ| and |∂ R T ∂ z log ρ| are about 5 times larger than the maximum |κ 2 z |, which reaches up to 2 times its value for a vertically isothermal disk.
We can expect this larger vertical shear than isothermal disks to be translated into larger local linear growth rates and an overall stronger instability, as suggested by the relation α ∝ q 2 in Manger et al. (2021).For instant cooling, a local Boussinesq linear analysis predicts that the local growth rate of the fastest-growing mode is & Brandenburg 1998;Klahr 2023).The 2D distribution of this quantity in our models is shown in Fig. 1, where it can be seen that the fastest growth is expected above the irradiation surface.This does not necessarily mean larger unstable regions than in isothermal disks, since these not only depend on the vertical shear but also on the stabilizing vertical buoyancy, which is also affected by the temperature transition.We can quantify this in terms of the squared vertical Brunt-Väisälä frequency corresponding to vertical buoyant oscillations, (e.g., Rüdiger et al. 2002).This quantity is larger than in isothermal disks at the region of maximal shear and smaller above (see Fig. 2), due to the larger and smaller vertical entropy gradient in these locations, respectively.The combined effect of these quantities on the disk stability is studied in Section 6.On the other hand, local analyses such as this are not enough to predict the nonlinear evolution of the instability, which is essentially nonlocal due to the vertical mass advection along several scale heights and the cooling due to radiative transfer over a broad range of lengthscales.
The evolution and saturation of the VSI in its nonlinear regime are studied in the next section.

Rad-HD simulations' setup
Starting from the described initial conditions, we ran Rad-HD simulations in order to study the growth and saturation of the VSI.A list of the presented runs and their corresponding parameters is shown in Table 1.Simulations are labeled as, dgN dg cN c _N θ where the numbers N dg and N c indicate the employed f dg and ĉ/c values, respectively, while N θ is the resolution in the θ-direction.We employed four different resolutions, N r × N θ = 240 × 256, 480 × 512, 960 × 1024, and 1920 × 2048, which correspond to approximately 25, 50, 100 and 200 cells per scale height, respectively.This discretization is chosen in such a way that in all cases the cell aspect ratio r∆θ/∆r is almost exactly 1.
In every run considered in the next sections we used ĉ/c = 10 −4 .This value is still high enough to avoid introducing unphysical phenomena in our simulations, as justified in Appendices B and C via analytical computations and testing with varying ĉ/c, respectively.In particular, it is shown there that the radiative cooling time is unaffected for high enough ĉ/c, which is crucial for this work given the dependence on the VSI on that timescale.
The Rad-HD equations are solved using third-order Runge-Kutta time integration for the gas fields and the IMEX method derived from it in Melon Fuksman & Mignone (2019) for the radiation fields.Gas and radiation fluxes are computed using the HLLC solvers in Toro (2009) and Melon Fuksman & Mignone (2019), respectively.We employed in every case third-order WENO spatial reconstruction (Yamaleev & Carpenter 2009) with the modifications introduced in Mignone (2014) for spherical grids.In all cases, Dirichlet boundary conditions are imposed by fixing both radiation and HD fields in the ghost cells to their corresponding values in the initial hydrostatic configurations.

Growth phases
We begin by characterizing the flows produced by the VSI in its different growth stages.In all of our simulations, vertically elongated VSI modes can be first seen above the irradiation surface beginning at the smallest radii and occupying the entire radial domain within ∼ 10 orbits, after which they propagate toward the midplane and occupy the entire domain.This time sequence holds some similarity with the occurrence of surface modes in the vertically isothermal disk model by Barker & Latter (2015).Such modes, which dominate the early growth phases in that work, have maximum amplitude close to the vertical boundaries, where the local growth rate predicted by linear theory is maximum (Urpin & Brandenburg 1998;Klahr 2023).In our disk models, the local growth rate is instead maximum right above the irradiation surface (Fig. 1), and thus VSI modes first appear in those locations before occupying the entire irradiated regions.Thus, both phenomena might result from the time ordering determined by the local growth rate distribution.
The resulting vertical gas velocity and mass flux in our highest-resolution simulations are shown in Figure 3 after 50 and 300 orbits (we define one orbit as the rotation period T 5.5 at 5.5 au).For f dg = 10 −3 , as shown in the velocity snapshots at t = 50 T 5.5 , the instability first produces "finger" modes (Nelson et al. 2013) characterized by an approximately vertically antisymmetric v z , in which case no mass is transported through the midplane.After ∼ 100 orbits, the velocity distribution transitions into vertically symmetric "body" modes occupying the entirety of the vertical domain, in which case mass is transported upward and downward at the midplane.On the contrary, for f dg = 10 −4 we observe no transition between finger and body modes.Unlike in the previous case, the VSI is almost entirely suppressed at the middle layer, where meridional mass fluxes are ∼ 10 2 times lower than for f dg = 10 −3 .A localized VSI suppression at the midplane is also observed in the β-cooling simulations by Fukuhara et al. (2023).As later described in Section 6, this reduced activity results from the stabilizing effect of buoyancy for long enough cooling times.
Gas flows are in every case fastest at the disk upper layers, where the vertical shear is maximum and the local instability criterion is met in all cases (Section 6).The gas r.m.s.velocity in the saturated phase is maximum above the irradiation surface (see Fig. 4), reaching at the highest resolution up to ∼ 50 and 40 m s −1 , corresponding to Mach numbers M ≈ 0.07 and 0.05, for f dg = 10 −3 and 10 −4 , respectively.At that resolution, the VSI modes produce maximum vertical velocities of ∼ 200 m s −1 (M ∼ 0.24) and ∼ 100 m s −1 (M ∼ 0.12) for f dg = 10 −3 and 10 −4 , re- Notes.Dust-to-gas mass ratio of small grains (f dg ), resolution (Nr × N θ ), approximate number of cells per scale height (H/∆r), ratio between the reduced and real values of the speed of light (ĉ/c), 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 cooling time Ωt cool and the ratio t cool t crit between the cooling time and the local critical value for instability (Equation ( 14)) computed in Section 6, both shown at the maximum height in the domain (z/r ≈ 0.3) and the height z = ΓH 2 determining the vertically global stability criterion in Equation ( 15).
spectively.At the middle layer, r.m.s.velocities are about 3 times lower than at the upper layers for f dg = 10 −3 , whereas for f dg = 10 −4 a reduction of 2 orders of magnitude occurs due to the suppression of the VSI in that region.Some level of convergence can be seen in the r.m.s.velocities at different resolutions in Fig. 4, but this is not true for other quantities.As discussed in the next sections, significant changes in the gas flow are seen when increasing the resolution, and likely our highest resolution is not enough to reach convergence.In particular, the width of the body modes consistently decreases with increasing resolution, as shown in Fig The different growth phases of the instability can be determined by evaluating the time evolution of the volumeaveraged meridional kinetic energy, 5 for varying dust content and resolution.For f dg = 10 −3 , two distinct growth phases can be seen: after ∼ 50 − 100 orbits, the kinetic energy plateaus following the saturation of the finger modes, after which it starts increasing again as body modes start growing and it stalls when these saturate.Similar behavior is seen, for example, in Stoll & Kley (2014), Manger et al. (2021), andFlores-Rivera et al. (2020).The saturation of the body modes occurs at t ∼ 125 − 160 orbits in all cases except for N θ = 256, in which case body modes take up to ∼ 400 orbits to grow due to the higher numerical diffusion.On the contrary, for f dg = 10 −4 , we observe a single growth phase, as no clear transition occurs between dominant modes due to the reduced VSI activity close to the midplane.
It can be seen in Fig. 5 that the saturation value of ⟨ϵ k ⟩ is smaller for increasing resolution for f dg = 10 −3 .This result highlights a general feature of global VSI models discussed, for example, in Barker & Latter (2015), which is the fact that the smallest scales in which the VSI can develop are limited by the viscous length (Latter & Papaloizou 2018), which in inviscid models such as ours is given entirely by numerical diffusion.On the contrary, the saturated value of ⟨ϵ k ⟩ for f dg = 10 −4 reaches convergence for N θ ≥ 512 after ∼ 30 orbits, when the saturation of the VSI modes occurs above the irradiation surface.However, changes in the gas flow can still be seen for increasing resolution, not only in the wavelength of the body modes but also in the increased growth of parasitic Kelvin-Helmholtz (KH) eddies, visible in Fig. 3 as ripples in between the vertical flows and further studied in Paper II.As discussed in that work, the saturation of the VSI may result from the transfer of kinetic energy from the VSI flows into KH eddies.This is consistent with the earlier VSI saturation for increasing resolution, and may also explain the decreasing saturated value of ⟨ϵ k ⟩ with resolution for f dg = 10 −3 if the Kelvin-Helmholtz instability (KHI) limits the maximum energy of the saturated body modes.This effect may be enhanced by the favored growth of the KHI when numerical diffusion is reduced, and also by the increase with resolution of the number of KH-unstable interfaces between vertical flows, given that larger resolutions lead to smaller VSI wavelengths (Figs.A.2 and C.3).On top of this, some KH eddies are subject to baroclinic amplification, as also analyzed in Paper II.The fast vortices produced by this process are responsible for the maximum velocities in Figure 5, which explains why the maximum velocities increase with resolution despite the overall decrease of the kinetic energy.

Vertical shear evolution
As the VSI modes grow, the produced vertical flows tend to remove the initially unstable velocity stratification that originated the instability.Even though the initial baroclinic state of the disk is maintained by the balance of stellar irradiation and radiative cooling, the vertical shear is no longer determined by Equation ( 5) after a few orbits, since it is affected by the departure of the gas velocity from its hydrostatic equilibrium distribution.A comparison of the vertical shear distribution in the initial state and after 300 orbits is shown in Fig. 6 next to the corresponding distributions at that time of the specific angular momentum, j z = R 2 Ω.In these snapshots it can be seen that the VSI flows generate regions of approximately uniform j z , where consequently the vertical shear is reduced with respect to its initial state, as As explored in Paper II, these bands begin forming due to the angular momentum redistribution produced by the motion of VSI-destabilized gas parcels across constant-j z surfaces (see also James & Kahn 1970).The subsequent reduction of the vertical shear causes the baroclinic torque to prevail over the centrifugal torque.Once the vertically moving flows are connected in upper disk regions, angular momentum can be transported in the sense of circulation favored by the baroclinic torque, until uniform-j z regions are formed (see Paper II and Klahr et al. 2023).
The formation of reduced shear bands is harder to see for increasing resolution, due to the increased development of KH eddies (see Paper II), which produce perturbations in the shear distribution.Conversely, the radial shear is increased with respect to the quasi-Keplerian initial state (not shown here), since v ϕ = j z /R ∼ R −1/2 for a quasi-Keplerian flow and ∼ R −1 for constant j z .In each case, the shear in all directions is always maximum in between the uniform-j z zones.
Overall, the increase of the shear in between uniform-j z regions, together with the extra shear produced by KH eddies, result in an increase of the time-and volume-averaged absolute value of the vertical shear with respect to the initial state in all regions where the VSI is active, as shown in Fig. 7.For f dg = 10 −4 , this occurs everywhere except in the VSI-inactive middle region, where the shear remains unchanged.The growth of this quantity increases with resolution, almost by one order of magnitude in our highestresolution runs.This effect is likely enhanced by the favored growth of small-scale structures such as KHI eddies as numerical diffusion is reduced (Paper II).Conversely, equivalent averages of the signed ∂ z v ϕ remain unchanged with respect to the initial state.Shear fluctuations therefore largely exceed their mean value resulting from hydrostatic equilibrium a result of the nonlinear evolution of the

Reynolds stresses
We now turn to the transport of angular momentum in the disk, which depends on the time-averaged Reynolds stress produced by the VSI.Since our disk model is axisymmetric, we replaced the ϕ-averages normally employed in this calculation with time averages between 150 and 300 orbits (after the saturation of the body modes for N θ ≥ 512), taking them as representative of an ensemble average over fluctuations.We thus computed the ϕ-row of the Reynolds stress tensor as where ⟨•⟩ t indicates time average and δ denotes fluctuations with respect to the mean.The components of this tensor can be turned into stress-to-pressure ratios which can be used in effective 1D disk evolution models (e.g., Burn et al. 2022).However, as with other forms of disk turbulence (Klahr & Bodenheimer 2003), the VSI-induced angular momentum transport cannot be modeled by means of an isotropic viscosity coefficient ν = αc 2 s /Ω K (Shakura & Sunyaev 1973), since the vertical stresses produced by the VSI are orders of magnitude larger than the radial stresses, as shown in Stoll et al. (2017) and next in this section.Despite this, the R-component of S ϕ can still be used to compute the flux of vertically integrated mass and angular momentum, and therefore it can still be parameterized as an α-viscosity model implementable in 1D models.This can be seen by taking the azimuthal average of the angular momentum conservation equation as done, for example, in Balbus & Papaloizou (1999), which in axisymmetric disks can be replaced with the mentioned average over fluctuations.This results in the following expression: where ⟨•⟩ denotes the chosen averaging procedure, and where we assumed ⟨δv ϕ δρ⟩ = 0.The second and third terms in this equation account for the transport of angular momentum due to gas advection and VSI-induced mixing, respectively.In quasi steady-state accretion, the projection of the average mass flux perpendicular to the surfaces of constant specific angular momentum is thus determined by the divergence of RS ϕ .Since such surfaces are approximately vertical close to the midplane, Equation ( 7) can be used to compute the vertically integrated mass flux F Σ via an integration in z.In absence of disk winds, in which case both ⟨ρv z ⟩ and S ϕz are null at z = ±∞, this results in (see, e.g., Jacquet 2013), where Ω is evaluated at the midplane.In turn, the column density evolves according to and so the radial transport of vertically averaged mass and angular momentum in steady state depends solely on S ϕR .In Fig. 8 we show pressure-weighted radial averages of S ϕR /⟨p⟩ t and |S ϕz |/⟨p⟩ t .In every case |S ϕz |/⟨p⟩ t is between 1 and 2 orders of magnitude larger than S ϕR /⟨p⟩ t , which evidences that the mixing of angular momentum in the disks occurs predominantly in the vertical direction, as obtained in Stoll et al. (2017).In particular, S ϕz is positive for z > 0 and negative for z < 0, which means that angular momentum is transported away from the midplane.In the 2D profiles of ||S ϕ /⟨p⟩ t || in Fig. 9, we note that this quantity is maximum in between the approximately uniform-j z bands described in Section 3.2 (we can see this because some body modes do not significantly migrate radially throughout the snapshots used to compute time averages).It is then likely that the mixing of angular momentum in the interfaces between such bands, enhanced by the triggering of the KHI in those locations, is what causes them to merge and grow over time (Paper II), which also explains why the radial wavelength of the saturated modes is generally largest at the disk surface layers (see Appendix D), where the mixing is maximal and the bands can merge most efficiently.
We can see in Fig. 8 that stresses are maximum above the irradiation surface in all cases, where shear is maximum.For f dg = 10 −4 , stresses are null close to the midplane.The obtained values decrease with increasing resolution as in Stoll & Kley (2014) and Manger et al. (2021), although they show early signs of convergence for the two largest resolutions.It is however unclear whether convergence would be reached for even higher resolution, and, even in that case, whether that would still happen in 3D.The observed stress reduction with resolution is possibly a consequence of the kinetic energy loss to KH eddies (Paper II), which limits the growth of the VSI modes and therefore also their produced transport of angular momentum.
To define an effective Reynolds α value that can be translated into effective eddy viscosity models of the form S ϕR = −ν R ρR ∂Ω ∂R (i.e., the ϕR-component of the strain rate tensor of a fluid with viscosity ν R ), while reproducing the usual parameterization S ϕR = αp for thin, vertically isothermal disks, we define or equivalently, ν R = 2 3 α R c 2 s Ω −1 K (see, e.g., Armitage 2022), where Ω K (R) is the Keplerian angular velocity (the extra 2 3 factor accounts for the fact that R ∂Ω ∂R ≈ − 3 2 Ω K in thin, vertically isothermal disks).We have verified that the graphs of α R defined in this way are practically indistinguishable from those of S ϕR /⟨p⟩ t , and so the S ϕR /⟨p⟩ t averages shown in Fig. 8 can be interpreted as α R .To obtain an effective 1D viscous evolution model similar to, e.g., Lynden-Bell & Pringle (1974), we can then replace S ϕR = α R ⟨p⟩ t in Equation (8), which leads to where ⟨α R ⟩ p is the pressure-weighted z-average of α R (the corresponding expression in Jacquet (2013) has a multiplying factor of 3 instead of 2 due to their different definition of α).In these expressions we assumed a vertically uniform squared isothermal sound speed c 2 s = p/ρ, which is a valid approximation given the temperature distribution in the disk middle layer.Values of ⟨α R ⟩ p computed via interpolation along constant-R lines are shown in Fig. 10.These are in the range ∼ 10 −6 −10 −5 and ∼ 10 −7 −10 −6 for f dg = 10 −3 and f dg = 10 −4 .These values cannot be compared with other works, since most commonly the spherical α r = S ϕr /⟨p⟩ t is computed instead of the cylindrical α R .We thus show 1D averages of α r = in Fig. 11, which shows that α r ∼ 10 −5 − 10 −4 and ∼ 10 −7 − 10 −6 for f dg = 10 −3 and f dg = 10 −4 , respectively.These values are respectively about 4 and 2 times larger than the α R averages in Fig. 10.This difference despite the small z/r values is somewhat surprising, but it simply results from the anisotropy of the angular momentum transport, i.e., from the fact that |S ϕz | ≫ |S ϕR |.Thus, care must be taken in effective vertically integrated models, in which case ⟨α R ⟩ p should be used instead of the averaged α r (see Equation ( 11)).
The vertical component of the Reynolds stress can also be parameterized via an effective α z number.This value is likely representative of the dust vertical diffusivity D via D = α z H 2 Ω given that, unlike in the radial direction, turbulent mixing transports angular momentum and dust particles in the same direction (i.e., away from the midplane), following both the negative angular momentum and dust density gradients.The ϕz-component of the strain rate tensor in cylindrical coordinates is of the form S ϕz = −ν z ρR ∂Ω ∂z , and so in order to parameterize the eddy viscosity as Even though ∂Ω ∂z changes with time, we used for this calculation its initial equilibrium value, as that is the one that would be used in a viscous prescription.Resulting pressureweighted α z averages along constant-θ lines are shown in Fig. 12.In every case, our obtained α z values are between 2 and 4 orders of magnitude larger than α R .For f dg = 10 −4 , α z decreases from values up to 10 −2 to almost-zero at the midplane, where the instability is suppressed.On the other hand, for f dg = 10 −3 at our highest resolution, we obtain α z values that consistently stay between 1 − 2 × 10 −2 , being approximately constant and equal to 2 × 10 −2 close to the midplane, and thus the vertical mixing can be effectively modeled with a constant α z of that value.A similar constancy of α z is seen in Dullemond et al. (2022) and likely also in Stoll et al. (2017), as suggested by their obtained linearity of S ϕz /ρ with z for z/H ≲ 1.

Meridional circulation
We now turn to the average mass fluxes in our simulations.In all runs with f dg = 10 −3 , we obtain an average mass inflow onto the star close to the midplane and an outflow in upper layers (Fig. 13), also seen in Stoll & Kley (2016); Stoll et al. (2017); Manger & Klahr (2018); Pfeil & Klahr (2021).Above |z| ≈ H, the gas moves away from the star with average radial velocities between 10 −5 and 10 −3 au yr −1 (Fig. 13).This so-called meridional circulation results from the vertical transport of angular momentum away from the midplane as determined by Equation ( 7), which causes gas flows at the midplane to fall to lower orbits around the star and surface gas layers to move to upper orbits due to their angular momentum excess.The same phenomenon occurs in isotropic α-viscosity models (e.g., Urpin 1984   the difference that in those cases the circulation pattern is inverted, and the radial velocity as a function of z is shaped as a concave downward parabola.Outward flows produced by meridional circulation can lead to large-scale outward transport of dust particles (Takeuchi & Lin 2002;Stoll & Kley 2016), which has been proposed (e.g., Ciesla 2007) to contribute to the puzzling observed presence of high-temperature minerals in outer regions of protoplanetary disks (Olofsson et al. 2009) and our own solar system (Brownlee et al. (2006), Bryson & Brennecka (2021), and references therein) that likely formed in hotter, inner regions.
A meridional circulation can also be seen in the dustdepleted case, albeit with much smaller mass fluxes, as shown in Fig. 14.In that case, the gas flux is negative close to the irradiation surface and positive for higher altitudes, while at the midplane it oscillates around zero.This occurs because a net upward vertical angular momentum trans- port is only produced at the VSI-active disk surface layers, as shown, e.g., in Figs. 8 and 9. Therefore, outward flows can still occur in upper layers even if the VSI is suppressed at the disk middle layer.Disk regions with an unstable surface layer and a VSI-inactive middle layer can be expected for sufficient depletion of small dust (Section 6 and Paper II).

Thermal evolution
Unlike in isothermal and β-cooling simulations, the average temperature in our simulations is not forced to remain close to its initial value, which allows us to quantify the temperature perturbations produced by the VSI.The ratio ∆T /T 0 = (T − T 0 )/T 0 comparing the temperature distributions T and T 0 , computed respectively at t = 300 and 4 orbits (the latter being after the initial system relaxation but before VSI modes grow), is shown in Fig. 15 for the maximum-resolution runs.
In all of our simulations, the temperature perturbations produced by the VSI are rather negligible.For f dg = 10 −3 , some heating occurs at the disk middle layer due to the radiative dissipation of the mechanical energy of the VSI modes.Thermal equilibrium is reached when radiative cooling balances the VSI heating rate.In our highest-resolution run, this occurs at approximately 250 orbits, i.e., about 100 orbits after the saturation of the body modes, as shown, for instance, by the midplane temperature perturbation profiles in Fig. 16.The temperature perturbations are maximum at the maximum-shear layers in between uniform-j z regions (see Fig. 6), with maximum ∆T /T 0 ∼ 5 × 10 −3 in such locations at the maximum resolution.This pattern is in part produced by the maximum dissipation occurring in such regions, as the heating rate is proportional to α R (Equation (13) later in this section), which reaches a maximum in between uniform-j z bands.On the other hand, the observed temperature perturbations are also affected by the transport of gas along the vertically oriented entropy gradient, given that also vertical velocities are maximum in those regions (Fig. 3).The energy dissipation in those locations may also be connected to the turbulent heating produced by the KHI at the regions of maximum shear and even by the baroclinic deceleration of KH eddies, as discussed in Paper II.
Above the irradiation surface, vertically elongated temperature perturbations cannot longer be seen due to the shorter thermal relaxation timescale in that region (Section 6).Instead, we obtain an overall temperature decrease that becomes larger with radius.This is caused by the slight vertical expansion of the disk due to the heating of the middle layer, which leads to a density increase in the surface layers and thus to additional shadowing of stellar irradiation.Also visible in those regions are faint "shadows" of the reduced density regions produced by small-scale amplified eddies (Paper II).
For f dg = 10 −4 , temperature perturbations are confined to a region near the irradiation surface, where both positive and negative temperature perturbations are produced by the vertical transport of entropy across the temperature transition (Fig. 1).Unlike in the previous case, the average heating produced by the VSI only becomes evident when radially averaging the values of ∆T /T 0 , as shown in Fig. 17.In that case, the VSI produces above the irradiation surface a maximum temperature increase about 10 times smaller than the maximum increase for f dg = 10 −3 in the middle layer.We also see a slight decrease in the midplane temperature due to the disk's reconfiguration over time that is not compensated by VSI heating, which however is reduced with increasing resolution.The overall smaller average ∆T /T 0 in this case results both from the weaker VSI activity in the entire domain (e.g., Figs. 4 and 8) and the lower optical depth in the region where heating takes place (Figs. 1 and 18).Conversely, the relatively larger heating of the middle layer for f dg = 10 −3 results from the fact that the disk is both VSI-active in that region and vertically optically thick.
By Reynolds-averaging the HD equations, the average local heating power density can be computed as where ν R and ν z are the effective eddy viscosities computed in the previous section.We have verified that the term proportional to ν z in this expression can be safely neglected, as it can be shown taking Ω ≈ Ω K that the ratio between the z-and R-contributions Resulting r-averages of the heating rate Q + /(ρϵ) in units of Ω are shown in Fig. 18.For f dg = 10 −3 , this quantity grows with height and is minimum at the midplane, where the disk's baroclinicity is zero and consequently the VSI driving is minimum (Fig. 8).This explains the observed temperature variation, which is locally minimum at the midplane, reaches a maximum a few scale heights above the midplane, and is finally overtaken by the balance of irradiation and cooling at the optically thin surface layers.For f dg = 10 −4 , the heating rate is only nonzero at the VSI-active surface layers, which is not enough to heat the midplane.Due to the vertical dependence of α R , our obtained heating distributions resemble those in the simulations of MRI-heated disks by (Flock et al. 2013), which are also locally minimum at the midplane and peak a few scale heights above it.This observation contradicts the usual assumption resulting from α-viscosity models that turbulent heating is maximum at the midplane, which is commonly used to compute temperature distributions via 1D radiative transfer (e.g., Nakamoto & Nakagawa 1994;Pfeil & Klahr 2019).The vertical temperature distributions resulting from this method peak at the midplane (e.g., Mori et al. 2019), as opposed to our obtained distributions (Fig. 17).

Stability analysis
In this section we analyze the applicability of local and global stability criteria in our disk models.Via a heuris-tic analysis, Urpin (2003) proposed a local Richardson-like criterion in which the VSI is subject to the instability condition where t rel is the local thermal relaxation timescale, while N z is the vertical Brunt-Väisälä frequency (see Fig. 2).A linear analysis in Klahr (2023) for arbitrary cooling times similar to that in Goldreich & Schubert (1967) shows that the linear growth rate decreases as ∼ 1/t rel for t rel < t crit , where the critical cooling time t crit is defined as the right-hand side of Equation ( 14).On the other hand, the vertically global instability criterion by Lin & Youdin (2015) reads where H and Γ are respectively the gas pressure scale height and heat capacity ratio.This expression, which corresponds to Equation ( 14) evaluated at a height z = ΓH/2, is derived in Lin & Youdin (2015) under the assumptions that (I) the disk is vertically isothermal with a temperature distribution depending on the cylindrical radius as T ∝ R q , and (II) that β = t rel /Ω K is vertically uniform, where Ω K is the Keplerian frequency.Most stability studies found in the literature (e.g., Flock et al. 2017b;Dullemond et al. 2022) rely on this criterion, some of them relying on its local application to produce stability maps (e.g., Malygin et al. 2017;Pfeil & Klahr 2019).However, the cited validity conditions cease to be satisfied above the irradiation surface, where both the temperature and the relaxation timescale abruptly vary (Figs. 1 and 19) and the vertical shear does no longer solely depend on the radial temperature gradient (Section 2.4).We therefore analyze the disk stability in terms of the local critical time t crit defined as the righthand side of Equation ( 14) and later discuss in Section 7.5 the local applicability of Equation ( 15).
In the simulations presented in this work, the thermal relaxation time t rel is entirely determined by radiative transfer, and equals 1 the local radiative cooling time, t cool .We estimate this timescale in Appendix B via a local lineal analysis assuming small temperature perturbations with a given wavelength λ on a uniform background, obtaining with where |k| = 2π/λ, λ P = (ρ d κ d P ) −1 and λ R = (ρ d κ d R ) −1 are respectively the Planck-and Rosseland-averaged mean free 1 In reality t rel = Γt cool if one includes the effect of PdV work in the relaxation process (see, e.g., Goldreich & Schubert 1967;Lesur et al. 2023;Klahr 2023), but this widely used approximation is accurate enough for the purpose of a comparison with tcrit, especially considering that the suppression of the VSI for t rel > tcrit as a function of t rel is not abrupt (Lin & Youdin 2015;Manger et al. 2021).
Article number, page 14 of 26 Melon Fuksman, Flock & Klahr: VSI in irradiated PPDs I. Angular momentum transport and turbulent heating .This expression evidences a dependence of the cooling mechanism on λ: temperature perturbations occurring over much smaller lengthscales than λ R and λ P (λ P λ R k 2 ≫ 1) get damped at the scale-independent rate t −1 cool,thin corresponding to thermal emission in an optically thin medium.On the other hand, optically thick perturbations (λ P λ R k 2 ≪ 1) cool down through radiative diffusion at a rate t −1 cool,thick , which increases proportionally to k 2 .The adiabatic limit t cool → ∞ is thus recovered for λ R k → 0. For constant opacities, Equation ( 16) tends to the corresponding expression in Lin & Youdin (2015).In Appendix B it is shown that this expression is exact for any treatment of gray radiative transfer as long as a local first-order perturbation analysis is applicable, but different cooling time prescriptions used in the literature (e.g., Stoll & Kley 2016;Malygin et al. 2017;Dullemond et al. 2022) only differ slightly and should yield equivalent results.In the same appendix, we show that the resulting cooling times are not affected by the RSLA for large enough ĉ, in particular for our employed value of ĉ = 10 −4 .
To estimate the characteristic wavelength of the VSI flows, we computed the distance between consecutive sign changes of v θ at fixed θ values and extrapolate the resulting k distributions to the entire domain.This computation is detailed in Appendix D, where it is shown that kH = O(10) and λ P λ R k 2 ≫ 1 everywhere, i.e., induced temperature perturbations of the same wavelength as the velocity flows are optically thin.This is true at all times, and since the optically thin cooling timescale is wavelength-independent, the estimated t cool values are unchanged through time.This validates the employment of β-cooling prescriptions in even higher-resolution axisymmetric simulations or 3D simulations using the same disk models, with the difference that such an approach cannot reproduce the (rather negligible) disk's thermal evolution (Section 5).Further differences can be expected if large-scale optically thick structures are eventually formed (e.g., vortices in 3D).
Resulting slices of t cool /t crit and t cool Ω at r = 5.5 au are shown in Fig. 19 for each f dg .In all cases, the cooling times are almost constant below the irradiation surface, which would validate in our model the application of a vertically global criterion with uniform cooling time (Equation ( 15)) in that region.Above that surface, t cool decreases up to 2 orders of magnitude despite the surface layers having much smaller optical depths.The reason for this is that the ratio ρϵ/a R T 4 , which stems from the balance between the gas internal energy and the dust emissivity (Appendix B), decreases much faster with height than the increase in λ P (Equation ( 17)).For f dg = 10 −3 , the resulting t cool is everywhere below t crit , which is consistent with the obtained instability in the whole domain.
For f dg = 10 −4 , the local instability criterion is met at the surface layers, which explains our obtained localized VSI activity in those regions.On the other hand, the increase of t cool below the irradiation surface forms two layers above the midplane where t cool > t crit , which explains our obtained VSI suppression in the middle layer.A similar phenomenon is seen in the β-cooling HD simulations by Fukuhara et al. (2023), who obtained that a stable midplane layer is formed in cases where a local application of Equation (15) predicts stability up to heights larger than ∼ H.In our case, Equation 14 predicts stability up to ∼ 3H, and thus our observed VSI suppression is consistent with that work.
Close to the midplane, the instability condition is again met.The reason for this is that both the vertical shear and the squared vertical buoyancy frequency vanish with height, the first as R∂ z Ω ∝ z and the second as N 2 z ∝ z 2 , and so the local critical cooling time close to the midplane diverges as 1/|z|.According to our local analysis, we should then expect a localized increased VSI growth in the unstable region close to the midplane.However, in Fig. 3 we instead see slow (v z ≪ c s ) vertical modes with similar amplitude in the midplane region and in the inactive zones above and below.One possible explanation for this is the fact that the local linear VSI growth rate equals the local vertical shear (e.g., Urpin & Brandenburg 1998;Klahr 2023), which is very small close to the midplane.However, under this criterion, and provided the numerical dissipation in our code is low enough, we should still see some growth occurring after several hundreds of orbits, as it can be seen in the local growth rate distribution shown in Fig. 1.Another explanation can be given if the maximum amplitude of the saturated VSI perturbations is limited by the local vertical shear.A hint in this direction is given by the vertically global model by Lin & Youdin (2015), which shows that the VSI growth rates are reduced if the maximum height of the domain, and therefore also the z-averaged vertical shear, are decreased.This hypothesis is also supported by the approximate relation ⟨α r ⟩ ρ ∝ q 2 in the global vertically isothermal model by Manger et al. (2021), which shows that the radial Reynolds stress in the saturated stage increases for larger average vertical shear.Lastly, resolution might also play a role in limiting VSI growth close to the midplane.The linearly unstable directions, strictly confined between the constant-R and constant-j z surfaces (e.g., Knobloch & Spruit 1982;Klahr et al. 2023), are for small |z| only slightly slanted with respect to the vertical direction, and thus large resolutions might be needed to reproduce linear growth in such directions.This can be investigated in high-resolution box simulations of the midplane region.
However useful, this picture has the limitation that VSI modes are nonlocal.Thus, another possibility is that the amplitude of the flows at the middle layer does not only depend on the in situ driving of the VSI but also on the excitation of inertial modes in unstable layers extending to otherwise inactive regions, as it likely occurs in the T states in Fukuhara et al. (2023).This, together with the fact that VSI modes can in general occupy regions of largely different cooling times, makes it difficult to apply local stability criteria in more complex t cool distributions to predict the saturated behavior of the VSI, for instance if the unstable region in this example is enlarged.Still, our analysis shows that evaluating which specific regions satisfy Equation 14gives a good insight into their expected stability (see also Pfeil & Klahr 2021).

Angular momentum distribution
The formation of bands of approximately uniform j z is likely enhanced by the imposed axisymmetry of our models.In 3D simulations, we would expect constant-j z zones to be unstable to non-axisymmetric modes (e.g., Papaloizou & Pringle 1984), in which case it remains to be seen if the VSI can still enforce the formation of reduced shear bands.Surprisingly, this is likely the case in a number of published works: in the 3D simulations by Richard et al. (2016), the vertical component of the vorticity, ω z , proportional to ∂ R j z for axisymmetric flows, peaks in between VSI modes, indicating abrupt jumps in j z .In between those jumps, ω z /Ω is reduced with respect to its Keplerian value of 1/2, although this reduction is below the value ∆ω z /Ω = −1/2 corresponding to constant j z .More evidence showing abrupt jumps in j z is given by the simulations by Manger & Klahr (2018), in which the azimuthally averaged critical function L, approximately proportional to 1/∂ R j z , has peaks of similar width as the VSI modes.On top of this, the 3D Rad-HD simulations of irradiated disks by Flock et al. (2020) exhibit a formation of approximately uniform-j z bands which is most predominant in the maximum-shear upper layers, as shown in Appendix A (Fig. A.1).It is then possible that non-axisymmetric modes, while preventing the formation of uniform-j z zones, may still permit the creation of regions with reduced ∇j z (see Paper II), although further thorough investigation is required to test this hypothesis.

Reynolds stress comparison
Our obtained vertically averaged α r values for the nominal dust-to-gas ratio f dg = 10 −3 , where the instability is not suppressed close to the midplane, are slightly lower than the values α r ≈ 10 −4 reported in Flock et al. (2020) with H/R = 0.09 − 0.14 and Manger et al. (2021) and Barraza-Alfaro et al. (2021) with H/R = 0.1.On the other hand, Pfeil & Klahr (2021) and Manger et al. (2020) obtained α r ∼ 10 −5 respectively in 2D and 3D simulations with H/R ∼ 0.05.Considering the dependence of α r with H/R studied in Manger et al. (2020Manger et al. ( , 2021)), who propose α r ∝ q 2 (H/R) 2.6 , our α r values obtained with H/R ∼ 0.05−0.06are consistent with those obtained in 3D models.The consistency with vertically isothermal simulations is likely a result of the pressure weighting of α r , which highlights its value at the vertically isothermal middle layer.Additionally, the consistency with 3D models points to the fact that the saturated flows produced by the VSI maintain some degree of axisymmetry, as shown, for example, in the high-resolution simulations in Flock et al. (2020); Barraza-Alfaro et al. (2021) (see, however, Richard et al. 2016;Manger & Klahr 2018), and thus it is likely that the mixing of angular momentum in 3D is predominantly axisymmetric.We defer such verification to follow-up works on this topic.Still, differences in the spatial distribution of averaged velocities and stresses in 2D and 3D simulations are seen in Stoll & Kley (2014); Stoll et al. (2017); Pfeil & Klahr (2021).

Meridional circulation and accretion
Using the obtained mean gas flows, we can compute the total mass flux delivered to the star by integrating ⟨ρv⟩ t onto any surface crossing our domain vertically.When doing so integrating the r-flux along constant-r surfaces, we obtain that the positive and negative contributions cancel each other out almost exactly.More precisely, the θaveraged radial mass flux takes small positive values between 2 and 3 orders of magnitude smaller than the maximum values shown in Fig. 14.This occurs in contradiction with Equation ( 11), which for our simulation parameters predicts that the positivity of ⟨α R ⟩ p and its approximate constancy (Fig. 10) are enough to produce an average mass inflow.The reason for this is that our imposed boundary conditions do not allow for either outflow or inflow of mass, and thus the inward-directed gas flow is forced at the inner radial boundary to turn away from the midplane and merge with the outward-directed meridional flow, while the inverse process occurs at the outer radial boundary.This increases the outward flux, creating an equilibrium configuration in which outgoing and ingoing flows balance each other while no mass escapes the domain.The extent to which this affects the resulting gas structure in our simulations and all mentioned works is unknown.This drawback could be solved by including the inner rim in the domain and allowing gas to escape the radial boundary toward the star, ideally including magnetic effects in well-ionized regions (Flock et al. 2017a).
The ideal MHD simulations of MRI turbulence by Fromang et al. ( 2011) and Flock et al. (2011) show similar downward-concave parabola patterns as in Figure 13, with the difference that, unlike both in Flock et al. (2011) and the present paper, in Fromang et al. (2011) the radial flux remains everywhere positive and there is no average inflow close to the midplane.Given the expected quenching of the MRI by nonideal effects in several disk regions (e.g., Lesur et al. 2023), global nonideal MHD simulations are required to better characterize the radial transport of solids in protoplanetary disks.

Turbulent heating
In our nominal model, heating due to turbulent viscosity is only marginal in comparison with stellar heating, and the average temperature increase produced by the VSI is at most on the order of 0.1 K.This is comparable to the prediction by (Pfeil & Klahr 2019), obtained via a 1D vertical diffusion model, for our obtained ⟨α R ⟩ p values of O(10 −5 ).If we assume a constant α R of that magnitude, we can expect from that work that the VSI can only produce a significant temperature increase well inside of 1 au.However, is it unclear whether the VSI can still produce such α R values in such dense regions, as the denser the disk is, the more the maximum VSI-unstable wavelengths predicted by linear theory get shifted to lower values, as determined by Equation ( 16).In our disk model, as detailed in Paper II, the VSI should be suppressed at the midplane for our maximum average wavenumbers (kH ≈ 70, see Fig. C.3) inside of ∼ 2 au, where only smaller-wavelength modes can grow.Given our obtained trend of decreasing VSI strength with decreasing dominant wavelength (see Figs. 8 and 18), we can expect in that region a suppression of the VSI at the midplane due to the reduction of the maximum VSIunstable wavelengths.It therefore appears unlikely that the VSI can in general lead to any significant heating, since that would simultaneously require (I) large enough vertical optical depths to retain the produced heat and (II) modes with small radial wavelengths producing larger levels of turbulence than predicted by simulations.This hypothesis can be further studied in high-resolution Rad-HD simulations of the inner disk inside of ∼ 1 au, keeping in mind that close enough to the star the MRI is expected to become predominant and suppress the VSI (Lesur et al. 2023).
The model implemented in Pfeil & Klahr (2019) also predicts that larger α R values of ∼ 10 −4 − 10 −3 can lead to an increase of up to a few tens of K outside of 1 au.Such high α R values can only be expected for steeper temperature profiles than in this work (T ∝ R −1 , Pfeil & Klahr 2021;Manger et al. 2021) or for large aspect ratios at a few au of up to 0.1 (Manger et al. 2020), both of which are hard to justify in passive irradiated disks.In this context, it is also unclear whether steeper temperature gradients, predicted by models of viscous-dominated heating in optically thick accretion disks (q = − 3 4 , Pringle 1981), can be maintained solely by VSI-induced turbulence.Moreover, the cited α-disk model in Pfeil & Klahr (2019) assumes maximum heating at the midplane, whereas our maximum heating regions occur in upper regions of lower optical depth, and thus the temperature increase in that work can be regarded as an upper limit in the case of VSIinduced heating.Further explorations are however required to test whether any heating can be produced by the VSI under a plausible choice of parameters.
The fact that the temperature distributions in our nominal disk are at all times minimum close to the r-boundaries (Fig. 16) evidences that the obtained temperature increase is not only limited by vertical diffusion but also by radiative cooling through the radial boundaries.This is a numerical effect produced by the employed Dirichlet boundary conditions, which leads to outflow of radiation as soon as the temperature inside of the domain exceeds that at the boundary.Larger domains could be used in the future to explore whether this effect limits the maximum temperature perturbations.
Determining the equilibrium temperatures produced by the VSI is also made difficult by the lack of convergence of the temperature distributions with resolution due to the reduction of the VSI strength, as shown in Fig. 17.However, the average temperatures obtained with our two highest resolutions are rather similar, with differences of ∼ 10%.Whether such resolutions are enough to accurately predict the temperature increase produced by the VSI can be explored in the future via even higher-resolution simulations.However, given our obtained trend of decreasing heating with increasing resolution, this effect should not alter the conclusion that the heating produced by the VSI is most likely negligible.

Local and vertically global stability criteria
In Fig. 19 it can be seen that the validity conditions of the global stability criterion by Lin & Youdin (2015), namely that the disk is vertically isothermal and that the cooling time is vertically uniform, are approximately met in the entire extent of the disk middle layer.Furthermore, the locations z = ±ΓH/2 at which t crit equals its value in the global stability criterion (Equation ( 15)) are well inside the middle layer (Fig. 19).This suggests that we can expect the global stability criterion to be applicable in the disk's middle layer, with a stability condition depending on the critical time t glob crit = t crit (ΓH/2).This is in fact verified in our disk models, since t cool is smaller and larger than t glob crit for f dg = 10 −3 and 10 −4 , respectively (see Fig. 19).Further evidence for this is given by the vertically isothermal β-cooling simulations by Fukuhara et al. (2023), which only show a stable middle layer if the height ∆L s /2 of the stability region predicted by the local application of Equation ( 15) is at least ∼ H.Given the functional form of t cool (z) in that work, ∆L s /2 must be larger than z = ΓH/2 to verify t cool ≫ t crit at z = ±ΓH/2.Thus, the midplane stability for ∆L s /2 ≳ H in that work is consistent with the vertically global criterion by Lin & Youdin (2015), even though, unlike these authors, Fukuhara et al. (2023) did not consider a vertically uniform cooling time.These results indicate that the cited vertically global criterion should suffice to predict the stability of regions below the irradiation surface.
The vertical constancy of t cool in the middle layer is not guaranteed if the radial optical depth of the VSI modes is large enough that the cooling timescale is determined by radiative diffusion (Equation ( 16)), as shown in Fig. 2 of Flock et al. (2017b).It remains to be explored if global stability criteria are still applicable to such cases.
Lastly, we have not considered in this work the effect on the thermal relaxation timescale of the dust-gas thermal equilibration time (see, e.g., Barranco et al. 2018).This becomes important at the disk upper layers, where the dust and gas densities are low enough that t rel is increased due to their infrequent collisions, with a stabilizing effect in those regions (Pfeil & Klahr 2021).On the other hand, t rel may be reduced by molecular emissivity depending on the chemical composition, favoring the instability.These effects can be studied in future simulations decoupling the gas and dust temperatures, as done in Muley et al. (2023).The impact of these phenomena on the disk stability is explored in terms of both local and global criteria in Paper II.

Conclusions
In this work we studied the linear and nonlinear evolution of the VSI in protoplanetary disks, focusing on the transport of angular momentum, the produced temperature perturbations, and the applicability of local stability con-ditions.We employed for this axisymmetric two-moment radiation-hydrodynamical simulations with realistic temperature stratifications and cooling times, both of which result from the balance of stellar irradiation and radiative cooling, which are treated in our code using realistic tabulated dust opacities.We assumed typical parameters for disks around T Tauri stars considering two different values of the dust-to-gas mass ratio of small grains, namely a nominal case with f dg = 10 −3 and a dust-depleted case with f dg = 10 −4 , making our studied disk region between 4 and 7 au vertically optically thick and thin, respectively.Our main findings can be summarized as follows: 1.The temperature gradient in our disk models leads to a localized increase of the vertical shear in the temperature transition region between the midplane and the atmosphere.The radial and vertical components of the temperature gradient tend to increase and decrease the overall vertical shear, respectively.At its peak, the resulting vertical shear reaches up to 2 times its value in equivalent vertically isothermal models, which is translated into similarly larger local VSI growth rates than predicted by linear theory for vertically isothermal disks.The stabilizing buoyancy frequency is also maximum in that region, while above it decreases below its vertically isothermal value.
2. In our nominal case, we observe two distinct growth stages corresponding to the growth and saturation of finger and body modes occupying the entire domain.
Instead, in the dust-depleted case, the VSI is suppressed at the disk middle layer, and body modes never grow.This behavior results in a single growth phase of VSI modes localized at the disk surface layers.
3. In every case, the VSI reduces the vertical shear inside of bands of approximately uniform specific angular momentum.Conversely, shear is maximal at the interfaces between these bands.Even though this process is likely enhanced by the enforced axisymmetry, bands of reduced shear can also be seen in 3D simulations, suggesting that their formation may in some cases predominate over their destruction by nonaxisymmetric modes.
4. The time-and volume-averaged absolute value of the vertical shear increases with resolution as a result of the growth of small-scale velocity fluctuations, exceeding the mean shear value by almost one order of magnitude in our highest-resolution runs.This is not translated into an increase of the VSI strength, which is instead enforced by the disk's baroclinicity and limited by the VSI saturation mechanism.
5. The Reynolds stress is in every case maximum at the maximum shear regions above the irradiation surfaces and in between the approximately uniform-j z bands, with z components between 1 and 2 orders of magnitude larger than the R components.The pressure-weighted α r values of ∼ 10 −5 obtained in the nominal case coincide with reported values for 3D vertically isothermal simulations.This suggests that, even in 3D, the radial mixing of angular momentum may be predominantly axisymmetric.
6.The anisotropy of the Reynolds stress tensor causes differences between its radial components in cylindrical and spherical coordinates, with α R < α r by factors of ∼ 4 and ∼ 2 in the nominal and dust-depleted cases, respectively.Determining whether this is a general trend for different disk parameters requires further testing.While α R should be used to compute the vertically integrated mass flux, α r is usually reported in the literature.
7. As in previous works (e.g., Stoll & Kley 2016), the anisotropic angular momentum transport produced by the VSI produces a meridional circulation with inflow at the midplane and outflow in upper layers.This pattern only occurs at the disk's surface layer in the dustdepleted case.This circulation does not lead to accretion in our simulations due to our imposed no-outflow and no-inflow boundary conditions, which likely increases the resulting outflow velocities with respect to their values if mass could escape the domain.
8. The turbulent heating produced by the VSI in our models causes negligible temperature increases of 0.1% and 0.01% in the nominal and dust-depleted models, respectively.In the nominal case, the energy dissipated by the VSI modes heats the disk middle layer with a heating rate distribution peaking a few (∼ 3) scale heights above the midplane, contrarily to the heating distributions in α-viscosity models, which peak at the midplane.We conclude that the VSI is in general unlikely to produce any significant temperature increase, as that would either require it to produce larger heating rates than expected in dense disk regions inside of ∼ 1 au, where largely optically thick radial wavelengths cease to be VSI-unstable, or unrealistically large radial temperature gradients or aspect ratios leading to larger average α R values (10 −4 − 10 −3 ) than can be obtained in models of passive disks.Further explorations of the parameter space via high-resolution Rad-HD simulations of inner disk regions should either deny or confirm this hypothesis.
9. The Richardson-like local criterion by Urpin (2003) explains the instability regions observed in our simulations, showing that, despite VSI modes are vertically global, a local analysis can provide good insight into the stability of specific disk regions.On the other hand, the global criterion in Lin & Youdin (2015) derived for vertically isothermal disks with uniform cooling time can still be applied to the vertically isothermal disk middle layer, which justifies its general application to assess the stability close to the midplane as long as VSI modes are radially optically thin.
10.We verified via analytical computations and numerical testing that the RSLA is safely applicable to this problem, since for high enough ĉ < c the cooling times are not affected by the artificial reduction of the speed of light.This is important due to the fact that the VSI is sensitive in both its linear and nonlinear regimes to the precise value of t cool if this value is similar or larger than t crit .This also shows that IMEX schemes for Rad-HD, with generally good scalability properties in comparison with globally implicit methods, can be applied to this problem without introducing unphysical phenom-ena, which is rather useful given the high resolutions required to resolve the VSI modes.
Our stability analysis does not contemplate effects due to nonzero thermal coupling timescales between dust grains and gas particles and gas thermal emissivity.These can affect the stability of the low-density disk surface layers, in particular in outer disk regions currently observable with ALMA (e.g., Andrews et al. 2018).These effects are included in Paper II to study the stability of such regions.In run dg3c4_256 (top left), it can be seen that the transition between finger and body modes is still not complete.
In principle we can obtain all three roots of P by applying a cubic formula, but instead we can get accurate analytic expressions for these roots by considering the limits of small and large |s|.
We first consider the limit which also guarantees |s| ≪ ĉ/λ P , since our employed opacities verify κ d P > κ d R for all temperatures (Melon Fuksman & Klahr 2022) (this is also true in general for opacity laws of the form κ d ν ∝ ν n as long as n > 0.81).In that case, the characteristic equation P (s) = 0 loses its dependence on ĉ and becomes This expression can be simplified by noting that the term iξ ′ λ R k is much smaller than 1 for large λ (λ R k ≲ 1) and much smaller than the O(λ R λ P k 2 ) term for small λ (λ R k ≫ 1), and thus it can be neglected.On the other hand, the factor ξ − ξ ′ f only affects the value of s c if λ R λ P k 2 is not much larger than 1, which can only occur if f ≪ 1, in which case ξ ′ f ≈ f 2 ≪ 1 and ξ ≈ 1/3.Therefore we can safely replace ξ −ξ ′ f by 1/3, which leads to the simpler expression Using the k distributions in our simulations (Section 6), we estimated the relative error of this approximation to be of at most 0.01%.If this analysis is repeated in 3D, the same expressions are obtained with the difference that λ P λ R k 2 is replaced by a term of the same order of magnitude that tends to λ P λ R k 2 /3 wherever λ P λ R k 2 ≲ 1.Thus, the particular form of the O(λ P λ R k 2 ) terms for anisotropic transport does not affect the resulting cooling rate, which always tends to s thin for λ P λ R k 2 ≫ 1 and to Equation (B.9) otherwise.Consequently, Equation (B.9) still holds when considering multidimensional radiation transport.Analogously, we would have obtained the same result using any other closure than M1 or even solving the full frequency-integrated radiative transfer equation, since P ij r always tends to E r δ ij /3 for λ P λ R k 2 ≲ 1.Therefore, Equation (B.9) is valid in general for gray radiative transfer.As we show in Section B.4, this eigenvalue alone determines the radiative cooling rate, which can be computed as which is the constant-opacity limit of Equation ( 16).The wavelength-dependent behavior of t cool can be understood by noting that the value of k regulates how much radiation energy is transported away through advection before it can reheat the dust-gas mixture via absorption.This can be seen by comparing the relative size of the different terms in Eq. (B.2).For large k, E r quickly evolves into a state such that δE r ≈ ηδp/(1 + O(λ P λ R k 2 )) ≪ ηδp, which can be seen by zeroing the time derivatives of the radiation fields.This means that δE r can be dropped from the cooling source term of δp, which is proportional to (ηδp−δE r ).In this way, radiation transport limits the maximum value of δE r enough to cause temperature perturbations to be damped solely due to optically thin cooling, which is scale-independent.This is not verified for smaller k values if λ P λ R k 2 /3 is not much larger than 1, in which case (ηδp−δE r ) ≈ ηδp λ R λ P k 2 3+λ R λ P k 2 .For small enough k, diffusion overcomes optically thin cooling and becomes the sole relaxation mechanism.
The above derivation relies on the constancy of the opacity coefficients.If we consider instead that opacities depend on the local temperature, as we do in this work, the resulting eigenvalue corresponding to thermal relaxation is 3 ∇E r , it is straightforward to prove that the relative variation of s c due to that term is at most 4b R λ/L, where L is the variation lengthscale of E r .Given the initial assumption that λ ≪ L, that term can be neglected in all regimes and the resulting cooling time is given by Equation (16).

B.2. Consistency check and constraints on ĉ
We now test the consistency of the small |s c | condition assumed in Eq. (B.6).This limit can be better understood by considering the components of δE r and δF x r that evolve proportionally to e sct+ikx , which we denote as δE r and δF (B.12) meaning that if s c satisfies Eq. (B.6), the time derivatives of the radiation fields can be dropped.This means that, in that case, δE r and δF x r evolve quasi-statically, as their readjustment time is much shorter than |s c | −1 , which explains why ĉ disappears from the equations.In particular, if f ≪ 1, this leads to the usual expression for the radiative flux in the diffusion limit, F r = − 1 3ρχ R ∇E r , whose perturbation in this analysis becomes δF For ĉ = c, this equation is satisfied close to the midplane, where η(Γ − 1)λ R /λ P ∼ 10 −5 − 10 −4 , but not at the disk atmosphere.However, our expression for s c is still valid in that region, as we explain next.Far enough from the midplane, the timescale |s thin | −1 becomes shorter than the light-crossing time of a distance λ R , typically of a few hundred au in those regions.In such locations, λ R λ P k 2 → ∞ (otherwise λ ∼ λ R/P , which would render the uniform background assumption inapplicable), in which case s = s thin provided a large-k condition is satisfied independently from Eq. (B.13).To see this, we rearrange the characteristic equation as that is, as long as either the cooling timescale is much shorter than the light-crossing time of a mean free path or the gas-dust mixture is largely optically thin along a perturbation wavelength.For the wavelengths in this work (see Section 6), this relation holds in the entire domain for all considered ĉ values in Section B.5, except for ĉ/c = 10 −5 in the case f dg = 10 −3 , in which the left-hand side of Eq. (B.16) is of order 1 close to the midplane.We further study the effect of reducing ĉ on the cooling timescale in Section B.5.

B.3. Remaining eigenvalues
Having verified the validity of Eq. (B.9), it only remains to show that δp ∼ δp 0 e sct , for which we need to compute the order of magnitude of the contribution to δp of the remaining eigenvalues.These can be obtained with good accuracy by assuming |s| ≫ |s c |, which transforms the characteristic equation into  Since we only want to obtain the order of magnitude of the two solutions to this equation, s ± , it is enough to compute their values in the optically thick and thin limits.If λ R k ≪ 1, which also guarantees λ P k ≪ 1, the two solutions of these equations are while for λ P k ≫ 1, and therefore also λ R k ≫ 1, we have B.19) where the argument of the square root is always positive if f < 1.Therefore, |s ± | is at least of order ∼ max(λ −1 R ĉ, kĉ).We can then check the consistency of the assumption |s ± | ≫ |s c |, for which it is enough to verify that |s ± | ≫ |s thin |, since |s thin | ≥ |s c |.For λ R k ≪ 1, this is equivalent to Eq. (B.13), which as already mentioned is verified in that regime for large enough ĉ.The same is true for λ R k ≫ 1, in which case the resulting condition is ĉ/c ≫ η(Γ − 1)/(λ P k), which in that limit holds everywhere for ĉ = c.

B.4. Exponential solution
With the obtained eigenvalues, we can compute an expression for δp(t) and compare the relative size of the terms proportional to {e sct , e s+t , e s−t }.The leading orders of δp = (e At (δp 0 , 0, 0) ⊺ ) 1 can be evaluated in a rather straightforward way by applying the Cayley-Hamilton theorem, which allows us to compute e At as a linear combination of {I, A, A 2 } with time-dependent coefficients, where I is the 3 × 3 identity matrix.This results in expressions of the form if λ P k ≫ 1 , where O i (x) denotes a term of order x.For ĉ = c, the leading term in these equations is proportional to e sct , whereas the other terms are much smaller, both close to the midplane, where η ≪ 1, as well as away from the midplane, where λ P k is always much larger than 1.This proves that t −1 cool = |s c |. On the other hand, the relative size of these terms changes for ĉ < c, which together with the changes in the eigenvalues can affect the overall relaxation rate for sufficiently small ĉ.We study this effect in the next section.

B.5. Dependence on the reduced speed of light
In the previous sections we showed that t cool can be different from its physical value when ĉ/c < 1, either due to changes in the value of s c or in the relative size of the terms in Equations (B.20) and (B.21).We now evaluate the overall change in the cooling time when ĉ/c < 1 by numerically computing δp(t) in our domain, defining t cool (ĉ) as the e-folding time of δp, i.e., δp(t cool (ĉ))/δp(0) = e −1 .This definition coincides with the real cooling time when c = ĉ.We achieve this by numerically computing e At for different times in the interval [0, 2 t cool (c)] using the SciPy (Virtanen et al. 2020) expm matrix exponentiation routine, which follows the method in Al-Mohy & Higham (2010).The A matrix is in turn computed at each location using the k distributions obtained in Appendix D for N θ = 512.
We repeated this computation for different ĉ/c values in the interval

Fig. 1 .
Fig.1.Initial temperature and vertical shear distributions.The latter also correspond to the VSI growth rates predicted by linear theory for instant cooling, shown here in units of the inverse orbital time at 5.5 au (T −1 5.5 ).Dashed, white lines indicate the τP (Ts) = 1 surface for stellar photons computed using Planck-averaged opacities, solid white lines show the τP = 1 surfaces for vertically traveling infrared photons, and green dash-dotted lines correspond to the z = ±H surfaces.

Fig. 2 .
Fig. 2. Vertical slices at r = 5.5 au in the initial hydrostatic state.Top panel: Temperature profiles.Middle panels: R∂zΩ 2 profiles compared to the corresponding values for a vertically isothermal disk (Ω = Ωiso), and different terms contributing to the vertical shear in Eq. (5) for f dg = 10 −3 (second from the top) and 10 −4 (third from the top) in units of the local squared Keplerian frequency Ω 2 K .Bottom panel: Squared vertical Brunt-Väisälä frequency for both f dg values in units of Ω 2 K (solid lines) compared to their values for vertically isothermal disks (dotted lines).

Fig. 3 .
Fig. 3. Vertical flows produced by the VSI.Top row: vertical velocity distributions normalized by the isothermal sound speed cs = p/ρ in our highest-resolution simulations after 300 orbits (also shown after 50 orbits for f dg = 10 −3 ).Bottom row: vertical mass flux in the same snapshots, increased by a factor of 100 for f dg = 10 −4 .
Fig. 4. R.m.s.velocities computed between 150 and 300 orbits for different resolutions and dust content, averaged at constant r (left) and constant θ (right).

Fig. 6 .Fig. 7 .
Fig.6.Vertical shear and angular momentum.Left panels: comparison between the vertical shear at t = 0 and t = 300 orbits in the highest-resolution runs.Right panels: corresponding specific angular momentum distributions after 300 orbits.

Fig. 8 .
Fig. 8. Reynolds stress cylindrical components normalized by the average pressure, averaged over r at constant θ for varying resolution and dust content.

Fig. 9 .
Fig. 9. Norm of the ϕ-row of the Reynolds stress tensor normalized by the average pressure at the highest resolution.

Fig. 10 .
Fig. 10.Pressure-weighted z-averages of αR computed for different resolutions and f dg as in Equation (10).

Fig. 12 .
Fig. 12. Averages at constant θ of αz computed as in Equation (12) for varying resolution and f dg .

Fig. 13 .
Fig. 13.Absolute value of the average radial velocity at r = 5.5 au in our highest-resolution simulations.Blue and orange lines correspond respectively to outward and inward velocities, respectively.The pressure scale height is indicated for reference (gray dotted line).

Fig. 14 .
Fig. 14.Average mass flux in the radial direction at r = 5.5 au for varying resolution and dust content.

Fig. 15 .
Fig. 15.Relative temperature variation after 300 orbits at the maximum employed resolution.
Fig. 16.Midplane relative temperature perturbation at different times in run dg3c4_2048.

Fig. 17 .Fig. 18 .
Fig. 17.Radially averaged relative temperature variation for all employed resolutions and f dg values.Displayed values have been averaged in time between 285 and 300 orbits to get rid of fluctuations.The values for f dg = 10 −4 have been multiplied by 10 to facilitate the comparison.

Fig. 19 .
Fig. 19.Radiative cooling times normalized by the critical cooling time tcrit (top) and the orbital frequency (bottom) at r = 5.5 au for both f dg values.Also shown are the corresponding values for t = tcrit and the locations z = ±ΓH/2 (gray dotted lines) at which the local tcrit coincides with its value in the global stability criterion (Equation (15)).

Fig. A. 2 .
Fig. A.2. Velocity distributions normalized by the isothermal sound speed after 300 orbits for varying resolution and dust content.In run dg3c4_256 (top left), it can be seen that the transition between finger and body modes is still not complete.

Fig. B. 1 .
Fig. B.1.Relative variation of the cooling rate with respect to its value when c = ĉ as a function of ĉ/c for f dg = 10 −3 (top row) and 10 −4 (bottom row).

Fig
Fig. C.3.Wavenumber histograms (normalized using the local scale height) as a function of radius for all resolutions and f dg values, measured at z/r = 0.1 and 0.28 (top and bottom panels, respectively).

Fig. D. 1 .
Fig. D.1.Estimated distributions of the parameter λP λRk 2 for all resolutions and f dg values, indicating the optical depth of the VSI modes (λP λRk 2 ≫ 1 indicates optically thin cooling).