Issue 
A&A
Volume 682, February 2024



Article Number  A139  
Number of page(s)  26  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/202346554  
Published online  16 February 2024 
Vertical shear instability in twomoment radiationhydrodynamical simulations of irradiated protoplanetary disks
I. Angular momentum transport and turbulent heating
Max Planck Institute for Astronomy, Königstuhl 17, 69117 Heidelberg, Germany
email: fuksman@mpia.de
Received:
31
March
2023
Accepted:
9
November
2023
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 stratification, 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 VSIstability criteria.
Methods. We modeled the gasdust mixture in our disk models by means of highresolution axisymmetric radiationhydrodynamical simulations including stellar irradiation with frequencydependent opacities. This is the first study of this instability to rely on twomoment 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 implicitexplicit methods, thus avoiding the employment of expensive global matrix solvers. This is done at the cost of artificially reducing the speed of light, which, as we verified for this work, does not introduce unphysical phenomena.
Results. Given sufficient depletion of small grains (with a dusttogas 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 specific 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 VSIinduced 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 justifiable 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 dustdepleted 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 significant temperature increase since that would either require it to efficiently operate in largely optically thick disk regions or to produce larger levels of turbulence than predicted by models of passive irradiated disks.
Key words: hydrodynamics / instabilities / radiative transfer / methods: numerical / protoplanetary disks
© The Authors 2024
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model.
Open Access funding provided by Max Planck Society.
1. Introduction
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, largescale 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 largescale 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 currently 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 scaleindependent, 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 largerscale 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 (RadHD) codes. Moreover, unlike for isothermal and βcooling prescriptions, such codes can be used to quantify the turbulent heating produced in VSIunstable 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 RadHD 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 lowdensity 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. 2024, Paper II henceforth).
We employed resolutions of up to ∼200 cells per scale height, which are unprecedented in RadHD 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 RadHD method applied to the study of disk instabilities not to rely on globally implicit solvers.
This article is organized as follows. In Sect. 2 we outline our employed numerical scheme and describe our disk models. In Sect. 3 we characterize the velocity and vertical shear distributions in the different growth and saturation phases of the VSI, while in Sect. 4 we analyze the resulting redistribution of angular momentum in relation with the small and largescale flows in our simulations. In Sect. 5 we explore the thermal evolution of our simulated disks, and in Sect. 6 we analyze our obtained stability regions in terms of local criteria. Finally, in Sect. 7 we discuss our results, and in Sect. 8 we summarize our conclusions. Additional figures and calculations are included in the appendices.
2. Numerical method and disk models
2.1. Radiation hydrodynamics
We modeled the gasdust mixture in a circumstellar disk and its heating and cooling processes by means of the RadHD module (Melon Fuksman et al. 2021) implemented in the opensource PLUTO code (version 4.4, Mignone et al. 2007). In the chosen configuration for this work, the following system of equations is solved:
$$\begin{array}{cc}& \frac{\partial \rho}{\partial t}+\mathrm{\nabla}\xb7\left(\rho \mathbf{v}\right)=0\hfill \\ \hfill & \frac{\partial (\rho \mathbf{v})}{\partial t}+\mathrm{\nabla}\xb7\left(\rho \mathbf{v}\mathbf{v}\right)+\mathrm{\nabla}p=\mathbf{G}\rho \mathrm{\nabla}\mathrm{\Phi}\hfill \\ \hfill & \frac{\partial (E+\rho \mathrm{\Phi})}{\partial t}+\mathrm{\nabla}\xb7[(E+p+\rho \mathrm{\Phi})\mathbf{v}]=c\phantom{\rule{0.166667em}{0ex}}{G}^{0}\mathrm{\nabla}\xb7{\mathbf{F}}_{\mathrm{Irr}}\hfill \end{array}$$
$$\begin{array}{cc}& \frac{1}{\widehat{c}}\frac{\partial {E}_{r}}{\partial t}+\mathrm{\nabla}\xb7{\mathbf{F}}_{r}={G}^{0}\hfill \\ \hfill & \frac{1}{\widehat{c}}\frac{\partial {\mathbf{F}}_{r}}{\partial t}+\mathrm{\nabla}\xb7{\mathbb{P}}_{r}=\mathbf{G},\hfill \end{array}$$(1)
where ρ, p, and v are the gas density, pressure and velocity, while E_{r}, F_{r}, and ℙ_{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=\rho \u03f5+\frac{1}{2}\rho {\mathbf{v}}^{2}$ imposing an ideal equation of state to the internal energy density, namely $\rho \u03f5=\frac{p}{\mathrm{\Gamma}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 ${P}_{r}^{ij}={D}^{ij}{E}_{r}$, where the Eddington tensor D^{ij} is defined as ${D}^{\mathit{ij}}=\frac{1\xi}{2}\phantom{\rule{0.166667em}{0ex}}{\delta}^{\mathit{ij}}+\frac{3\xi 1}{2}{n}^{i}{n}^{j}$, with $\xi =\frac{3+4{f}^{2}}{5+2\sqrt{43{f}^{2}}}$, where n = F_{r}/F_{r}, f = F_{r}/E_{r}, 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
$$\begin{array}{c}\hfill {\mathbf{F}}_{\mathrm{Irr}}(r,\theta )=\pi {\left(\frac{{R}_{\mathrm{s}}}{r}\right)}^{2}{\displaystyle {\int}_{{\nu}_{\mathrm{min}}}^{{\nu}_{\mathrm{max}}}}\mathrm{d}\nu \phantom{\rule{0.166667em}{0ex}}{B}_{\nu}({T}_{\mathrm{s}})\phantom{\rule{0.166667em}{0ex}}{e}^{{\tau}_{\nu}(r,\theta )}\phantom{\rule{0.166667em}{0ex}}\widehat{\mathbf{r}},\end{array}$$(2)
where T_{s} is the star temperature, R_{s} is the star radius, B_{ν}(T)=(2hν^{3}/c^{2})/(e^{hν/kBT} − 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=\frac{\mu u}{{k}_{\mathrm{B}}}\frac{p}{\rho}$, 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 (Sect. 2.2), we estimate this value following the prescription in Melon Fuksman & Klahr (2022). In the RadHD simulations, we retrieve it from the initial conditions.
Optical depths are computed considering only opacities due to dust absorption, which at optical and nearinfrared 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 dusttogas 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 frequencyintegrated radiationmatter interaction terms, which in the gas comoving frame have expressions
$$\begin{array}{cc}& {\stackrel{\sim}{G}}^{0}={\kappa}_{\mathrm{P}}^{d}\phantom{\rule{0.166667em}{0ex}}{\rho}_{\mathrm{d}}({\stackrel{\sim}{E}}_{r}{a}_{\mathrm{R}}{T}^{4})\hfill \\ \hfill & \stackrel{\sim}{\mathbf{G}}={\chi}_{\mathrm{R}}^{d}\phantom{\rule{0.166667em}{0ex}}{\rho}_{\mathrm{d}}\phantom{\rule{0.166667em}{0ex}}{\stackrel{\sim}{\mathbf{F}}}_{r},\hfill \end{array}$$(3)
where tildes indicate comoving quantities, while a_{r} = 4σ_{SB}/c is the radiation constant, σ_{SB} is the StefanBoltzmann constant, and ${\kappa}_{\mathrm{P}}^{d}$ and ${\chi}_{\mathrm{R}}^{d}$ are, respectively, the Planck and Rosselandaveraged absorption and extinction coefficients. The laboratory frame coefficients (G^{0}, G) are obtained via a Lorentz transformation of $({\stackrel{\sim}{G}}^{0},\stackrel{\sim}{\mathbf{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 ${\chi}_{\mathrm{R}}^{d}={\kappa}_{\mathrm{R}}^{d}$, which is the Rosselandaveraged absorption opacity, although the derivations in Appendix B consider a general ${\chi}_{\mathrm{R}}^{d}$. Graphs of the frequencydependent and frequencyaveraged 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 $\widehat{c}<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 CourantFriedrichsLewy condition using the real speed of light, thus reducing the overall computational cost. The obtained solutions of the RadHD equations are in general unchanged for large enough $\widehat{c}$ values, as we later verify in our setup.
2.2. 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 selfshadowing (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.04 M_{⊙}. If a vertically integrated dusttogas 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 dusttogas 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 dustdepleted case with f_{dg} = 10^{−4}. If we assume a total dusttogas 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 smallsize population 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).
2.3. Hydrostatic distributions
In our RadHD 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 (Sect. 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.
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}_{5.5}^{1}$). Dashed, white lines indicate the τ_{P}(T_{s})=1 surface for stellar photons computed using Planckaveraged opacities, solid white lines show the τ_{P} = 1 surfaces for vertically traveling infrared photons, and green dashdotted lines correspond to the z = ±H surfaces. 
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}{\left(\frac{R}{4\phantom{\rule{0.166667em}{0ex}}\mathrm{au}}\right)}^{\frac{1q}{2}}$, 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, $\frac{1q}{2}$, is in every case positive and close to 1/4.
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 ${\mathrm{\Omega}}_{K}^{2}$. Bottom panel: squared vertical BruntVäisälä frequency for both f_{dg} values in units of ${\mathrm{\Omega}}_{K}^{2}$ (solid lines) compared to their values for vertically isothermal disks (dotted lines). 
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 Planckaveraged 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 traveling reprocessed photons, defined as the vertical τ_{P} = 1 surfaces computed with ${\kappa}_{\mathrm{P}}^{d}(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 Eq. (1),
$$\begin{array}{cc}\hfill T& ={\left(\frac{\mathrm{\nabla}\xb7{\mathbf{F}}_{\mathrm{Irr}}}{c{\rho}_{d}{\kappa}_{\mathrm{P}}^{d}(T){a}_{\mathrm{R}}}\right)}^{1/4}\hfill \\ \hfill & \approx {\left(\frac{{\kappa}_{\mathrm{P}}^{d}({T}_{\mathrm{s}})}{4{\kappa}_{\mathrm{P}}^{d}(T)}\right)}^{1/4}{\left(\frac{{R}_{\mathrm{s}}}{r}\right)}^{1/2}{e}^{{\tau}_{\mathrm{P}}({T}_{\mathrm{s}})/4}\phantom{\rule{0.166667em}{0ex}}{T}_{\mathrm{s}},\hfill \end{array}$$(4)
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 powerlaw 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).
2.4. 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 Eq. (1)) in cylindrical coordinates (R, ϕ, z), which leads to the following relations:
$$\begin{array}{cc}& {\kappa}_{z}^{2}\equiv R{\partial}_{z}{\mathrm{\Omega}}^{2}\hfill \\ \hfill & \phantom{\rule{1em}{0ex}}=\frac{{(\mathrm{\nabla}p\times \mathrm{\nabla}\rho )}_{\varphi}}{{\rho}^{2}}\hfill \\ \hfill & \phantom{\rule{1em}{0ex}}={(\mathrm{\nabla}\overline{T}\times \mathrm{\nabla}log\rho )}_{\varphi}\hfill \\ \hfill & \phantom{\rule{1em}{0ex}}={\partial}_{z}\overline{T}{\partial}_{\mathrm{R}}log\rho {\partial}_{\mathrm{R}}\overline{T}{\partial}_{z}log\rho ,\hfill \end{array}$$(5)
where ${\kappa}_{z}^{2}=\frac{1}{{R}^{3}}{\partial}_{z}{j}_{z}^{2}$, with an analogous expression for the epicyclic frequency ${\kappa}_{\mathrm{R}}^{2}=\frac{1}{{R}^{3}}{\partial}_{\mathrm{R}}{j}_{z}^{2}$, quantifies the vertical variation of the specific angular momentum j_{z} = R^{2}Ω, while Ω = v_{ϕ}/R is the orbital angular frequency and $\overline{T}=p/\rho $. This equation, which is a form of the socalled 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, ${\kappa}_{z}^{2}$ 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 ${\partial}_{z}\overline{T}$ becomes important, and the obtained vertical shear results from the balance of both terms on the righthand side of Eq. (5), which have opposite signs. In particular, ${\partial}_{z}\overline{T}{\partial}_{\mathrm{R}}log\rho $ is positive above the midplane and negative below, and the opposite is true for ${\partial}_{\mathrm{R}}\overline{T}{\partial}_{z}log\rho $. 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 ${\partial}_{z}\overline{T}{\partial}_{\mathrm{R}}log\rho $ and ${\partial}_{\mathrm{R}}\overline{T}{\partial}_{z}log\rho $ are about 5 times larger than the maximum ${\kappa}_{z}^{2}$, 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 fastestgrowing mode is ${\mathrm{\Gamma}}_{\mathrm{VSI}}=\frac{{\kappa}_{z}^{2}}{2\mathrm{\Omega}}={\partial}_{z}{v}_{\varphi}$ (Urpin & 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 BruntVäisälä frequency corresponding to vertical buoyant oscillations, ${N}_{z}^{2}=\frac{1}{\rho \mathrm{\Gamma}}\frac{\partial p}{\partial z}\frac{\partial log(p/{\rho}^{\mathrm{\Gamma}})}{\partial z}$ (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 Sect. 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.
2.5. RadHD simulations’ setup
Starting from the described initial conditions, we ran RadHD 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 $\widehat{c}/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.
Varying parameters of the presented simulations and derived quantities.
In every run considered in the next sections we used $\widehat{c}/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 $\widehat{c}/c$, respectively. In particular, it is shown there that the radiative cooling time is unaffected for high enough $\widehat{c}/c$, which is crucial for this work given the dependence on the VSI on that timescale.
The RadHD equations are solved using thirdorder RungeKutta 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 thirdorder 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.
3. VSI growth and saturation
3.1. 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 highestresolution simulations are shown in Fig. 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 Sect. 6, this reduced activity results from the stabilizing effect of buoyancy for long enough cooling times.
Fig. 3. Vertical flows produced by the VSI. Top row: vertical velocity distributions normalized by the isothermal sound speed ${c}_{\mathrm{s}}=\sqrt{p/\rho}$ in our highestresolution 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}. 
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 (Sect. 6).
The gas rms 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 ℳ ≈ 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} (ℳ ∼ 0.24) and ∼100 m s^{−1} (ℳ ∼ 0.12) for f_{dg} = 10^{−3} and 10^{−4}, respectively. At the middle layer, rms 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 rms 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. A.2.
Fig. 4. Rms velocities computed between 150 and 300 orbits for different resolutions and dust content, averaged at constant r (left) and constant θ (right). 
The different growth phases of the instability can be determined by evaluating the time evolution of the volumeaveraged meridional kinetic energy, $\u27e8{\u03f5}_{k}\u27e9=\u27e8\frac{\rho}{2}({v}_{r}^{2}+{v}_{\theta}^{2})\u27e9$, shown in Fig. 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), and FloresRivera 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.
Fig. 5. Time evolution of the volumeaveraged kinetic energy (top) and maximum gas velocity (bottom) in all runs with $\widehat{c}/c={10}^{4}$. 
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 KelvinHelmholtz (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 KelvinHelmholtz 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 KHunstable 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 Fig. 5, which explains why the maximum velocities increase with resolution despite the overall decrease of the kinetic energy.
3.2. 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 Eq. (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 ${\partial}_{z}{v}_{\varphi}=\frac{1}{R}{\partial}_{z}{j}_{z}$. As explored in Paper II, these bands begin forming due to the angular momentum redistribution produced by the motion of VSIdestabilized gas parcels across constantj_{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 uniformj_{z} regions are formed (see Paper II and Klahr et al. 2023).
Fig. 6. Vertical shear and angular momentum. Left panels: comparison between the vertical shear at t = 0 and t = 300 orbits in the highestresolution runs. Right panels: corresponding specific angular momentum distributions after 300 orbits. 
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 quasiKeplerian initial state (not shown here), since v_{ϕ} = j_{z}/R ∼ R^{−1/2} for a quasiKeplerian flow and ∼R^{−1} for constant j_{z}. In each case, the shear in all directions is always maximum in between the uniformj_{z} zones.
Overall, the increase of the shear in between uniformj_{z} regions, together with the extra shear produced by KH eddies, result in an increase of the time and volumeaveraged 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 VSIinactive 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 smallscale 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 VSI. As shown by the decrease with resolution of both the saturated ⟨ϵ_{k}⟩ values and the components of the Reynolds stress tensor (Sect. 4.1), the increasing amplitude of the shear fluctuations with resolution 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 (see Paper II).
Fig. 7. Same as Fig. 4, but showing instead constantr and constantθ averages of the averaged absolute value of the vertical shear between 150 and 300 orbits. The corresponding values at t = 0 are shown for comparison. 
4. Angular momentum transport
4.1. Reynolds stresses
We now turn to the transport of angular momentum in the disk, which depends on the timeaveraged 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
$$\begin{array}{c}\hfill {\mathbf{S}}_{\varphi}={\u27e8\delta {v}_{\varphi}\delta (\rho \mathbf{v})\u27e9}_{t}={\u27e8\rho {v}_{\varphi}\mathbf{v}\u27e9}_{t}{\u27e8{v}_{\varphi}\u27e9}_{t}{\u27e8\rho \mathbf{v}\u27e9}_{t},\end{array}$$(6)
where ⟨ ⋅ ⟩_{t} indicates time average and δ denotes fluctuations with respect to the mean. The components of this tensor can be turned into stresstopressure 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 VSIinduced angular momentum transport cannot be modeled by means of an isotropic viscosity coefficient $\nu =\alpha {c}_{\mathrm{s}}^{2}/{\mathrm{\Omega}}_{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 Rcomponent 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:
$$\begin{array}{c}\hfill {\partial}_{t}(\u27e8\rho \u27e9{j}_{z})+\u27e8\rho \mathbf{v}\u27e9\xb7\mathrm{\nabla}{j}_{z}+\mathrm{\nabla}\xb7(R{\mathbf{S}}_{\varphi})=0,\end{array}$$(7)
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 VSIinduced mixing, respectively. In quasi steadystate 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, Eq. (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
$$\begin{array}{c}\hfill {F}_{\mathrm{\Sigma}}\equiv {\displaystyle {\int}_{\infty}^{\infty}}\u27e8\rho {v}_{\mathrm{R}}\u27e9\mathrm{d}z\approx \frac{2}{{R}^{1/2}}\frac{\partial}{\partial R}\left(\frac{{R}^{1/2}}{\mathrm{\Omega}}{\displaystyle {\int}_{\infty}^{\infty}}{S}_{\varphi \mathrm{R}}\mathrm{d}z\right)\end{array}$$(9)
(see, e.g., Jacquet 2013), where Ω is evaluated at the midplane. In turn, the column density evolves according to
$$\begin{array}{c}\hfill {\partial}_{t}\mathrm{\Sigma}+\frac{1}{R}\frac{\partial (R{F}_{\mathrm{\Sigma}})}{\partial R}=0,\end{array}$$(9)
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 pressureweighted 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 uniformj_{z} bands described in Sect. 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.
Fig. 8. Reynolds stress cylindrical components normalized by the average pressure, averaged over r at constant θ for varying resolution and dust content. 
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.
Fig. 9. Norm of the ϕrow of the Reynolds stress tensor normalized by the average pressure at the highest resolution. 
To define an effective Reynolds α value that can be translated into effective eddy viscosity models of the form ${S}_{\varphi \mathrm{R}}={\nu}_{\mathrm{R}}\rho R\frac{\partial \mathrm{\Omega}}{\partial R}$ (i.e., the ϕRcomponent of the strain rate tensor of a fluid with viscosity ν_{R}), while reproducing the usual parameterization ${S}_{\varphi \mathrm{R}}=\alpha p$ for thin, vertically isothermal disks, we define
$$\begin{array}{c}\hfill {\alpha}_{\mathrm{R}}=\frac{3}{2}\frac{{S}_{\varphi \mathrm{R}}}{p}{\mathrm{\Omega}}_{K}{\left(R\phantom{\rule{0.166667em}{0ex}}\frac{\partial \mathrm{\Omega}}{\partial R}\right)}^{1},\end{array}$$(10)
or equivalently, ${\nu}_{\mathrm{R}}=\frac{2}{3}{\alpha}_{\mathrm{R}}{c}_{\mathrm{s}}^{2}{\mathrm{\Omega}}_{K}^{1}$ (see, e.g., Armitage 2022), where Ω_{K}(R) is the Keplerian angular velocity (the extra $\frac{2}{3}$ factor accounts for the fact that $R\frac{\partial \mathrm{\Omega}}{\partial R}\approx \frac{3}{2}{\mathrm{\Omega}}_{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., LyndenBell & Pringle (1974), we can then replace ${S}_{\varphi \mathrm{R}}={\alpha}_{\mathrm{R}}{\u27e8p\u27e9}_{t}$ in Eq. (8), which leads to
$$\begin{array}{c}\hfill {F}_{\mathrm{\Sigma}}\approx \frac{2}{{R}^{1/2}}\frac{\partial}{\partial R}\left({R}^{1/2}{\u27e8{\alpha}_{\mathrm{R}}\u27e9}_{\mathrm{p}}\frac{{c}_{\mathrm{s}}^{2}}{\mathrm{\Omega}}\mathrm{\Sigma}\right),\end{array}$$(11)
where ⟨α_{R}⟩_{p} is the pressureweighted zaverage 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}_{\mathrm{s}}^{2}=p/\rho $, which is a valid approximation given the temperature distribution in the disk middle layer.
Values of ⟨α_{R}⟩_{p} computed via interpolation along constantR 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}.
Fig. 10. Pressureweighted zaverages of α_{R} computed for different resolutions and f_{dg} as in Eq. (10). 
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}_{\varphi z}\gg {S}_{\varphi \mathrm{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 Eq. (11)).
Fig. 11. α_{r} values for varying resolution and dust content computed as S_{ϕr}/⟨p⟩_{t}, averaged at constant r (left) and constant θ (right). 
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 ϕzcomponent of the strain rate tensor in cylindrical coordinates is of the form ${S}_{\varphi z}={\nu}_{z}\rho R\frac{\partial \mathrm{\Omega}}{\partial z}$, and so in order to parameterize the eddy viscosity as ${\nu}_{z}={\alpha}_{z}{c}_{\mathrm{s}}^{2}{\mathrm{\Omega}}_{K}^{1}$ we define
$$\begin{array}{c}\hfill {\alpha}_{z}=\frac{{S}_{\varphi z}}{p}{\mathrm{\Omega}}_{K}{\left(R\phantom{\rule{0.166667em}{0ex}}\frac{\partial \mathrm{\Omega}}{\partial z}\right)}^{1}.\end{array}$$(12)
Even though $\frac{\partial \mathrm{\Omega}}{\partial 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 almostzero 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.
4.2. 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, 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 socalled meridional circulation results from the vertical transport of angular momentum away from the midplane as determined by Eq. (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; Jacquet 2013; Kley & Lin 1992; Philippov & Rafikov 2017), with 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 largescale 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 hightemperature 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.
Fig. 13. Absolute value of the average radial velocity at r = 5.5 au in our highestresolution simulations. Blue and orange lines correspond respectively to outward and inward velocities, respectively. The pressure scale height is indicated for reference (gray dotted line). 
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 transport is only produced at the VSIactive 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 VSIinactive middle layer can be expected for sufficient depletion of small dust (Sect. 6 and Paper II).
Fig. 14. Average mass flux in the radial direction at r = 5.5 au for varying resolution and dust content. 
5. 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 maximumresolution runs.
Fig. 15. Relative temperature variation after 300 orbits at the maximum employed resolution. 
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 highestresolution 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 maximumshear layers in between uniformj_{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} (Eq. (13) later in this section), which reaches a maximum in between uniformj_{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.
Fig. 16. Midplane relative temperature perturbation at different times in run dg3c4_2048. 
Above the irradiation surface, vertically elongated temperature perturbations cannot longer be seen due to the shorter thermal relaxation timescale in that region (Sect. 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 smallscale 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 VSIactive in that region and vertically optically thick.
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. 18. Radial average of the heating rate produced by the VSI in units of Ω for all employed resolutions and f_{dg} values. 
By Reynoldsaveraging the HD equations, the average local heating power density can be computed as
$$\begin{array}{cc}& {Q}^{+}={S}_{\varphi \mathrm{R}}R{\partial}_{\mathrm{R}}\mathrm{\Omega}{S}_{\varphi z}R{\partial}_{z}\mathrm{\Omega}\hfill \\ \hfill & \phantom{\rule{1em}{0ex}}\phantom{\rule{4pt}{0ex}}={\nu}_{\mathrm{R}}\rho \phantom{\rule{0.166667em}{0ex}}{(R{\partial}_{\mathrm{R}}\mathrm{\Omega})}^{2}+{\nu}_{z}\rho \phantom{\rule{0.166667em}{0ex}}{(R{\partial}_{z}\mathrm{\Omega})}^{2},\hfill \end{array}$$(13)
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 Rcontributions to Q^{+} is approximately $\frac{1}{36}\frac{{\alpha}_{z}}{{\alpha}_{\mathrm{R}}}{\left(\frac{z}{R}\right)}^{2}$, which is ≪1 despite α_{z} ≫ α_{R}. Resulting raverages 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 VSIactive 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 MRIheated 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).
6. Stability analysis
In this section we analyze the applicability of local and global stability criteria in our disk models. Via a heuristic analysis, Urpin (2003) proposed a local Richardsonlike criterion in which the VSI is subject to the instability condition
$$\begin{array}{c}\hfill {t}_{\mathrm{rel}}\lesssim \frac{{\partial}_{z}{v}_{\varphi}}{{N}_{z}^{2}},\end{array}$$(14)
where t_{rel} is the local thermal relaxation timescale, while N_{z} is the vertical BruntVä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 righthand side of Eq. (14). On the other hand, the vertically global instability criterion by Lin & Youdin (2015) reads
$$\begin{array}{c}\hfill {t}_{\mathrm{rel}}<\frac{q}{\mathrm{\Gamma}1}\frac{H}{R}\phantom{\rule{0.166667em}{0ex}}{\mathrm{\Omega}}_{K}^{1},\end{array}$$(15)
where H and Γ are respectively the gas pressure scale height and heat capacity ratio. This expression, which corresponds to Eq. (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 (Sect. 2.4). We therefore analyze the disk stability in terms of the local critical time t_{crit} defined as the righthand side of Eq. (14) and later discuss in Sect. 7.5 the local applicability of Eq. (15).
Fig. 19. Radiative cooling times normalized by the critical cooling time t_{crit} (top) and the orbital frequency (bottom) at r = 5.5 au for both f_{dg} values. Also shown are the corresponding values for t = t_{crit} and the locations z = ±ΓH/2 (gray dotted lines) at which the local t_{crit} coincides with its value in the global stability criterion (Eq. (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
$$\begin{array}{c}\hfill {t}_{\mathrm{cool}}={t}_{\mathrm{cool},\mathrm{thin}}+{t}_{\mathrm{cool},\mathrm{thick}}\end{array}$$(16)
with
$$\begin{array}{cc}& {t}_{\mathrm{cool},\mathrm{thin}}=\frac{{\lambda}_{\mathrm{P}}\zeta}{c}\hfill \\ \hfill & {t}_{\mathrm{cool},\mathrm{thick}}=\frac{3\zeta}{c\phantom{\rule{0.166667em}{0ex}}{\lambda}_{\mathrm{R}}{k}^{2}}\hfill \\ & \zeta =\frac{\rho \u03f5}{4\phantom{\rule{0.166667em}{0ex}}{a}_{\mathrm{R}}{T}^{4}+{b}_{\mathrm{P}}({a}_{\mathrm{R}}{T}^{4}{E}_{r})},\hfill \end{array}$$(17)
where k=2π/λ, ${\lambda}_{\mathrm{P}}=({\rho}_{d}{\kappa}_{\mathrm{P}}^{d}{)}^{1}$ and ${\lambda}_{\mathrm{R}}=({\rho}_{d}{\kappa}_{\mathrm{R}}^{d}{)}^{1}$ are respectively the Planck and Rosselandaveraged mean free paths, and ${b}_{\mathrm{P}}=\frac{dlog{\kappa}_{\mathrm{P}}^{d}}{dlogT}$. 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 scaleindependent rate ${t}_{\mathrm{c}\mathrm{o}\mathrm{o}\mathrm{l},\mathrm{t}\mathrm{h}\mathrm{i}\mathrm{n}}^{1}$ 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}_{\mathrm{c}\mathrm{o}\mathrm{o}\mathrm{l},\mathrm{t}\mathrm{h}\mathrm{i}\mathrm{c}\mathrm{k}}^{1}$, which increases proportionally to k^{2}. The adiabatic limit t_{cool} → ∞ is thus recovered for λ_{R}k → 0. For constant opacities, Eq. (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 firstorder 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 $\widehat{c}$, in particular for our employed value of $\widehat{c}={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 = 𝒪(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 wavelengthindependent, the estimated t_{cool} values are unchanged through time. This validates the employment of βcooling prescriptions in even higherresolution 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 (Sect. 5). Further differences can be expected if largescale 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 (Eq. (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} (Eq. (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 Eq. (15) predicts stability up to heights larger than ∼H. In our case, Eq. (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}_{z}^{2}\propto {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 zaveraged 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 constantR and constantj_{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 highresolution 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 Eq. (14) gives a good insight into their expected stability (see also Pfeil & Klahr 2021).
7. Discussion
7.1. 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 constantj_{z} zones to be unstable to nonaxisymmetric 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 ℒ, approximately proportional to 1/∂_{R}j_{z}, has peaks of similar width as the VSI modes. On top of this, the 3D RadHD simulations of irradiated disks by Flock et al. (2020) exhibit a formation of approximately uniformj_{z} bands which is most predominant in the maximumshear upper layers, as shown in Fig. A.1. It is then possible that nonaxisymmetric modes, while preventing the formation of uniformj_{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.
7.2. Reynolds stress comparison
Our obtained vertically averaged α_{r} values for the nominal dusttogas 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 BarrazaAlfaro 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. (2020, 2021), who propose α_{r} ∝ q^{2}(H/R)^{2.6}, our α_{r} values obtained with H/R ∼ 0.05 − 0.06 are 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 highresolution simulations in Flock et al. (2020), BarrazaAlfaro 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 followup 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, 2017), Pfeil & Klahr (2021).
7.3. 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 rflux along constantr 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 Eq. (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 inwarddirected gas flow is forced at the inner radial boundary to turn away from the midplane and merge with the outwarddirected 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 wellionized regions (Flock et al. 2017a).
The ideal MHD simulations of MRI turbulence by Fromang et al. (2011) and Flock et al. (2011) show similar downwardconcave parabola patterns as in Fig. 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.
7.4. 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 𝒪(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 VSIunstable wavelengths predicted by linear theory get shifted to lower values, as determined by Eq. (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 smallerwavelength 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 highresolution RadHD 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 viscousdominated heating in optically thick accretion disks ($q=\frac{3}{4}$, Pringle 1981), can be maintained solely by VSIinduced 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 rboundaries (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 higherresolution 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.
7.5. 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 (Eq. (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}_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}^{\mathrm{g}\mathrm{l}\mathrm{o}\mathrm{b}}={t}_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}(\mathrm{\Gamma}H/2)$. This is in fact verified in our disk models, since t_{cool} is smaller and larger than ${t}_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}^{\mathrm{g}\mathrm{l}\mathrm{o}\mathrm{b}}$ 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 Eq. (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 (Eq. (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 dustgas 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.
8. 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 conditions. We employed for this axisymmetric twomoment radiationhydrodynamical 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 dusttogas mass ratio of small grains, namely a nominal case with f_{dg} = 10^{−3} and a dustdepleted 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 dustdepleted 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 mechanism may in some cases predominate over their destruction by nonaxisymmetric modes.

4.
The time and volumeaveraged absolute value of the vertical shear increases with resolution as a result of the growth of smallscale velocity fluctuations, exceeding the mean shear value by almost one order of magnitude in our highestresolution 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 uniformj_{z} bands, with z components between 1 and 2 orders of magnitude larger than the R components. The pressureweighted α_{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 dustdepleted 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 nooutflow and noinflow 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 dustdepleted 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 VSIunstable, 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 highresolution RadHD simulations of inner disk regions should either deny or confirm this hypothesis.

9.
The Richardsonlike 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 $\widehat{c}<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 RadHD, with generally good scalability properties in comparison with globally implicit methods, can be applied to this problem without introducing unphysical phenomena, 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 lowdensity 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 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 t_{crit}, especially considering that the suppression of the VSI for t_{rel} > t_{crit} as a function of t_{rel} is not abrupt (Lin & Youdin 2015; Manger et al. 2021).
Acknowledgments
The research of J.D.M.F. and H.K. is supported by the German Science Foundation (DFG) under the priority program SPP 1992: “Exoplanet Diversity” under contract KL 1469/161/2. We thank our collaboration partners on this project in Kiel, Sebastian Wolf and Anton Krieger, under contract WO 857/171/2, for providing us with the employed tabulated opacity coefficients. M.F. acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No. 757957). All numerical simulations were run on the ISAAC and VERA clusters of the MPIA and the COBRA cluster of the Max Planck Society, all of these hosted at the MaxPlanck Computing and Data Facility in Garching (Germany). We thank the anonymous referee, whose comments helped us improve the quality of this work.
References
 AlMohy, A. H., & Higham, N. J. 2010, SIAM J. Matrix Anal. Appl., 31, 970 [CrossRef] [Google Scholar]
 Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41 [NASA ADS] [CrossRef] [Google Scholar]
 Armitage, P. J. 2022, ArXiv eprints [arXiv:2201.07262] [Google Scholar]
 Bai, X.N. 2015, ApJ, 798, 84 [Google Scholar]
 Bai, X.N., & Stone, J. M. 2013, ApJ, 769, 76 [Google Scholar]
 Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [Google Scholar]
 Balbus, S. A., & Papaloizou, J. C. B. 1999, ApJ, 521, 650 [Google Scholar]
 Barker, A. J., & Latter, H. N. 2015, MNRAS, 450, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Barranco, J. A., Pei, S., & Marcus, P. S. 2018, ApJ, 869, 127 [NASA ADS] [CrossRef] [Google Scholar]
 BarrazaAlfaro, M., Flock, M., Marino, S., & Pérez, S. 2021, A&A, 653, A113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Béthune, W., Lesur, G., & Ferreira, J. 2017, A&A, 600, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brownlee, D., Tsou, P., Aléon, J., et al. 2006, Science, 314, 1711 [CrossRef] [PubMed] [Google Scholar]
 Bryson, J. F. J., & Brennecka, G. A. 2021, ApJ, 912, 163 [NASA ADS] [CrossRef] [Google Scholar]
 Burn, R., Emsenhuber, A., Weder, J., et al. 2022, A&A, 666, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chandrasekhar, S. 1961, Hydrodynamic and Hydromagnetic Stability (Oxford: Clarendon Press) [Google Scholar]
 Ciesla, F. J. 2007, Science, 318, 613 [NASA ADS] [CrossRef] [Google Scholar]
 Cui, C., & Bai, X.N. 2022, MNRAS, 516, 4660 [NASA ADS] [CrossRef] [Google Scholar]
 Decampli, W. M., Cameron, A. G. W., Bodenheimer, P., & Black, D. C. 1978, ApJ, 223, 854 [NASA ADS] [CrossRef] [Google Scholar]
 Dullemond, C. P., Ziampras, A., Ostertag, D., & Dominik, C. 2022, A&A, 668, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Flock, M., Dzyurkevich, N., Klahr, H., Turner, N. J., & Henning, T. 2011, ApJ, 735, 122 [NASA ADS] [CrossRef] [Google Scholar]
 Flock, M., Fromang, S., González, M., & Commerçon, B. 2013, A&A, 560, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Flock, M., Nelson, R. P., Turner, N. J., et al. 2017a, ApJ, 850, 131 [Google Scholar]
 Flock, M., Fromang, S., Turner, N. J., & Benisty, M. 2017b, ApJ, 835, 230 [CrossRef] [Google Scholar]
 Flock, M., Turner, N. J., Nelson, R. P., et al. 2020, ApJ, 897, 155 [Google Scholar]
 FloresRivera, L., Flock, M., & Nakatani, R. 2020, A&A, 644, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fromang, S., Lyra, W., & Masset, F. 2011, A&A, 534, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fukuhara, Y., Okuzumi, S., & Ono, T. 2023, PASJ, 75, 233 [NASA ADS] [CrossRef] [Google Scholar]
 Gerbig, K., MurrayClay, R. A., Klahr, H., & Baehr, H. 2020, ApJ, 895, 91 [CrossRef] [Google Scholar]
 Goldreich, P., & Schubert, G. 1967, ApJ, 150, 571 [Google Scholar]
 Gressel, O., Ramsey, J. P., Brinch, C., et al. 2020, ApJ, 896, 126 [Google Scholar]
 Hartmann, L., Calvet, N., Gullbring, E., & D’Alessio, P. 1998, ApJ, 495, 385 [Google Scholar]
 Jacquet, E. 2013, A&A, 551, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 James, H. A., & Kahn, F. D. 1970, A&A, 5, 232 [NASA ADS] [Google Scholar]
 Johansen, A., Oishi, J. S., Mac Low, M.M., et al. 2007, Nature, 448, 1022 [Google Scholar]
 Johansen, A., Blum, J., Tanaka, H., et al. 2014, in Protostars and Planets VI, 547 [Google Scholar]
 Klahr, H. 2023, A&A, submitted [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Klahr, H. H., & Bodenheimer, P. 2003, ApJ, 582, 869 [NASA ADS] [CrossRef] [Google Scholar]
 Klahr, H., & Schreiber, A. 2020, ApJ, 901, 54 [NASA ADS] [CrossRef] [Google Scholar]
 Klahr, H., Baehr, H., & Melon Fuksman, J. D. 2023, ApJ, submitted [arXiv:2305.08165] [Google Scholar]
 Kley, W., & Lin, D. N. C. 1992, ApJ, 397, 600 [CrossRef] [Google Scholar]
 Knobloch, E., & Spruit, H. C. 1982, A&A, 113, 261 [NASA ADS] [Google Scholar]
 Krieger, A., & Wolf, S. 2020, A&A, 635, A148 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Krieger, A., & Wolf, S. 2022, A&A, 662, A99 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Latter, H. N., & Papaloizou, J. 2018, MNRAS, 474, 3110 [NASA ADS] [CrossRef] [Google Scholar]
 Lesur, G., Flock, M., Ercolano, B., et al. 2023, ASP Conf. Ser., 534, 465 [NASA ADS] [Google Scholar]
 Levermore, C. D. 1984, J. Quant. Spectr. Rad. Transf., 31, 149 [NASA ADS] [CrossRef] [Google Scholar]
 Lin, M.K. 2019, MNRAS, 485, 5221 [CrossRef] [Google Scholar]
 Lin, M.K., & Youdin, A. N. 2015, ApJ, 811, 17 [NASA ADS] [CrossRef] [Google Scholar]
 LyndenBell, D., & Pringle, J. E. 1974, MNRAS, 168, 603 [Google Scholar]
 Lyra, W., & Umurhan, O. M. 2019, PASP, 131, 072001 [Google Scholar]
 Malygin, M. G., Klahr, H., Semenov, D., Henning, T., & Dullemond, C. P. 2017, A&A, 605, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Manara, C. F., Rosotti, G., Testi, L., et al. 2016, A&A, 591, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Manara, C. F., Morbidelli, A., & Guillot, T. 2018, A&A, 618, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Manger, N., & Klahr, H. 2018, MNRAS, 480, 2125 [NASA ADS] [CrossRef] [Google Scholar]
 Manger, N., Klahr, H., Kley, W., & Flock, M. 2020, MNRAS, 499, 1841 [NASA ADS] [CrossRef] [Google Scholar]
 Manger, N., Pfeil, T., & Klahr, H. 2021, MNRAS, 508, 5402 [NASA ADS] [CrossRef] [Google Scholar]
 Marcus, P. S., Pei, S., Jiang, C.H., et al. 2015, ApJ, 808, 87 [NASA ADS] [CrossRef] [Google Scholar]
 Melon Fuksman, J. D., & Klahr, H. 2022, ApJ, 936, 16 [NASA ADS] [CrossRef] [Google Scholar]
 Melon Fuksman, J. D., & Mignone, A. 2019, ApJS, 242, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Melon Fuksman, J. D., Klahr, H., Flock, M., & Mignone, A. 2021, ApJ, 906, 78 [NASA ADS] [CrossRef] [Google Scholar]
 Melon Fuksman, J. D., Flock, M., & Klahr, H. 2024, A&A, 682, A140 (Paper II) [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mignone, A. 2014, J. Comput. Phys., 270, 784 [Google Scholar]
 Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228 [Google Scholar]
 Mori, S., Bai, X.N., & Okuzumi, S. 2019, ApJ, 872, 98 [NASA ADS] [CrossRef] [Google Scholar]
 Muley, D., Melon Fuksman, J. D., & Klahr, H. 2023, A&A, 678, A162 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nakamoto, T., & Nakagawa, Y. 1994, ApJ, 421, 640 [Google Scholar]
 Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610 [Google Scholar]
 Olofsson, J., Augereau, J. C., van Dishoeck, E. F., et al. 2009, A&A, 507, 327 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ormel, C. W., & Cuzzi, J. N. 2007, A&A, 466, 413 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Papaloizou, J. C. B., & Pringle, J. E. 1984, MNRAS, 208, 721 [NASA ADS] [CrossRef] [Google Scholar]
 Pfeil, T., & Klahr, H. 2019, ApJ, 871, 150 [NASA ADS] [CrossRef] [Google Scholar]
 Pfeil, T., & Klahr, H. 2021, ApJ, 915, 130 [NASA ADS] [CrossRef] [Google Scholar]
 Pfeil, T., Birnstiel, T., & Klahr, H. 2023, ApJ, 959, 121 [NASA ADS] [CrossRef] [Google Scholar]
 Philippov, A. A., & Rafikov, R. R. 2017, ApJ, 837, 101 [NASA ADS] [CrossRef] [Google Scholar]
 Pringle, J. E. 1981, ARA&A, 19, 137 [NASA ADS] [CrossRef] [Google Scholar]
 Richard, S., Nelson, R. P., & Umurhan, O. M. 2016, MNRAS, 456, 3571 [NASA ADS] [CrossRef] [Google Scholar]
 Rüdiger, G., Arlt, R., & Shalybkov, D. 2002, A&A, 391, 781 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schäfer, U., Johansen, A., & Banerjee, R. 2020, A&A, 635, A190 [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, IAU Symp., 55, 155 [Google Scholar]
 Simon, J. B., Bai, X.N., Flaherty, K. M., & Hughes, A. M. 2018, ApJ, 865, 10 [NASA ADS] [CrossRef] [Google Scholar]
 Skinner, M. A., & Ostriker, E. C. 2013, ApJS, 206, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Stoll, M. H. R., & Kley, W. 2014, A&A, 572, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Stoll, M. H. R., & Kley, W. 2016, A&A, 594, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Stoll, M. H. R., Kley, W., & Picogna, G. 2017, A&A, 599, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Svanberg, E., Cui, C., & Latter, H. N. 2022, MNRAS, 514, 4581 [NASA ADS] [CrossRef] [Google Scholar]
 Takeuchi, T., & Lin, D. N. C. 2002, ApJ, 581, 1344 [NASA ADS] [CrossRef] [Google Scholar]
 Toro, E. F. 2009, Riemann Solvers and Numerical Methods for Fluid Dynamics (Berlin Heidelberg: Springer) [Google Scholar]
 Townsend, A. A. 1958, J. Fluid Mech., 4, 361 [NASA ADS] [CrossRef] [Google Scholar]
 Urpin, V. 2003, A&A, 404, 397 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Urpin, V. A. 1984, Sov. Astron., 28, 50 [NASA ADS] [Google Scholar]
 Urpin, V., & Brandenburg, A. 1998, MNRAS, 294, 399 [Google Scholar]
 Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]
 Voelkel, O., Klahr, H., Mordasini, C., & Emsenhuber, A. 2022, A&A, 666, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Yamaleev, N. K., & Carpenter, M. H. 2009, J. Comput. Phys., 228, 4248 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Additional figures
Figure A.1 shows a constantϕ slice of the specific angular momentum in the 3D RadHD disk simulation by Flock et al. (2020) after 350 inner orbits of evolution. Approximately uniformj_{z} bands are formed predominantly at the disk upper layers, where, as in the present work, the perturbations produced by the VSI are maximal. These bands are approximately axisymmetric, and consequently they are also visible in azimuthal averages of the same quantity. Figure A.2 shows the velocity distributions in our simulations after 300 orbits for all of our employed resolutions.
Fig. A.1. Polar slice showing the specific angular momentum in the 3D disk simulation by Flock et al. (2020) after 350 orbits of evolution. 
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. 
Appendix B: Radiative cooling timescale
We estimate the thermal relaxation timescale due to radiative processes by computing the damping time of small temperature perturbations in a hydrostatic disk. In this analysis we neglect mass advection and assume a constant ρ and $\mathbf{v}={v}_{\varphi}\widehat{\mathit{\varphi}}$ distribution taken from the initial condition of each RadHD simulation. We can therefore write the temperature in terms of the pressure using the ideal gas law, which simplifies the form of the resulting equations and the following analysis. We write the pressure and the radiation fields as (p, E_{r}, F_{r})=(p_{0} + δp, E_{r0} + δE_{r}, F_{r0} + δF_{r}), where (p_{0}, E_{r0}, F_{r0}) are the hydrostatic unperturbed values and (δp, δE_{r}, δF_{r}) are timedependent perturbations. Inserting this expression into Equation (1) and keeping only linear terms in the perturbations, we obtain the following system of equations:
$$\begin{array}{cc}\hfill {\partial}_{t}\delta p& ={\kappa}_{P}^{d}{\rho}_{d}\phantom{\rule{0.166667em}{0ex}}c\phantom{\rule{0.166667em}{0ex}}(\mathrm{\Gamma}1)(\eta \phantom{\rule{0.166667em}{0ex}}\delta p\delta {E}_{r})\hfill \\ \hfill {\partial}_{t}\delta {E}_{r}& ={\kappa}_{P}^{d}{\rho}_{d}\phantom{\rule{0.166667em}{0ex}}\widehat{c}\phantom{\rule{0.166667em}{0ex}}(\eta \phantom{\rule{0.166667em}{0ex}}\delta p\delta {E}_{r})\widehat{c}\mathrm{\nabla}\xb7\delta {\mathbf{F}}_{r}\hfill \\ \hfill {\partial}_{t}\delta {\mathbf{F}}_{r}& =\widehat{c}\phantom{\rule{0.166667em}{0ex}}\mathrm{\nabla}\xb7\delta {\mathbb{P}}_{r}{\rho}_{d}{\chi}_{R}^{d}\widehat{c}\phantom{\rule{0.166667em}{0ex}}\delta {\mathbf{F}}_{r}\phantom{\rule{0.166667em}{0ex}},\hfill \end{array}$$(B.1)
where the quantity η = 4 a_{R}T^{4}/p is evaluated at the unperturbed state, while source terms are written neglecting terms of order v/c and opacity variations with temperature. The last of these equations considers a general total opacity coefficient ${\chi}_{R}^{d}$, while in this work we neglect scattering and assume ${\chi}_{R}^{d}={\kappa}_{R}^{d}$. We then simplify this expression by considering onedimensional radiation transport along a single direction n_{0} = F_{r0}/F_{r0}. In this way the radiation pressure tensor takes the simple form ${P}_{r}^{xx}=\xi {E}_{r}$, where x is the only spatial coordinate. This simplification does not alter the resulting cooling rates, as later justified in Section B.1.
To evaluate the dependence of cooling rates on the perturbation lengthscale λ, we consider perturbations of the form $(\delta \widehat{p}(k,t),\delta \widehat{{E}_{r}}(k,t),\delta \widehat{{F}_{r}^{x}}(k,t))\phantom{\rule{0.166667em}{0ex}}{e}^{\mathit{ikx}}$, where k ∈ ℝ is the wavenumber, which verifies k=2π/λ. Assuming that the background fields are uniform in space, or equivalently, that their variation lengthscales are much larger than λ, Eq. (B.1) becomes the following homogeneous system of ordinary differential equations:
$$\begin{array}{c}\hfill \frac{\partial}{\partial t}\left(\begin{array}{c}\delta p\\ \delta {E}_{r}\\ \delta {F}_{r}^{x}\end{array}\right)=A\left(\begin{array}{c}\delta p\\ \delta {E}_{r}\\ \delta {F}_{r}^{x}\end{array}\right)\phantom{\rule{0.166667em}{0ex}},\end{array}$$(B.2)
where the matrix A is defined as
$$\begin{array}{c}\hfill A=\left(\begin{array}{ccc}\eta c\phantom{\rule{0.166667em}{0ex}}(\mathrm{\Gamma}1)/{\lambda}_{P}& c\phantom{\rule{0.166667em}{0ex}}(\mathrm{\Gamma}1)/{\lambda}_{P}& 0\\ \eta \widehat{c}/{\lambda}_{P}& \widehat{c}/{\lambda}_{P}& ik\widehat{c}\\ 0& i\widehat{c}k(\xi {\xi}^{\prime}f)& \widehat{c}(1/{\lambda}_{R}+i{\xi}^{\prime}k)\end{array}\right),\end{array}$$(B.3)
while ${\lambda}_{P}=({\rho}_{d}{\kappa}_{P}^{d}{)}^{1}$ and ${\lambda}_{R}=({\rho}_{d}{\chi}_{R}^{d}{)}^{1}$ are respectively the Planck and Rosseland mean free paths and ${\xi}^{\prime}=\mathrm{d}\xi /\mathrm{d}f=2f/\sqrt{43{f}^{2}}$, which is always positive for f > 0 and vanishes for f = 0.
Given an initial temperature perturbation, which in this analysis takes the form (δp_{0}, 0, 0), the timedependent solution of Eq. (B.2) can be obtained via matrix exponentiation as e^{At}(δp_{0}, 0, 0)^{⊺}. The matrix A has three different eigenvalues, {s_{c}, s_{+}, s_{−}}, and therefore all perturbations evolve as a linear combination of {e^{sct}, e^{s+t}, e^{s−t}}. We derive expressions for these eigenvalues in Sections B.1–B.3 and show in Section B.4 that δp ∼ e^{sct}, where s_{c}≪s_{±}. The cooling time is therefore given by ${t}_{\mathrm{c}\mathrm{o}\mathrm{o}\mathrm{l}}^{1}={s}_{\mathrm{c}}$, while the other two eigenvalues correspond to the faster adjustment of the radiation fields. These conclusions remain unaltered when $\widehat{c}<c$ for high enough $\widehat{c}$. We study the applicability of the RSLA in Section B.5, where we use this analysis to obtain constraints to the values of $\widehat{c}$ that can be employed without introducing unphysical results.
B.1. Eigenvalues and cooling rate
We start by computing the eigenvalues of A, whose characteristic polynomial is
$$\begin{array}{cc}\hfill P(s)& =det(As{I}_{3})\hfill \\ \hfill & =(s{s}_{\mathrm{thin}})[(s+\frac{\widehat{c}}{{\lambda}_{P}})(s+\frac{\widehat{c}}{{\lambda}_{R}}+i\widehat{c}{\xi}^{\prime}k)\hfill \\ \hfill & +{\widehat{c}}^{2}{k}^{2}(\xi {\xi}^{\prime}f)]{s}_{\mathrm{thin}}\frac{\widehat{c}}{{\lambda}_{P}}(s+\frac{\widehat{c}}{{\lambda}_{R}}+i\widehat{c}{\xi}^{\prime}k)\phantom{\rule{0.166667em}{0ex}},\hfill \end{array}$$(B.4)
where I_{3} is the rank3 identity matrix and
$$\begin{array}{c}\hfill {s}_{\mathrm{thin}}=\frac{\eta c(\mathrm{\Gamma}1)}{{\lambda}_{P}}\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(B.5)
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
$$\begin{array}{c}\hfill s\ll \frac{\widehat{c}}{{\lambda}_{R}}\phantom{\rule{0.166667em}{0ex}},\end{array}$$(B.6)
which also guarantees $s\ll \widehat{c}/{\lambda}_{P}$, since our employed opacities verify ${\kappa}_{P}^{d}>{\kappa}_{R}^{d}$ for all temperatures (Melon Fuksman & Klahr 2022; this is also true in general for opacity laws of the form ${\kappa}_{\nu}^{d}\propto {\nu}^{n}$ as long as n > 0.81). In that case, the characteristic equation P(s)=0 loses its dependence on $\widehat{c}$ and becomes
$$\begin{array}{cc}& (s{s}_{\mathrm{thin}})[(1+{\xi}^{\prime}{\lambda}_{R}k)+{\lambda}_{P}{\lambda}_{R}{k}^{2}(\xi {\xi}^{\prime}f)]\hfill \\ \hfill & \phantom{\rule{2em}{0ex}}={s}_{\mathrm{thin}}(1+i{\xi}^{\prime}{\lambda}_{R}k)\phantom{\rule{0.166667em}{0ex}},\hfill \end{array}$$(B.7)
whose solution is
$$\begin{array}{c}\hfill s={s}_{\mathrm{thin}}\frac{(\xi {\xi}^{\prime}f){\lambda}_{R}{\lambda}_{P}{k}^{2}}{1+i{\xi}^{\prime}{\lambda}_{R}k+(\xi {\xi}^{\prime}f){\lambda}_{R}{\lambda}_{P}{k}^{2}}\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(B.8)
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 𝒪(λ_{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
$$\begin{array}{c}\hfill {s}_{c}={s}_{\mathrm{thin}}\phantom{\rule{0.166667em}{0ex}}\frac{{\lambda}_{R}{\lambda}_{P}{k}^{2}}{3+{\lambda}_{R}{\lambda}_{P}{k}^{2}}\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(B.9)
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 𝒪(λ_{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 frequencyintegrated radiative transfer equation, since ${P}_{r}^{ij}$ 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 t_{cool} = s_{c}^{−1}. Specifically,
$$\begin{array}{c}\hfill {t}_{\mathrm{cool}}=\frac{\rho \u03f5}{4\phantom{\rule{0.166667em}{0ex}}{a}_{R}{T}^{4}}\phantom{\rule{0.166667em}{0ex}}\frac{3+{\lambda}_{R}{\lambda}_{P}{k}^{2}}{{\lambda}_{R}{\lambda}_{P}{k}^{2}}\phantom{\rule{0.166667em}{0ex}}\frac{{\lambda}_{P}}{c}\phantom{\rule{0.166667em}{0ex}},\end{array}$$(B.10)
which is the constantopacity limit of Equation (16). The wavelengthdependent 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 dustgas 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 + 𝒪(λ_{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 scaleindependent. This is not verified for smaller k values if λ_{P}λ_{R}k^{2}/3 is not much larger than 1, in which case $(\eta \delta p\delta {E}_{r})\approx \eta \delta p\frac{{\lambda}_{R}{\lambda}_{P}{k}^{2}}{3+{\lambda}_{R}{\lambda}_{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
$$\begin{array}{cc}\hfill {s}_{c}& =\frac{c(\mathrm{\Gamma}1)}{{\lambda}_{P}(1+{\lambda}_{P}{\lambda}_{R}{k}^{2}/3)}\hfill \\ \hfill & \times [\frac{{\lambda}_{P}{\lambda}_{R}{k}^{2}}{3}(\eta +{b}_{P}\frac{{a}_{R}{T}^{4}{E}_{r}}{p})i{b}_{R}\frac{{\lambda}_{P}\mathbf{k}\xb7{\mathbf{F}}_{r}}{p}]\phantom{\rule{0.166667em}{0ex}},\hfill \end{array}$$(B.11)
where ${b}_{P}=\frac{dlog{\kappa}_{P}^{d}}{dlogT}$ and ${b}_{R}=\frac{dlog{\chi}_{R}^{d}}{dlogT}$. This expression only depends on b_{P} far from the midplane, where E_{r} ≠ a_{R}T^{4}. On the other hand, the rightmost term proportional to b_{R} only appears in the diffusion regime, in which case it introduces an oscillatory component. Since in that case ${\mathbf{F}}_{\mathbf{r}}=\frac{{\lambda}_{R}}{3}\mathrm{\nabla}{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 $\widehat{c}$
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 $\delta {F}_{r}^{x}$ that evolve proportionally to e^{sct + ikx}, which we denote as $\delta {\overline{E}}_{r}$ and $\delta {\overline{F}}_{r}^{x}$. Those components must also satisfy Eq. (B.2), which leads to
$$\begin{array}{cc}\hfill ({s}_{c}/\widehat{c}+1/{\lambda}_{P})\delta {\overline{E}}_{r}& =(\eta /{\lambda}_{P})\delta pik\delta {\overline{E}}_{r}\hfill \\ \hfill ({s}_{c}/\widehat{c}+1/{\lambda}_{R})\delta {\overline{F}}_{r}^{x}& =ik(\xi {\xi}^{\prime}f)\delta {\overline{E}}_{r}i{\xi}^{\prime}k\delta {\overline{F}}_{r}^{x}\hfill \end{array}$$(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, $\delta {\overline{E}}_{r}$ and $\delta {\overline{F}}_{r}^{x}$ evolve quasistatically, as their readjustment time is much shorter than s_{c}^{−1}, which explains why $\widehat{c}$ disappears from the equations. In particular, if f ≪ 1, this leads to the usual expression for the radiative flux in the diffusion limit, ${\mathbf{F}}_{r}=\frac{1}{3\rho {\chi}_{R}}\mathrm{\nabla}{E}_{r}$, whose perturbation in this analysis becomes $\delta {\overline{F}}_{r}^{x}=i({\lambda}_{R}k/3)\delta {\overline{E}}_{r}$.
Using the obtained expression for s_{c} (Eq. (B.9)), we can see that Eq. (B.6) is satisfied for all k as long as ${s}_{\mathrm{thin}}\ll \widehat{c}/{\lambda}_{R}$, which leads to the equivalent condition
$$\begin{array}{c}\hfill \eta (\mathrm{\Gamma}1)\frac{{\lambda}_{R}}{{\lambda}_{P}}\ll \frac{\widehat{c}}{c}\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(B.13)
For $\widehat{c}=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 lightcrossing 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 largek condition is satisfied independently from Eq. (B.13). To see this, we rearrange the characteristic equation as
$$\begin{array}{c}\hfill s={s}_{\mathrm{thin}}(1\phi (s))\phantom{\rule{0.166667em}{0ex}},\end{array}$$(B.14)
with
$$\begin{array}{c}\hfill \phi (s)={(1+\frac{s{\lambda}_{P}}{\widehat{c}}+\frac{(\xi {\xi}^{\prime}f){\lambda}_{R}{\lambda}_{P}{k}^{2}}{1+s{\lambda}_{R}/\widehat{c}+i{\xi}^{\prime}k{\lambda}_{R}})}^{1}\phantom{\rule{0.166667em}{0ex}},\end{array}$$(B.15)
from which we can see that s = s_{thin} is obtained as long as φ(s_{thin}) ≪ 1. We can then combine both conditions, which leads to the conclusion that Eq. (B.9) is valid provided
$$\begin{array}{c}\hfill min(\frac{\eta (\mathrm{\Gamma}1)({\lambda}_{R}/{\lambda}_{P})}{c/\widehat{c}},\phantom{\rule{0.166667em}{0ex}}\phi ({s}_{\mathrm{thin}}))\ll 1\phantom{\rule{0.166667em}{0ex}},\end{array}$$(B.16)
that is, as long as either the cooling timescale is much shorter than the lightcrossing time of a mean free path or the gasdust 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 $\widehat{c}$ values in Section B.5, except for $\widehat{c}/c={10}^{5}$ in the case f_{dg} = 10^{−3}, in which the lefthand side of Eq. (B.16) is of order 1 close to the midplane. We further study the effect of reducing $\widehat{c}$ 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
$$\begin{array}{c}\hfill (\frac{s{\lambda}_{P}}{\widehat{c}}+1)(\frac{s{\lambda}_{R}}{\widehat{c}}+1+i{\xi}^{\prime}{\lambda}_{R}k)+(\xi {\xi}^{\prime}f){\lambda}_{R}{\lambda}_{P}{k}^{2}=0\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(B.17)
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
$$\begin{array}{c}\hfill {s}_{+}\approx \frac{\widehat{c}}{{\lambda}_{P}},\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}{s}_{}\approx \frac{\widehat{c}}{{\lambda}_{R}}\phantom{\rule{0.166667em}{0ex}},\end{array}$$(B.18)
while for λ_{P}k ≫ 1, and therefore also λ_{R}k ≫ 1, we have
$$\begin{array}{c}\hfill {s}_{\pm}\approx \frac{\widehat{c}}{2}(\frac{1}{{\lambda}_{R}}+\frac{1}{{\lambda}_{P}}+i{\xi}^{\prime}k)\pm i\widehat{c}k\sqrt{\xi {\xi}^{\prime}f{({\xi}^{\prime})}^{2}/4}\phantom{\rule{0.166667em}{0ex}},\end{array}$$(B.19)
where the argument of the square root is always positive if f < 1. Therefore, s_{±} is at least of order $\sim max({\lambda}_{R}^{1}\widehat{c},k\widehat{c})$. 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 $\widehat{c}$. The same is true for λ_{R}k ≫ 1, in which case the resulting condition is $\widehat{c}/c\gg \eta (\mathrm{\Gamma}1)/({\lambda}_{P}k)$, which in that limit holds everywhere for $\widehat{c}=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 CayleyHamilton theorem, which allows us to compute e^{At} as a linear combination of {I, A, A^{2}} with timedependent coefficients, where I is the 3 × 3 identity matrix. This results in expressions of the form
$$\begin{array}{c}\hfill \delta p=\delta {p}_{0}[{e}^{{s}_{\mathrm{c}}t}+{\displaystyle \sum _{i}}{\mathcal{O}}_{i}\left(\eta \frac{c}{\widehat{c}}\right){e}^{{s}_{i}t}]\end{array}$$(B.20)
if λ_{R}k ≲ 1, and
$$\begin{array}{c}\hfill \delta p=\delta {p}_{0}[{e}^{{s}_{\mathrm{c}}t}+{\displaystyle \sum _{i}}{\mathcal{O}}_{i}\left(\frac{\eta}{{({\lambda}_{P}k)}^{2}}\frac{c}{\widehat{c}}\right){e}^{{s}_{i}t}]\phantom{\rule{0.166667em}{0ex}},\end{array}$$(B.21)
if λ_{P}k ≫ 1 , where 𝒪_{i}(x) denotes a term of order x. For $\widehat{c}=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}_{\mathrm{c}\mathrm{o}\mathrm{o}\mathrm{l}}^{1}={s}_{c}$. On the other hand, the relative size of these terms changes for $\widehat{c}<c$, which together with the changes in the eigenvalues can affect the overall relaxation rate for sufficiently small $\widehat{c}$. 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 $\widehat{c}/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 $\widehat{c}/c<1$ by numerically computing δp(t) in our domain, defining ${t}_{\mathrm{cool}}(\widehat{c})$ as the efolding time of δp, i.e., $\delta p({t}_{\mathrm{cool}}(\widehat{c}))/\delta p(0)={e}^{1}$. This definition coincides with the real cooling time when $c=\widehat{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 AlMohy & 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 $\widehat{c}/c$ values in the interval [10^{−5}, 10^{0}]. The obtained distribution of the relative error $\delta {t}_{\mathrm{cool}}^{1}(\widehat{c})/{t}_{\mathrm{cool}}^{1}(c)$, with $\delta {t}_{\mathrm{cool}}^{1}(\widehat{c})={t}_{\mathrm{cool}}^{1}(\widehat{c}){t}_{\mathrm{cool}}^{1}(c)$, is shown in Fig. B.1 for all considered $\widehat{c}/c$ and f_{dg} values. For the two highest $\widehat{c}/c$ values, we also show in Fig. B.2 the signed relative variation $\delta {t}_{\mathrm{cool}}^{1}(\widehat{c})/{t}_{\mathrm{cool}}^{1}(c)$, which shows that the cooling timescale does not always decrease with $\widehat{c}$, since the efolding time depends on the contribution of the different terms in Eqs. (B.20) and (B.21). The computed relative errors grow for decreasing $\widehat{c}$ and stay on the order of 1% and below for all parameters except for $(\widehat{c}/c,{f}_{\mathrm{dg}})=({10}^{5},{10}^{3})$, in which case values of 10% are reached in a significant portion of the domain covering the disk midplane. For those parameters, the obtained cooling times are always lower than those with $\widehat{c}=c$. We can then expect significant deviations in the resulting VSI evolution in that case.
Fig. B.1. Relative variation of the cooling rate with respect to its value when $c=\widehat{c}$ as a function of $\widehat{c}/c$ for f_{dg} = 10^{−3} (top row) and 10^{−4} (bottom row). 
Fig. B.2. Signed relative variation of the cooling rate with respect to its value with $c=\widehat{c}$ for $\widehat{c}/c={10}^{4},{10}^{5}$ and both considered f_{dg} values. 
An analytical constraint for the minimum value of $\widehat{c}/c$ that does not alter t_{cool} can be obtained by noting that, in the midplane region where the largest errors occur for $\widehat{c}/c={10}^{5}$, the most stringent condition for the value of $\widehat{c}/c$ among those derived in Sections B.2B.4 is given by Eq. (B.16), which guarantees that s_{c} is not affected by the reduction of $\widehat{c}$. This condition is verified except for $(\widehat{c}/c,{f}_{\mathrm{dg}})=({10}^{5},{10}^{3})$, in which case the lefthand side of the inequality is of order 1 close to the midplane, which is consistent with the error distributions in Fig. B.1.
We can now use these results to estimate the maximum possible reduction of $\widehat{c}$ that does not affect the dynamics of the VSI. In Skinner & Ostriker (2013), it was argued that the RSLA is applicable as long as all relevant timescale hierarchies are unaffected by such reduction. However, it is not enough for this analysis to compare the order of magnitude of t_{cool} and Ω^{−1} (or any other dynamical time) to guarantee the independence of the growth of VSI modes on $\widehat{c}$, since if t_{cool} > t_{crit}, the linear growth rate and the strength of the saturated instability depend on the precise value of t_{cool} (Lin & Youdin 2015; Manger & Klahr 2018; Klahr 2023). It is therefore necessary that t_{cool} is unchanged by the reduction of $\widehat{c}$, which is approximately the case if $\widehat{c}/c\ge {10}^{4}$, as already shown. It remains to be evaluated if $\widehat{c}$ is much larger than the maximum gas velocity, which in our runs is on average at most 7% of the isothermal sound speed. Since $\widehat{c}={10}^{4}c$ is always between 30 and 60 times larger than that speed, we can expect this value to be high enough not to introduce unphysical results. However, this can only be assessed through testing, as we do in Appendix C.
Appendix C: RSLA tests
Fig. C.1. Average kinetic energy (top) and maximum velocity (bottom) for different $\widehat{c}/c$ values at a resolution of 480 × 512. Solid and dashed lines correspond to f_{dg} = 10^{−3} and 10^{−4}, respectively. 
Fig. C.2. Same as Fig. 8 but comparing the Reynolds stress components for different $\widehat{c}/c$ values at a resolution of 480 × 512. 
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). 
To validate the nominal choice of $\widehat{c}/c={10}^{4}$ used in all simulations shown so far in this work, we tested the results of Appendix B by recomputing some quantities related to the evolution of the VSI for N_{θ} = 512 using different $\widehat{c}/c$ values spanning from 10^{−5} to 10^{−3}, with an additional run with $\widehat{c}/c={10}^{2}$ for f_{dg} = 10^{−3} (see Table 1). We studied differences in the time evolution of the VSI by computing the kinetic energy and the maximum velocity as a function of time as in Section 3.1 (Fig. C.1), and we evaluated changes in the angular momentum transport in the saturated stage by computing the components of the Reynolds stress tensor as in Section 4.1 (Fig. C.2). For f_{dg} = 10^{−4}, all curves match almost perfectly except for the random oscillations in v_{max}, as expected from the analysis in Appendix B. For f_{dg} = 10^{−3}, we do notice that, both during the growth of the finger modes and the body modes, the runs with $\widehat{c}/c={10}^{5}$ have a slightly smaller kinetic energy than the other runs, whereas the saturated steady state seems identical in all cases. The curves for all other $\widehat{c}/c$ values overlap except during the growth of the body modes, where, however, no identifiable trend as a function of $\widehat{c}/c$ is followed. A lower kinetic energy with $\widehat{c}/c={10}^{5}$ is also observed when comparing runs dg3c4_1024 and dg3c5_1024, with much smaller relative differences between the curves (not shown here). It must be noted that the value of $\widehat{c}/c$ not only may affect the cooling time, but it also determines the maximum time step used for the time evolution of radiation fields (Melon Fuksman et al. 2021), and thus some variations are also expected due to the accumulation of small numerical differences and possibly even different amounts of numerical diffusion. On the other hand, all curves overlap almost perfectly in the different saturated stages of the VSI. From these results, we conclude that our chosen value of $\widehat{c}/c={10}^{4}$ is high enough not to introduce unphysical effects in any stage of the VSI. On top of this, it even seems that $\widehat{c}/c={10}^{5}$ is enough to describe the steady state of the saturated VSI, and only introduces small variations in the growth rates. However, this result only applies to our employed setup, and similar testing is needed in general.
Appendix D: Wavelength estimation
In this appendix we describe our procedure for the estimation of the wavelengths used in the computation of radiative cooling timescales. We achieve this by computing the distance between two consecutive sign changes of v_{θ} at fixed θ. Given that wavelengths are typically larger above the irradiation surface than below (see, e.g., Fig. 3), we carry out this computation for two different heights on each side of the surface, namely θ − π/2 = 0.1 and 0.28. This is possible also in the VSIsuppressed region for f_{dg} = 10^{−4}, given that in that case we still obtain slow vertical modes below the irradiation surface (see Fig. 3). To improve the statistics of our estimation, we accumulate the λ values computed in this way every ∼0.8 orbits between t = 150 and 300 T_{5.5}. The obtained wavelengths are therefore the predominant ones after saturation, and not necessarily the fastest growing ones in the linear stage. 2D histograms of kH = 2πH/λ and radial distance to the star are shown in Fig. C.3, where we used a locally defined scale height $H={c}_{s}{\mathrm{\Omega}}_{K}^{1}$ with ${c}_{s}^{2}=p/\rho $.
At the highest resolution, maximum wavenumbers up to kH ∼ 100 occur in every case close to the inner radial boundary, where the resolution is maximal. At that resolution, we obtain typical values of kH between ∼20 and ∼80 below the irradiation surface, whereas above it values with kH ∼ 10 also occur. The dispersion of kH values grows in all cases with resolution and with distance from the midplane, since it is mainly caused by KH eddies. In Lin & Youdin (2015) it was proposed that the maximum wavenumbers, if not limited by numerical diffusion, are limited by the effective turbulent viscosity following the approximate relation (kH)^{2} ≲ qH/α, where α parameterizes the produced viscosity as ν = αc_{s}H. Even though the Reynolds stresses produced by the VSI are largely anisotropic (Section 4.1) we can test this prescription using α ≈ α_{r}, since the wave vector is mainly radial. For our simulation parameters and the obtained ⟨α_{r}⟩_{p} values, this gives the condition kH ≲ 70, which is similar to the obtained wavelengths. However, this analysis is not very precise, as larger values of up to ∼280 can be obtained if ⟨α_{R}⟩_{p} is used instead.
Fig. D.1. Estimated distributions of the parameter λ_{P}λ_{R}k^{2} for all resolutions and f_{dg} values, indicating the optical depth of the VSI modes (λ_{P}λ_{R}k^{2} ≫ 1 indicates optically thin cooling). 
As a side note, Fig. C.3 shows discontinuous wavenumber jumps in some cases, for which different explanations can be found in the literature. On the one hand, Stoll & Kley (2014) proposed that these can result from the lengthscale cutoff imposed on the VSI modes by the employment of a finitesized grid with finite resolution. This would be in principle supported by the variation in the discontinuous jump location as resolution is changed (Fig. C.3). On the other hand, the simulations in Svanberg et al. (2022) show wavelength jumps at fixed radial locations, which were proposed by the authors to correspond to turning points where a linear analysis predicts the radial wavenumber to go to zero. However, we have not seen signs of such turning points in r − t maps of the angular momentum at different heights at any of our resolutions (not shown here), although this could be due to the relatively smaller radial extent of our domain.
Since we are only interested in obtaining a rough estimate of k to be used in Eq. (B.10), we fit the measured wavelength distributions as a function of r using a thirdorder polynomial. The resulting k(r) functions, shown in Fig. C.3 as kH(r), are then taken as representative values of k in the entire regions above and below z/r ∼ 0.2, where the wavelength transition approximately occurs. A 2D distribution of k is then obtained by connecting the fitted values in both regions via a linear taper function in the region limited by z/r = 0.2 ± 0.02. As shown in Fig. D.1, the resulting λ_{R}λ_{P}k^{2} values are always much larger than 1 for all of our employed resolutions, meaning that temperature perturbations cool down at the optically thin rate even in vertically optically thick regions of the disk (see Fig. 1). This also justifies the rather approximate way in which we computed k, which introduces negligible errors in the estimated cooling times, given that the minimum estimated λ_{R}λ_{P}k^{2} values are 𝒪(10), and therefore $\frac{3+{\lambda}_{R}{\lambda}_{P}{k}^{2}}{{\lambda}_{R}{\lambda}_{P}{k}^{2}}\approx 1$
All Tables
All Figures
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}_{5.5}^{1}$). Dashed, white lines indicate the τ_{P}(T_{s})=1 surface for stellar photons computed using Planckaveraged opacities, solid white lines show the τ_{P} = 1 surfaces for vertically traveling infrared photons, and green dashdotted lines correspond to the z = ±H surfaces. 

In the text 
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 ${\mathrm{\Omega}}_{K}^{2}$. Bottom panel: squared vertical BruntVäisälä frequency for both f_{dg} values in units of ${\mathrm{\Omega}}_{K}^{2}$ (solid lines) compared to their values for vertically isothermal disks (dotted lines). 

In the text 
Fig. 3. Vertical flows produced by the VSI. Top row: vertical velocity distributions normalized by the isothermal sound speed ${c}_{\mathrm{s}}=\sqrt{p/\rho}$ in our highestresolution 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}. 

In the text 
Fig. 4. Rms velocities computed between 150 and 300 orbits for different resolutions and dust content, averaged at constant r (left) and constant θ (right). 

In the text 
Fig. 5. Time evolution of the volumeaveraged kinetic energy (top) and maximum gas velocity (bottom) in all runs with $\widehat{c}/c={10}^{4}$. 

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

In the text 
Fig. 7. Same as Fig. 4, but showing instead constantr and constantθ averages of the averaged absolute value of the vertical shear between 150 and 300 orbits. The corresponding values at t = 0 are shown for comparison. 

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

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

In the text 
Fig. 10. Pressureweighted zaverages of α_{R} computed for different resolutions and f_{dg} as in Eq. (10). 

In the text 
Fig. 11. α_{r} values for varying resolution and dust content computed as S_{ϕr}/⟨p⟩_{t}, averaged at constant r (left) and constant θ (right). 

In the text 
Fig. 12. Averages at constant θ of α_{z} computed as in Eq. (12) for varying resolution and f_{dg}. 

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

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

In the text 
Fig. 15. Relative temperature variation after 300 orbits at the maximum employed resolution. 

In the text 
Fig. 16. Midplane relative temperature perturbation at different times in run dg3c4_2048. 

In the text 
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. 

In the text 
Fig. 18. Radial average of the heating rate produced by the VSI in units of Ω for all employed resolutions and f_{dg} values. 

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

In the text 
Fig. A.1. Polar slice showing the specific angular momentum in the 3D disk simulation by Flock et al. (2020) after 350 orbits of evolution. 

In the text 
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. 

In the text 
Fig. B.1. Relative variation of the cooling rate with respect to its value when $c=\widehat{c}$ as a function of $\widehat{c}/c$ for f_{dg} = 10^{−3} (top row) and 10^{−4} (bottom row). 

In the text 
Fig. B.2. Signed relative variation of the cooling rate with respect to its value with $c=\widehat{c}$ for $\widehat{c}/c={10}^{4},{10}^{5}$ and both considered f_{dg} values. 

In the text 
Fig. C.1. Average kinetic energy (top) and maximum velocity (bottom) for different $\widehat{c}/c$ values at a resolution of 480 × 512. Solid and dashed lines correspond to f_{dg} = 10^{−3} and 10^{−4}, respectively. 

In the text 
Fig. C.2. Same as Fig. 8 but comparing the Reynolds stress components for different $\widehat{c}/c$ values at a resolution of 480 × 512. 

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

In the text 
Fig. D.1. Estimated distributions of the parameter λ_{P}λ_{R}k^{2} for all resolutions and f_{dg} values, indicating the optical depth of the VSI modes (λ_{P}λ_{R}k^{2} ≫ 1 indicates optically thin cooling). 

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