Thermal Structure Determines Kinematics: Vertical Shear Instability in Stellar Irradiated Protoplanetary Disks

Turbulence is crucial for protoplanetary disk dynamics, and Vertical Shear Instability (VSI) is a promising mechanism in outer disk regions to generate turbulence. We use Athena++ radiation module to study VSI in full and transition disks, accounting for radiation transport and stellar irradiation. We find that the thermal structure and cooling timescale significantly influence VSI behavior. The inner rim location and radial optical depth affect disk kinematics. Compared with previous vertically-isothermal simulations, our full disk and transition disks with small cavities have a superheated atmosphere and cool midplane with long cooling timescales, which suppresses the corrugation mode and the associated meridional circulation. This temperature structure also produces a strong vertical shear at $\mathrm{\tau_*}$ = 1, producing an outgoing flow layer at $\tau_*<1$ on top of an ingoing flow layer at $\tau_* \sim 1$. The midplane becomes less turbulent, while the surface becomes more turbulent with effective $\alpha$ reaching $\sim10^{-2}$ at $\tau_* \lesssim$1. This large surface stress drives significant surface accretion, producing substructures. Using temperature and cooling time measured/estimated from radiation-hydro simulations, we demonstrate that less computationally-intensive simulations incorporating simple orbital cooling can almost reproduce radiation-hydro results. By generating synthetic images, we find that substructures are more pronounced in disks with larger cavities. The higher velocity dispersion at the gap edge could also slow particle settling. Both properties are consistent with recent Near-IR and ALMA observations. Our simulations predict that regions with significant temperature changes are accompanied by significant velocity changes, which can be tested by ALMA kinematics/chemistry observations.


INTRODUCTION
Turbulence in protoplanetary disks plays significant roles in planet formation, such as determining mass accretion, angular momentum transport, and dust dynamics.In the disk midplane beyond 0.1 au, where the magnetorotational instability (MRI) is suppressed due to low ionization rates, turbulence can be generated by hydrodynamic instabilities (Turner et al. 2014;Armitage 2020;Lesur et al. 2023).One promising candidate among these instabilities is the vertical shear instability (VSI) Corresponding author: Shangjia Zhang zhangs17@unlv.nevada.edu(Nelson et al. 2013;Stoll & Kley 2014;Barker & Latter 2015;Umurhan et al. 2016;Lesur et al. 2023).The VSI is driven by the vertical differential rotation of the disk (dv ϕ /dZ ̸ =0).It occurs in baroclinic disks, characterized by non-parallel density and pressure gradients (Lesur et al. 2023;Klahr et al. 2023).Such a configuration is often satisfied in stellar irradiated protoplanetary disks, where the temperature decreases away from the central star in the radial direction.Most studies on the VSI assume a vertically constant temperature, which is thought to be a valid assumption near the disk midplane (e.g., Calvet et al. 1991;Chiang & Goldreich 1997).
A distinctive feature of the VSI is the presence of the corrugation mode in the meridional plane, which involves large-scale gas circulations in the vertical di-rection, while remaining confined in the radial direction (Nelson et al. 2013;Lyra & Umurhan 2019).As a result, anisotropic turbulence arises in the Z-ϕ and R-ϕ stresses, with the former typically reaching magnitudes of 10 −2 and the latter ranging from 10 −4 to 10 −3 (Stoll & Kley 2014, 2016;Stoll et al. 2017;Flock et al. 2017;Manger & Klahr 2018;Manger et al. 2020Manger et al. , 2021;;Pfeil & Klahr 2021).The amplitude of turbulence increases with the aspect ratio (i.e., temperature) of the disk (Manger et al. 2020(Manger et al. , 2021)).Simulations with mm-sized dust particles show that they can be stirred up very high above the midplane due to the strong turbulence in the vertical direction (Stoll et al. 2017;Flock et al. 2017Flock et al. , 2020;;Blanco et al. 2021;Dullemond et al. 2022).High resolution simulations show that strong shears between neighbouring bands of these meridional circulations can generate vortices in the R-Z plane (Flores-Rivera et al. 2020;Klahr et al. 2023;Melon Fuksman et al. 2023a).Threedimensional simulations also demonstrate they can generate vortices and zonal flows in the R-ϕ plane due to the Kelvin-Helmholtz instability (Flock et al. 2017;Manger & Klahr 2018;Flock et al. 2020;Blanco et al. 2021;Pfeil & Klahr 2021).These zonal flows can possibly explain some of the ubiquitous substructures observed by ALMA dust continuum observations (e.g., Andrews et al. 2018;Long et al. 2018;van der Marel et al. 2019;Cieza et al. 2021;Blanco et al. 2021;Andrews 2020;Bae et al. 2023).
The role of thermodynamics is crucial in determining whether the vertical shear instability (VSI) can occur (Lin & Youdin 2015;Lyra & Umurhan 2019;Lesur et al. 2023).Specifically, the global mode of VSI requires a fast cooling timescale, where the cooling time normalized to the orbital time should be less than a threshold (Lin & Youdin 2015), where Ω K is the Keplerian frequency, h/r is the disk aspect ratio, q is the radial temperature power-law index, and γ g is the adiabatic index.Analytical studies considering dust-gas coupling and dust evolution have identified regions where the VSI can operate (Malygin et al. 2017;Pfeil & Klahr 2019;Fukuhara et al. 2021).Disks with globally uniform cooling times larger than this critical cooling time do not develop VSI (Manger et al. 2021).However, short length-scale perturbations may still grow after evolving for a very long time (Klahr et al. 2023;Pfeil et al. 2023).In cases where the midplane has a long cooling timescale, yet the atmosphere has a short cooling time, studies using vertically isothermal simulations demonstrate the persistence of the VSI in the disk atmosphere (Pfeil & Klahr 2021;Fukuhara et al. 2023).The VSI in the unstable disk atmosphere can even penetrate the stable midplane, as long as the VSI-stable layer is less than two gas scale heights and the VSI-unstable layer is thicker than two gas scale heights (Fukuhara et al. 2023).Additionally, with equivalent short cooling times, radiation hydrodynamic simulations produce similar results to vertically isothermal simulations (Stoll & Kley 2014;Flock et al. 2017).
While VSI aligns with certain observational facets, such as a low level of R-ϕ turbulence, it also has tensions with some other facets.
• A possible mechanism for substructures.The (sub-)mm dust continuum emission traces ∼0.1-10 mm dust particles residing in the disk midplane, where substructures have been detected in a majority of observed disks when sufficient resolution is achieved (Bae et al. 2023).These substructures predominantly manifest as gaps and rings, although occasional arcs and spirals have been observed as well.VSI can generate density perturbations and vortices owing to the zonal flow led by the corrugation mode.However, these perturbation are typically small and requires very high resolution and sensitivity to detect them (Blanco et al. 2021).
• Overpredicted turbulence α Z .In the disk's vertical direction, edge-on and inclined disks exhibit remarkably thin dust emission layers, which can be translated to a low value of α Z /St (< 10 −2 , Pinte et al. 2016;Doi & Kataoka 2021;Villenave et al. 2020Villenave et al. , 2022;;Sierra & Lizano 2020;Ueda et al. 2020Ueda et al. , 2021)), where α Z is the turbulence in the Z-ϕ plane, and St is the Stokes number that characterizes the gas-dust coupling.An exception is that the inner ring of HD 163296, which has α Z /St > 1 (Doi & Kataoka 2021).Adopting a typical Stokes number, the α Z is estimated to be ≲ 10 −4 in most cases.However, VSI generates very strong turbulence in the vertical direction (∼10 −2 , e.g., Stoll & Kley 2014), with the correlation at neighbouring radii, which conflicts with many observations (Dullemond et al. 2022).
• Yet-to-be-detected VSI corrugation mode.Gas line emission observations demand higher sensitivities, yet data have been accumulating from recent and ongoing ALMA large programs such as, MAPS ( Öberg et al. 2021).With the aid of channel maps derived from these observations, we can measure the 3D velocity and temperature structure of the disk for a large sample of disks (e.g., Miotello et al. 2023;Pinte et al. 2023).Barraza-Alfaro et al. (2021) predicts that the alternating blueshiftedredshifted corrugation mode can be observed in the CO channel maps given very high spectral resolution (50 m s −1 at ALMA band 6).However, there is no firm detection of this pattern so far.
While disk thermodynamics plays a key role in VSI, most previous studies focused on simple thermodynamics such as locally isothermal or orbital cooling treatments.The vertical thermal structure is also underexplored.With a more self-consistent treatment, we can provide a more robust model and improve connections between observations and theory.
To that end, we employ a self-consistent radiationhydrodynamics approach with temperature-dependent DSHARP opacity (Birnstiel et al. 2018).Unlike previous simulations that utilized locally isothermal equations of state or flux-limited diffusion approximation with constant opacity, we utilize the Athena++ (Stone et al. 2020) implicit radiation module (Jiang et al. 2014;Jiang 2021).This module incorporates angle-dependent radiative transfer equations with implicit solvers to accurately model the disk radiation transport.The module can capture both optically thin and thick regimes and shadowing and beam crossing accurately.Additionally, we incorporate stellar irradiation using longcharacteristic ray tracing as a heating source.
Recently, Melon Fuksman et al. (2023b,a) independently study VSI in irradiated protoplanetary disks using M1 method (Melon Fuksman et al. 2021).While we focus on the outer disk beyond 20 au and they focus on 4-7 au in the inner disk, the results for our fiducial model are consistent with their dust depleted disk models, which show quiescent midplane and turbulent atmosphere.
The paper layout is the following.In Section 2, we discuss the numerical setup of the simulations.In Section 3, we present results of radiation hydrodynamic simula-tions, accompanied with a series of pure hydro simulations.In Section 4, we discuss the formation of substructures, and observational/modeling prospect.Finally, we conclude the paper in Section 5.

Disk Model Setup
In our study, we explored both full disks and transition disks with varying inner disk truncation radii or cavity sizes, denoted as r cav .Our simulations were performed in spherical polar coordinates r, θ, ϕ, while most of the disk structure was set up in cylindrical coordinates R, Z, ϕ.
The gas surface density profile follows a power-law with an exponential cutoff, consistent with viscous evolution models (Lynden-Bell & Pringle 1974;Hartmann et al. 1998) and observational constraints (e.g., Miotello et al. 2023).The gas surface density is given by: where Σ g,0 is the gas surface density at a reference radius of R 0 = 1 au.In our fiducial models, Σ g,0 is set to 178 g cm −2 , following Zhu et al. (2012).
Regarding the initial temperature structure, we assumed a radial power-law profile as the initial condition: where q is set to -0.5.The reference temperature T (R 0 ) is given by: Here, f accounts for the flaring of the disk (e.g., Chiang & Goldreich 1997;D'Alessio et al. 1998;Dullemond et al. 2001), and we used a value of f = 0.1 in our initial conditions.The stellar luminosity L * is assumed to be 1 L ⊙ , and σ b represents the Stefan-Boltzmann constant.
Then the hydrostatic equilibrium in the R − Z plane requires the initial density profile at the disk midplane to be (e.g., Nelson et al. 2013) where at the initial condition, the gas surface density and the volume density is related by The midplane radial density profile power-law index p can be related to the surface density power-law index r (Σ g ∝ R r ) and temperature power-law index q, by p = r − q/2 − 3/2.Since we adopted q = −0.5 and r = −1, p = -2.25.The temperature is used to calculate the gas scale height h = c s /Ω K , where c s = (P/ρ) 1/2 is the isothermal sound speed, and Ω K = (GM * /R 3 ) 1/2 is the Keplarian orbital frequency.We adopt M * = 1 M ⊙ .
In the vertical direction and the azimuthal velocity GM * /R (e.g., Nelson et al. 2013).We can see that the vertical shear rate dv ϕ /dZ is non-zero as long as q ̸ = 0.The other two velocity components v R and v Z are set to be zero at the initial condition.We do not consider the self-gravity of the disk in this paper, which should be a valid assumption for most of the Class II disks (Miotello et al. 2023).
In the initial condition (Equation 3), we also assumed the temperature to be vertically isothermal.This assumption is valid near the midplane.However, we will demonstrate that the quasi-steady state of our simulations exhibits a cool midplane and a superheated atmosphere which is consistent with classical analytical calculations (Calvet et al. 1991;Chiang & Goldreich 1997;D'Alessio et al. 1998), Monte Carlo radiative transfer calculations (Pascucci et al. 2004;Pinte et al. 2009), previous radiation hydrodynamic simulations (Flock et al. 2013(Flock et al. , 2017(Flock et al. , 2020;;Kuiper et al. 2010;Kuiper & Klessen 2013), and recent ALMA CO observations (Law et al. 2022(Law et al. , 2023)).We will also show that this vertical temperature gradient, together with the varying local orbital cooling time, is crucial for the gas kinematics.
The full disk corresponds to a value of r cav equal to 3 solar radii (3r ⊙ -rad), which represents the magnetic spherical accretion truncation radius (e.g., Hartmann et al. 2016).The total gas mass is approximately 0.01 solar masses (0.007 solar masses for 54au-rad).The transition disks have r cav values of 18 au (18au-rad, fiducial model) and 54 au (54au-rad).Additionally, we considered a case where the gas surface density is reduced to 1% of the fiducial value, i.e., Σ g,0 is set to 1.78 g cm −2 (18au-lowdens-rad).In this low-density disk scenario, the total gas mass is approximately 10 −4 solar masses.For transition disk with r cav = 54 au, we used a tanh profile to make a smooth transition at the inner gap so that it satisfies the Rayleigh criterion (e.g., Yang & Menou 2010) at the initial condition.We summarize all our models in Table 1.
Limited by the computational cost, the inner boundary of the disk cannot be too small.Thus, we set the inner boundary at 21.6 au.For r cav = 3 r ⊙ , and 18 au disks, the simulation inner boundaries are beyond the cavity sizes.To mimic the optical depth effect of these disks, we artificially added optical depth between 3 r ⊙ /18 au and the simulation's inner boundary at 21.6 au (see Equation A8 in Appendix A).However, this preset optical depth cannot adjust its vertical structure selfconsistently and can lead to discontinuities at the simulation's inner boundary.The direct irradiation on disk inner cavity is also related to the shadowing effect in Dullemond et al. (2001);Jang-Condell & Turner (2012, 2013); Siebenmorgen & Heymann (2012); Zhang et al. (2021b).Therefore, what occurs near the inner boundary might not be reliable.

Radiation Hydrodynamics
We introduce the numerical setup for the radiation hydrodynamic simulations using the frequency-integrated (gray) radiation module (Jiang 2021).We detail the implementation of the stellar irradiation and unit conversion in Appendix A.
We adopted the Courant-Friedrichs-Lewy (CFL) number to be 0.4, and used second order Van Leer time integrator (vl2), second order spatial reconstruc-tion, and HLLC Riemann Solver.We adopted adiabatic index γ g = 1.4.We discretized the radial direction into 1568 cells, logarithmically spaced from 0.54 to 8 times the reference radius (r 0 = 40 au, so 21.6 au to 320 au from inner and outer boundaries).The polar direction was divided into 1536 cells, covering a range from 0.383 to 2.76 radians (68 • above and below the midplane).For our fiducial model, this amounts to 45 cells per scale height at r 0 (40 au).For the hydro boundary conditions, we used outflow at the inner boundary, and copied initial conditions for outer, upper and lower boundaries.As for radiation boundary conditions, light beams can freely transport out of the domain.If the beam points inward the computational domain, the radiation is assumed to have the background temperature (10 K = 1.63×10 −3 T 0 , where T 0 is the temperature in code unit), which is a typical temperature of molecular clouds.We adopted periodic boundary condition in the azimuthal ϕ-direction.
The radiation transport uses discrete ordinate, where rays are discretized into different angles.We used the discretization better suited for curvilinear coordinates (angle flag = 1) and set nzeta = 2, npsi = 2, where nzeta represents angles from 0 to π/2 in ζ direction and npsi represents 0 to π in ψ direction.Here ζ and ψ are the polar and longitudinal angles with respect to the local coordinate, so these angles can point to different directions at different spatial locations (Jiang 2021).There are 16 angles in total.To test the convergence of the temperature calculated by different numbers of rays, we froze the hydrodynamics (assuming a static disk) and tried nzeta = 4, npsi = 4, and nzeta = 8, npsi = 8.We found convergence when nzeta = 4 and npsi = 4. Since the temperature difference is already small between nzeta = npsi = 2 and nzeta = npsi = 4, we adopted the former to save computational time.
We ran simulations with cfl rad = 0.3 (cfl rad is an additional factor multiplied in front of the CFL number to help convergence for implicit method), reduced speed of light R = 4 × 10 −3 (Zhang et al. 2018a;Zhu et al. 2020), and error limit = 10 −3 for the first 10 orbits (P in , the orbital period of the inner boundary at 0.54 r 0 or 21.6 au) to approach the quasi-steady state.Otherwise the iteration times or errors were extremely large.Then we restarted the simulation and changed cfl rad and R back to 1. Thus, we do not use reduced speed of light for the longer time evolution for our simulations.We also changed the error limit to 10 −5 after the restart.We note that after the restart it only took one iteration to reach the error limit and the typical error was only 10 −7 -10 −6 .

Dust Opacity Setup
While the radiation module can treat isotropic scattering properly, we neglected dust scattering and only considered absorption opacity in this paper to better compare with previous isothermal simulations.We also used the frequency-integrated (gray) radiation transport, but we note that multi-group radiation module is also available (Jiang 2022).Both dust scattering and multi-frequency radiative transfer will be considered in a future publication.
We used the DSHARP composition (Birnstiel et al. 2018) and a power law MRN dust size distribution (n(a) ∝ a −3.5 , Mathis et al. 1977).The minimum grain size a min = 0.1 µm and maximum grain size a max = 1 mm.In our fiducial models, we assumed that only small grains determine the temperature distribution due to their high opacity at the peak of the stellar spectrum; therefore, we considered grains sized between 0.1 and 1 µm, which account for f s =0.02184 of the total dust mass.The mass ratio between all the dust and gas was assumed to be 1/100.Then we calculated the Planck and Rosseland mean opacities normalized to the total dust mass (κ P,d , and κ R,d ) at various disk temperatures and fitted by univariate spline functions labeled as solid and dashed curves, respectively in Figure 1.We also calculated the stellar Rosseland mean opacity normalized to the total dust mass (κ * ,d = 3995 cm 2 g −1 ) at the solar temperature labeled as the star legend.These opacities can be simply converted to the ones normalized to gas (κ P,g , κ R,g , and κ * ,g ) by multiplying the dust to gas mass ratio, which is 0.01 for all models.Disk opacities are inputs for the radiation module, whereas the stellar opacity is used in the stellar irradiation as an extra heating source term (see Appendix A).

Comparison with RADMC-3D
We froze the hydrodynamics and used the initial conditions of full and transition disk models to test the temperature calculation of our irradiation implementation.Then we compared the results with Monte Carlo radiative transfer code RADMC-3D (Dullemond et al. 2012) using the same density structure.In RADMC-3D, we also used the same DSHARP opacity with the same dust properties and dust scattering turned off.When the disk is mostly optically thin to the stellar irradiation, the difference is at most ∼7%, such as in the transition disk with r cav = 54 au (shown in Figure 2).
The difference comes from the frequency integrated radiative transfer in Athena++ module and the more accurate multi-frequency treatment in RADMC-3D.This difference has been extensively studied in Kuiper et al. (2010); Kuiper & Klessen (2013), where they found that The temperature-dependent dust opacities adopted for all the radiation-hydro models.The solid line indicates the Planck opacity of the disk, whereas the dashed line indicates the Rosseland mean opacity of the disk.The wavelength-dependent DSHARP opacity (Birnstiel et al. 2018) is convolved at different temperatures to obtain temperature-dependent mean opacities.The stellar temperature is assumed to be 1 T⊙.Its Planck opacity and Rosseland mean opacity are assumed be the same and marked by the star.
the temperature calculated by the gray radiation transfer can be underestimated in the optically thick regime and overestimated in the optically thin regime for the stellar irradiation, which is consistent with our results.This is because the region near the midplane is optically thick to the stellar irradiation, so the heating comes from the τ * = 1 (τ * is the stellar optical depth integrated in the radial direction, see Equation A8) surface in the atmosphere (Calvet et al. 1991;Chiang & Goldreich 1997).However, even though the stellar spectrum peaks at optical to UV wavelengths, the continuum stellar spectrum still has a set of τ * = 1 surfaces for each frequency instead of a single one.The τ * = 1 surfaces at longer wavelengths can penetrate deeper and transport more energy to the midplane, thus increasing the temperature at the optically thick region (Kuiper et al. 2010).They also demonstrated that even within the single frequency radiation-hydro framework, the temperature calculation can be as accurate as the multi-frequency one by integrating the multi-frequency stellar irradiation to mimic continuous τ * = 1 surfaces ("hybrid method" therein).This treatment has also been implemented and tested for our problem and can be used in our future projects.For the full disk model (r cav = 3r ⊙ ) in the current paper, the temperature at the midplane can be underestimated by 40%, due to the very high optical depth.Therefore, processes that are sensitive to the absolute temperature values (e.g., chemistry) need the hybrid method (Kuiper et al. 2010;Kuiper & Klessen 2013) or the multi-band radiation module (Jiang 2022) in the optically thick regime.

Pure Hydro Simulations with Different Levels of Simplifications
Since most of the understanding on VSI was from previous vertically isothermal simulations and linear theory, we also ran various pure hydro simulations with different levels of simplifications to compare with our rad-hydro simulations.Namely, they are (a) vertically isothermal simulations with adiabatic EoS and instant cooling (β = 10 −6 ), (b) vertically varying background temperature with adiabatic EoS and instant cooling, and (c) varying background temperature with adiabatic EoS and local orbital cooling.Their model names are also listed in Table 1.
(a) Vertically Isothermal Simulations.We tried to compare the vertically isothermal simulations using the same disk aspect ratio (h/r) as the rad-hydro simulations, as the Reynolds stress is dependent on h/r (Manger et al. 2020) and also temperature power-law index q (Manger et al. 2021).Since radiation hydro simulation will adjust the temperature to reach hydrostatic equilibrium from the initial condition, the disk scale height can change from initial condition depending on the radial optical depth of the star.The midplane and atmosphere also have different temperatures.Thus, we measured the midplane temperature and surface density at r 0 (40 au), and radial power-law density and temperature indices (p and q) of the rad-hydro simulations and put them as initial conditions for these vertically isothermal simulations.The midplane temperatures for four radiation hydro simulations are shown in Figure 3, along with their initial condition (dashed line).
(b) Background Temperature with Isothermal EoS.We used the R-Z two dimensional background temperature averaged between t = 1000-1200 P in (500-700 P in for r cav = 54 au) from rad-hydro simulations and fixed them throughout the simulation.The gas density will adjust according to the temperature profile after the simulation begins.
(c) Background Temperature with Local Oribital Cooling.For these simulations, we used the R-Z background temperature as (b), but used adiabatic EoS with local orbital cooling, where the cooling means that the temperature will be relaxed to the background temperature in a dimensionless cooling time β (the cooling time normalized by the Keplerian orbital frequency).We used the simple optically thin and thick cooling times (Flock et al. 2017), where q is the midplane temperature power-law slope in the radial direction shown in Table 1.We assumed k = 10 as h/10 is a typical length scale measured in Lin & Youdin (2015).However, since the optically thin term (β thin ) dominates in most of the region, this adoption of h and k is not critical.By combining realistic vertical temperatures and location dependent cooling times, these simulations should have the closest thermodynamical properties compared with rad-hydro simulations, as we will demonstrate in the next section.

Overview
The significant difference between the rad-hydro simulation (top panel) and classical vertically isothermal simulation (bottom panel) can be demonstrated in Figure 4, where we show the line integral convolution1 (LIC, similar to Flores-Rivera et al. 2020) of our fiducial models (r cav = 18 au), 18au-rad and 18au-iso at t = 1000 P in , color-coded by the meridional velocity, v 2 R +v 2 Z 1/2 .In the vertically and locally isothermal simulation (bottom panel), the classical corrugation mode (the radially narrow, vertically extended circulation pattern) is clear.In the rad-hydro simulation (top panel), the disk is separated in two parts, the cool midplane and the superheated atmosphere.In the cool midplane, the velocity and turbulence levels are low, whereas the superheated atmosphere is more turbulent.The boundary between the cool midplane and superheated atmosphere exists a strong shear that leads to many small-scale vortices.The global circulation pattern that can be easily identified in locally and vertically isothermal simulation is replaced by turbulence on smaller scales.
Next, we analyze all of our radiation models in detail accompanied by pure-hydro simulations with various levels of simplifications.These quantities are taken either at t = 1000 P in or time-averaged values between 1000-1200 P in .In the case of the r cav = 54 au transition disk, the flow at the cavity edge will reach the critical condition for the Rayleigh stability criterion and become unstable.A giant vortex develops at ∼ 800 orbits, which should break into smaller vortices in realistic 3D disks.Therefore, we analyze this particular model from 500 to 700 P in (t=500 P in for the snapshot) to avoid this unphysical feature in 2D.

Thermal Structure Determines Kinematics
We will use four rad-hydro simulations to demonstrate that τ * = 1 surface (τ * : radial stellar optical depth, see Equation A8) sets disk temperature and equivalent local orbital cooling structures.Subsequently, these two structures determine the disk kinematics.
The 2D (R-Z) snapshots of gas density (left panels), temperature (middle panels), and the vertical velocity (right panels) for four models are shown in Figure 5.The gas density does not deviate significantly from the initial conditions.Gas scale heights are represented by the white contours overlaid on gas densities, calculated by assuming that vertical density follows a Gaussian profile (see Equation 6) and that the temperature is taken at the midplane.In cases where r cav = 3r ⊙ and 18 au (3r ⊙ -rad and 18au-rad), the midplane temperature is lower than the initial condition, whereas the atmosphere temperature is higher, resulting in smaller effective gas scale heights in the midplane (see Table 1 and Figure 3) and larger effective gas scale heights in the atmosphere.In contrast, for r cav = 54 au and 18 au, low-density cases (54au-rad and 18au-lowdens-rad), temperatures are higher than the initial condition, leading to larger effective gas scale heights (see Table 1).The black contours represent τ * = 1 surfaces where stellar irradiation intercepts the disk, setting the two temperature structure of the disk (cool midplane and superheated atmosphere).Note that the low density model (18au-lowdens-rad) lacks this surface, meaning that the entire disk is optically thin to stellar irradiation.
The locations of the τ * = 1 surfaces are shown in the temperature panels (middle panels) in Figure 5.We observe sharp decreases in temperature below these surfaces for the first three models.The low-density model (18au-lowdens-rad) is nearly vertically isothermal because the entire disk is optically thin to stellar irradiation.The white contours on top of the temperature maps indicate where the cooling time, estimated using Equation 8, equals the critical cooling time from Equation 1, denoted as β c in each panel.The radial temperature gradient q is estimated by fitting a power-law to the midplane temperature profile (see Table 1 and Figure 3).The r cav = 3 r ⊙ and 18 au models exhibit cooling times exceeding the critical cooling time in the cool midplane indicating stability, while the superheated atmosphere has cooling times less the critical cooling time, indicating that the region is unstable to the VSI.The r cav = 54 au and low-density models do not have cooling times exceeding β c , indicating instability throughout the domain.
The vertical velocity (v Z /c s ) structure in the right panels of Figure 5 also reflects temperature and cooling time structures.For the r cav = 3r ⊙ and 18 au models, the vertical velocity is below 1% of the local sound speed in the cool midplane and above 10% of the local sound speed in the superheated atmosphere.The separation occurs around a gas scale height and slightly below the τ * = 1 surface.In the atmosphere, the vertical velocity is still vertically extended and radially narrow, similar to the n=1 corrugation mode, but the upper and lower disk vertical velocities tend to have opposite signs, differing from the classical corrugation mode.We show in Appendix D (Figure 22) that these velocities tend to be anti-correlated.The classical n=1 corrugation mode only dominates when the disk has a vertically constant temperature and a short cooling time (β < β c ), which is the case inside the inner cavity of the r cav = 54 au model and throughout the low-density model.
For the full disk model (r cav = 3 r ⊙ ), we observe strong density and velocity perturbations near the inner boundary.This is a simulation artifact because our simulation domain cannot extend to 3 r ⊙ due to the computational cost associated with the very large dynamical range (see the end of Section 2.1).
To quantify the vertical structure of the disk in these radiative-hydrodynamic models, we present timeaveraged vertical profiles of gas density, temperature, the radially integrated stellar optical depth (τ * ), vertically integrated disk optical depth (τ P,d ), and cooling time (β) at r = 80 au in Figure 6.The densities of the r cav = 3 r ⊙ and 18 au models (3r ⊙ -rad and 18au-rad) exhibit two Gaussian distributions, one concentrated at the midplane and the other more extended in the atmosphere.These correspond to the cool midplane and the warmer atmosphere, as shown in the temperature profiles.The temperatures remain almost constant in these two regions, with the transition occurring between 0.1 to 0.3 radians, roughly equivalent to 2-4 gas scale heights.For the r cav = 54 au and low-density models (54au-rad and 18au-lowdens-rad), two-Gaussian profiles are not as clear, given their smoother temperature transition.The low-density model is nearly vertically isothermal, with a slight temperature drop at the midplane.The transition between optically thin and thick stellar irradiation occurs at approximately 0.17 radians for all three models at 80 au (indicated by the vertical dashed lines).The low-density case is optically thin to stellar irradiation.This model only contains 10 −4 M ⊙ , which falls at the lower end of protoplanetary disk masses.Another way to achieve a very optically thin disk is to modify our assumptions regarding the fiducial small grain fraction and the dust-to-gas mass ratio.If the disk is entirely depleted of small particles or has an extremely low dust-to-gas mass ratio, it can become optically thin to stellar irradiation.The τ P,d panel indicates that all models are optically thin to the disk's emission.The dimensionless cooling times for these models range from 10 −3 to 10. Their critical cooling times, denoted by the Figure 5. From left to right: the gas density (ρ), temperature (T), and vertical velocity (vZ ) of four radiation-hydro models at t = 1000 Pin (at t =500 Pin for rcav = 54 au, 54au-rad) in the meridional (R-Z) plane.From top to bottom: the full disk model with rcav = 3 r⊙ (3r⊙-rad), the transition disks with rcav = 18 au (18au-rad), rcav = 54 au (54au-rad), and 1% of the fiducial density (18au-lowdens-rad).The white contours on the density maps mark the one gas scale height.The black contours are the locations where the stellar optical depth in the radial direction reaches unity (the last model, 18au-lowdens-rad, has τ * < 1 for the whole disk).The βc on the temperature maps is the critical cooling time for VSI, represented by the white contours.The region enclosed by the contour near the midplane has β > βc, so the VSI should not be operating according to linear analysis.The last two models, 54au-rad and 18au-lowdens-rad, have β < βc for the whole domain.horizontal dashed lines, differ due to variations in gas scale heights and radial temperature gradients.For the r cav = 3 r ⊙ and 18 au models, the transition between the VSI unstable atmosphere and the VSI stable midplane occurs around 0.15 radians, approximately 2-3 gas scale heights.The r cav = 54 au and low-density models have cooling times smaller than the critical cooling times throughout their vertical extent, suggesting that the VSI should operate along the entire vertical extent.
To indicate locations where VSI is active or inactive in the R-Z plane, we present R-ϕ (T R,ϕ , left panels) and Z-ϕ (T Z,ϕ , middle panels) Reynolds stresses, normalized by the time-averaged pressure, along with root mean square velocity normalized by the local sound speed (right panels) for four radiation models in Fig- (Nelson et al. 2013;Stoll & Kley 2014;Flock et al. 2017).The values are calculated between t = 1000-1200 P in , except for the 54au-rad simulation which was taken between t = 500-700 P in .In all three columns, brighter colors indicate higher turbulence values.In the two stress columns, red colors denote positive values, while blue colors indicate negative values.The 3r ⊙ -rad model shows low levels of α R (≲ 10 −5 ) and α Z (≲ 10 −4 ) in the midplane, but they increase from the inner to the outer disk.The stress becomes much larger around the τ * = 1 surface and approaches 10 −2 before reversing sign to negative values.Overall, α Z is larger than α R near the midplane.The root mean square velocity is at the percent level of the sound speed near the midplane, approaching a fraction of the sound speed in the atmosphere, with the strongest values at the transition region (see Section 3.3).The 18au-rad model exhibits similar behavior.Its midplane turbulence becomes slightly higher and the low turbulence region is more confined to the midplane.The 54au-rad model shows VSI-like anisotropic turbulence between α R and α Z (similar to 18au-lowdens-rad) in the cavity and at the ring until ∼ 90 au, where α R is ∼ 10 −4 and α Z is ∼ 10 −2 .At the outer disk (100-160 au), turbulence levels become much higher (≳ 10 −2 ) in both stress components, which are much stronger than the midplane regions in all other models.This suggests that the gap edge in transition disks can be highly turbulent due to direct stellar irradiation.Beyond 160 au, however, the 18au-rad model has higher turbulence levels in the midplane, since the more turbulent atmosphere has enough energy to disturb the midplane due to lower density contrast at the outer disk (also see Figure 4).For the 18au-lowdens-rad model, the turbulence structures are very similar to those of vertically and locally isothermal VSI simulations.The α R is ∼ 10 −4 and α Z is ∼ 10 −2 throughout the disk, while the root mean square velocity is at the percent level of the local sound speed globally, yet still higher than the midplane values for 3r ⊙ -rad and 18au-rad models.

Accretion, Zonal Flow, and Vertical Shear Rate
The stress structure determines the disk's accretion structure.We can reveal this relation by averaging the angular momentum equation in the azimuthal direction (e.g., Turner et al. 2014;Lesur 2021;Rabago & Zhu 2021;Zhu et al. 2023) and obtain where Since our simulations are in 2D, we calculate time averaged instead of azimuthally averaged quantities.If we adopt a smooth disk structure (Equation 7), the change of ⟨v ϕ ⟩ is small in Z direction, so the second term in the first line can be neglected.We note that in our rad hydro models, the shear at the transition region between the atmosphere and midplane can be large, but since ⟨ρv Z ⟩ is still small compared to ⟨ρv R ⟩, the second term is still less than the first term.If we also assume the disk reaches a steady state (left hand side is zero), the accretion structure is only determined by the derivatives of the stresses.That is 2 , Then we can use Figure 7 and Equation 11 to explain the time-averaged radial velocity (v R /c s ) in Figure 8 shown on the left panels.Since in all our models the vertical gradient of the stress is greater than the radial gradient among the transition region between midplane and atmosphere, and T Z,ϕ is greater than T R,ϕ , we can just use the second term on the right hand side of Equation 11 and the center column of Figure 7 to understand the accretion structure.The radial velocity in the midplane is very small for both r cav = 3 r ⊙ and 18 au models due to the small stress and stress gradient there.At the boundary between the cool midplane and superheated atmosphere, an outgoing flow is on top of an ingoing flow, resembling layered accretion.In Figure 7, the region of the ingoing flow is aligned with the sharp transition of negative stress to positive stress above the midplane, whereas the region of the outgoing flow is aligned with the transition of positive stress to negative stress in the upper atmosphere.The sign and magnitude of ⟨ρv R ⟩ can be perfectly calculated from Equation 11.These regions also align with the regions that have the strongest shear (right panels).In contrast, the midplane shear rate is at least an order of magnitude smaller.The r cav = 54 au and low-density models have ingoing flow in the midplane, and outgoing flow in the atmosphere, which is consistent with previous analytical studies on VSI operating isothermal disks (Stoll et al. 2017;Rabago & Zhu 2021).Such an accretion structure is expected when the turbulence levels between R − ϕ and Z − ϕ directions are anisotropic, with the former smaller than the latter by around two orders of magnitude.These vertical shear rates for these two models are more uniformly distributed in the vertical direction, due to their more vertically uniform stress profiles which result from smoother temperature profiles.Compared to r cav = 3 r ⊙ and 18 au models, the shear rates in these two models are larger in the midplane.
The direction of the radial velocity and shear rate can be intuitively understood by examining the azimuthal velocity shown in the middle panels of Figure 8 (azimuthal velocity subtracted by the local Keplerian velocity).For the first two models that have temperature stratification, the midplane has faster rotational velocity than the atmosphere.At the transition region near the midplane, the gas loses its angular momentum due to the shear, thus moving inward, whereas the gas near the atmosphere gains its angular momentum and moves outward.In the bottom two models, which lack a strong vertical temperature structure, the azimuthal velocity has less vertical dependence, resulting in significantly smaller shear and radial velocity.Additionally, all four models exhibit zonal flows with perturbations to the azimuthal velocity.The r cav = 3 r ⊙ and 18 au models have relatively smaller scale perturbations on several au, whereas the 54 au model has zonal flows on tens of au.The low density model's azimuthal velocity perturbation is small and follows the corrugation pattern for vertical velocity.These zonal flows are associated with substructures in gas (see Section 4.1) and can be observed in near-infrared scattered light (see Section 4.2).
To be more quantitative, we present vertical profiles of various time-averaged velocity fields at r = 80 au in Figure 9. Regarding radial velocity, the outgoing and ingoing flows for the r cav = 3 r ⊙ and 18 au models are approximately 5% of the local sound speed, with the outgoing flow having a higher magnitude than the in-going flow.The r cav = 54 au model shows ingoing flow in the midplane and outgoing flow above the midplane, both at around 2% of the sound speed.In contrast, the low-density model displays ingoing flow in the midplane at less than 1% of the sound speed, and outgoing flow above the midplane, which gradually increases to 4% of the sound speed at 0.5 radians.The time-averaged vertical velocities in the first three models near the midplane are at the level of 1-2% of the sound speed, while the lowdensity model has negligible velocity.As for azimuthal velocity, all models show sub-Keplerian motion due to the radial pressure gradient.The r cav = 3r ⊙ and 18 au models exhibit a sharp transition in azimuthal velocity around 0.2 radians, attributed to the steep vertical temperature gradient between the cool midplane and the superheated atmosphere.In the same region, there is a smooth transition in azimuthal velocity for the r cav = 54 au model.The shear rate for the r cav = 3 r ⊙ and 18 au models peaks at 0.2 radians, reaching 0.15 Ω −1 K,0 , where Ω −1 K,0 is the time unit, or the inverse of Keplerian orbital frequency at unit radius r 0 = 40 au.The r cav = 54 au model also exhibits a slight increase in the same region.In contrast, the low-density model displays a linear increase in vertical shear above and below the midplane.The root mean square of the full velocity is also highest around 0.2 radians for the two models.At the midplane, the r cav = 3 r ⊙ model is relatively quiet, with velocities at approximately 2% of the sound speed.The 18 au model and the low-density model exhibit similar values at the midplane, approximately 10% of the sound speed.In contrast, the 54 au model is much more turbulent, with velocities exceeding 20% of the sound speed even at the midplane.
We also integrated ⟨ρv R ⟩ along the vertical direction to obtain (time-averaged and vertically integrated) radial mass accretion rates ( Ṁacc = 2πR ⟨ρv R ⟩dZ) as functions of R and show them in Figure 10 (upper panel).We also show the vertically integrated α R parameter in the lower panel, defined as Ṁacc is associated with the radial gradient of α int , since if we integrate Equation 11 along Z and assume that ⟨v ϕ ⟩ equals to the midplane Keplerian speed v K and does not change with Z, Equation 11 can be written as   where the T Z,ϕ is typically small at the boundaries due to the small densities at the disk surface.
The fluctuation is strong in the radial direction with the accretion rate frequently changing signs on au scale (Figure 10 upper panel), whereas the integrated α int has a lower variability (Figure 10 lower panel).For r cav = 3 r ⊙ model, the accretion rate is negative (ingoing) in the inner disk until 35 au, then it changes signs rapidly until reaching positive (outgoing) 10 −9 M ⊙ yr −1 in the outer disk.Except near the inner boundary, the α int increases from ≲ 10 −6 at 50 au to ∼ 10 −2 beyond 200 au.The r cav = 18 au model has positive (outgoing) 10 −9 M ⊙ yr −1 inside 25 au, but it can be affected by the setup for inner boundary.Then the accretion rate becomes negative until 100 au, and becomes positive in the outer disk.The α int increases from ∼ 10 −4 to ∼ 10 −2 from inner to outer disk.The r cav = 54 au model has the highest magnitude of accretion rate, reaching positive 10 −8 M ⊙ yr −1 around 100 au and positive 10 −8 M ⊙ yr −1 around 200 au.The low density model has almost zero accretion rate (< 10 −11 M ⊙ yr −1 ; the M acc for the low density model is multiplied by 100 to show its value) and α int between 10 −5 -10 −4 throughout the disk.We want to emphasize that even though the fiducial 18au-rad model has lower stress values than the equivalent vertically isothermal model 18au-rad-lowdens at the midplane, the integrated α int can be still larger due to the layered accretion at the disk surface.
We quantify the stress levels in Figure 11.These stresses exhibit fluctuations with respect to time, radius, and θ, thus we compute the average of these quantities between 60-100 au to represent the stress at 80 au.The time average spans 200 P in , consistent with previous figures.We present the stresses for three rad-hydro simulations (r cav = 3 r ⊙ , 18 au, and 54 au models in panels a, b and d), alongside three pure-hydro simulations with various assumptions (panels c, e, and f, which we will discuss more in the next section, Section 3.4).Solid lines represent α Z , while dashed lines represent α R .
For the full disk model (3r ⊙ -rad), both α R and α Z peak around 0.2 radians, coinciding with the location of the strongest shear and reaching values on the order of 10 −2 .In the midplane, α R remains less than 10 −6 , whereas α Z is around 10 −4 .The 18au-rad model exhibits similar behavior, with α Z and α R both reaching values of approximately 10 −2 at the transition region between the midplane and the atmosphere.For both components, α remains around 10 −4 to 10 −3 in the midplane.In the case of the 54au-rad model, vertical turbulence α Z consistently remains at a higher level, around 10 −3 to 0.4 radians, whereas α R can be lower by an order of magnitude, approximately 10 −3 .For classical locally- isothermal and vertically isothermal VSI (18au-iso), the anisotropy between these two components are more pronounced, with α Z (∼ 10 −2 ) exceeding α R (∼ 10 −4 ) by two orders of magnitude, consistent with previous studies (e.g., Stoll et al. 2017).
In the radial direction, not only the integrated turbulence α int increases with radius as shown in Figure 10, but also the midplane turbulence.This can be seen in 2D velocity maps such as Figures 4 and 5, and turbulence map in Figure 7.We quantify this by plotting the time-averaged v 2 Z in Figure 12.On the left panel we show four rad-hydro models.For 3r ⊙ -rad and 18au-rad, the vertical velocities are low in the inner disk ∼ 1% of the local sound speed at ∼ 50 au and increases to more than 10% of the sound speed in the outer disk at ∼ 200 au.For the 54au-rad model, the turbulence level is high across the disk (≳ 10% of the sound speed), where the perturbation is aligned with the zonal flow.The trend can be explained by the fact that as the effective scale height becomes larger in the outer disk, the density contrast between the midplane and atmosphere becomes smaller.At ∼ 200 au, the turbulent flow at the atmosphere already has enough energy to disturb the quiet midplane.For the 18au-lowdens-rad model, the turbulence is almost constant across the disk, on the order 10 1 10 2 10 3 10 4 10 5 0 10 5 10 4 10 3 10 2 10 1 3r -rad a) 18au-rad b) 18au-bkgT-bkgCool c) 0.5 0.0 0.5 /2 -[rad] 10 1 10 2 10 3 10 4 10 5 0 10 5 10 4 10 3 10 2 10 1 . The turbulence stresses for Z-ϕ and (solid lines) R-ϕ (dashed lines) components normalized by the averaged gas pressure at 80 au along the vertical direction.The values are calculated between t = 1000-1200 Pin, and averaged from r = 60-100 au.From left to right, they are 3r⊙-rad, 18au-rad, 18au-bkgT-bkgCool, 54au-rad, 18au-iso, and 18au-bkgT models, respectively.
of 10% of the sound speed, since there is no thermal stratification.
In a recent paper by Melon Fuksman et al. (2023b), their f dg = 10 −4 model (f dg is the dust to gas mass ratio) also has a similar temperature and cooling time stratification to our fiducial model, which leads to a quiet midplane and turbulent atmosphere.We find consistent results except a higher level of turblence in the midplane.Their midplane turbulence is below 10 −7 whereas ours is between 10 −5 -10 −4 .This difference might be explained by multiple factors.First, their inner disk has a lower h/r, and a lower h/r leads to a lower turbulence (Manger et al. 2020(Manger et al. , 2021)).In Figure 12, we also show that the turbulence increases with R since the energy contrast between the midplane and atmosphere becomes lower, which points to a lower turbulence value for the inner disk.They also have 200 cells per scale height resolution which can better resolve lower turbulence levels.

Good Approximation: Background Temperature with Local Orbital Cooling
As rad-hydro simulations differ significantly from classical VSI shown in vertically and locally isothermal simulations, we attempt to use pure-hydro simulations with varying levels of assumptions for comparison with radhydro simulations.We use Figure 13 to illustrate that pure hydro simulations with background temperature and local orbital cooling provide a good approximation for rad-hydro simulations.
Figure 13 displays (from left to right) vertical velocity snapshots, time-averaged radial velocity, azimuthal velocity subtracted from Keplerian velocity, and vertical shear rate.From top to bottom panels, we have the rad-hydro fiducial model (18au-rad), the vertically and locally isothermal model (18au-iso), the locally isothermal and background temperature model (18au-bkgT), and the orbital cooling and background temperature model (18au-bkgT-bkgCool).The rad-hydro model has Z /c 2 s that represents the turbulence level in the midplane for rad hydro models (3r⊙-rad, 18au-rad,54au-rad, and 18au-lowdens-rad) on the left panel and fiducial models (18au-rad, 18au-iso, 18au-bkgT, and 18au-bkgT-bkgCool) on the right panel.
been introduced in Figures 5 and 8; we retain them since these panels can be directly compared with the isothermal model in the second row, representing a classical VSI picture.The isothermal model (18au-iso) also resembles the low-density rad-hydro model.We note that the radial temperature gradient and h/r in 18au-iso model are measured in the midplane of rad-hydro simulations for a close comparison.In the pure-hydro simulation that incorporates the background temperature of the rad-hydro simulation (third row, 18au-bkgT), we observe the disruption of the delicate n=1 corrugation mode.Layered accretion and strong shear in the temperature transition region become evident, but the midplane vertical velocity remains higher than in rad-hydro simulations, as this region is still VSI unstable due to the almost zero cooling time (β = 10 −6 ).The zonal flow in the azimuthal velocity also differs from the rad-hydro simulation.
By incorporating the estimated cooling time from the rad-hydro simulation in the bottom panels (18au-bkgT-bkgCool), all four fields closely resemble the rad-hydro simulations.Therefore, we demonstrate that computationally inexpensive pure-hydro simulations can capture crucial features of rad-hydro simulations.In future, one can prescribe a temperature structure or obtain a self-consistent temperature structure by iterating the temperature calculated from Monte Carlo Radiative Transfer code such as RADMC-3D (Dullemond et al. 2012) and enforcing vertically hydrostatic equilibrium (Bae et al. 2019;Ueda et al. 2019;Zhang et al. 2021b).Then the cooling structure can be estimated using Equation 8.
It is important to note that they are not identical; an apparent difference is that the zonal flow in the purehydro model has narrower length scales than in the radhydro simulations.This discrepancy could be attributed to the static temperature and estimated cooling profiles in pure-hydro simulations, while rad-hydro simulations undergo secular evolution of temperature and cooling times.Additionally, the cooling is no longer local, since the radiative cooling can be affected by other parts of the disk.Their Reynolds stresses are similar but not identical.
In Figure 11, the 18au-bkgT-bkgCool and 18au-bkgT models exhibit shapes more similar to the rad-hydro model (18au-rad) but still display some differences.Their α Z and α R exhibit similar magnitudes, indicating that the turbulence becomes much more isotropic than in the isothermal model (18au-iso).However, 18au-bkgT-bkgCool shows lower turbulence levels, and the sign of Z-ϕ stress differs from 18au-rad beyond 0.3 radians and within 0.1 radians.Additionally, α R is smaller in the midplane.The α Z of the 18au-bkgT model has a different sign than 18au-rad between 0.2-0.3radians.Overall, values from rad-hydro models and those pure-hydro models that adopted the rad-hydro thermal structures have smaller than 10 −3 α Z in the midplane, except in the case of a large cavity (54 au).They also exhibit more isotropic turbulence than the isothermal one.The low-density case is not presented here, but its stress profile is similar to that of an isothermal simulation with a larger h/r than the 18au-iso  model due to its almost isothermal profiles and low equivalent cooling time (Figures 5 to Figure 8).
In the radial direction, the turbulence levels between the pure hydro and rad-hydro simulations are similar but not identical.On the right panel of Figure 12, we show pure hydro models along with the fiducial radhydro model.18au-bkgT and 18au-bkgT-bkgCool models follow the fiducial model's trend, but the former over-estimates the turbulence whereas the latter underestimate the turbulence.The local orbital cooling prescription strongly under-predict the turbulence level within 70 au.The vertically isothermal model has almost constant turbulence level around the disk similar to that of the low density rad-hydro model.

Gas Substructures
Related to the zonal flows shown in Figure 13, gas substructures can also develop depending on the inner cavity size as shown in Figure 14.The full disk model (3r ⊙ -rad) preserves the initial condition except at the inner disk due to the boundary effect.The 18au-rad model has perturbation on several au scale, whereas the 54au-rad model has perturbation on tens of au scale.18au-bkgT-bkgCool and 18au-bkgT models have similar perturbations as the 18au-rad, whereas the isothermal model 18au-iso keeps the initial condition.
Figure 15 shows the time evolution of the surface density.The full disk model (3r ⊙ -rad) and the isothermal model (18au-iso) do not show evident substructures.For the rest of the models, substructures can form and propagate to the outer disk.for r cav = 18 au models, some of these rings form at the inner boundary, but other rings can form in the middle of the disk and move to the outer disk at a lower speed.The time evolutions between 18au-rad, 18au-bkgT, and 18au-bkgT-bkgCool are not identical.The rings in 18au-bkgT-bkgCool model move at a lower speed than those in 18au-rad.For the 54au-rad model, the perturbation becomes much stronger after 700 orbits due to the vortex formation around 100 au in R-θ plane.However, in 3D simulations, vortices tend to develop in R-ϕ plane.Future 3D studies will unveil a more realistic structure of this model.These gas substructures can also possibly lead to dust substructures, but this needs to be tested in future studies that includes dust particles.Overall, our limited sample of rad-hydro models shows a tentative trend that transition disks with larger cavity sizes are more prone to develop zonal flows and substructures.

Observational and Modeling Prospect
While we focus on 2D modeling, the thermal structure is calculated self-consistently so we can make some predictions for axis-symmetric disks.In Figure 16 we show the near-infrared scattered light polarized intensity at λ=1.63 µm (H-band) using RADMC-3D for faceon 54au-rad, 18au-rad, and 18au-iso models in linear scale from 0 to maximum value.We used the same DSHARP opacity (Figure 1), assuming that the dust to gas ratio is 0.01 and small grains (0.1-1 µm) account for 0.02184 of the total dust mass.The temperature is directly taken from rad-hydro simulations.For the isothermal model (18au-iso), we run thermal Monte Carlo (radmc3d mctherm) to calculate temperature in the r − θ plane.The RADMC-3D calculated the full Stokes image using the scattering matrix.The polarized intensity is (Q 2 + U 2 ) 1/2 .For the large cavity transition disk model 54au-rad (taken at t = 500 P in ), we can clearly observe the inner rim and rings at around 100 au.If we take the snapshot at a later time, the ring structure becomes more evident, consistent with the surface density perturbation in Figure 15.However, we need to test whether this perturbation is still large in 3D simulations in future.For 18au-rad model, several rings can also be seen close to the inner rim, resembling Figure 14, but the length scale is smaller than the disk with a larger cavity, making their substructures more difficult to be observed.In contrast, we can only see the inner rim of the vertically and locally isothermal model (18au-iso), meaning that the outer disk is in the shadow (i.e., this is a self-shadowed disk as defined in Garufi et al. 2018Garufi et al. , 2022)).While we only have limited numbers of rad-hydro models with varying cavity sizes, we find the tendency that disks with larger cavity sizes can produce wider rings and can be easier to be observed in near-infrared scattered light images.This is consistent with the finding in current scattered light disk demographics which shows that ring structures are predominantly found in disks with weak Near-IR excess (Benisty et al. 2023), where weak Near-IR excess is often interpreted as no inner disk.
As shown in Figure 12, the turbulence level in the vertical direction appears to be stronger in the region where the disk is directly exposed to stellar irradiation.For example, 54au-rad and 18au-rad-lowdens models that have less attenuation of the stellar irradiation have higher turbulence values than the 3 ⊙ -rad and 18au-rad models.This suggests that dust turbulent diffusion could be also strong at the cavity edge.In HD 163296, Rosotti et al. (2020); Doi & Kataoka (2021, 2023) have found that the α/St value is higher for the inner ring at 68 au, which is more exposed to stellar irradiation than the outer ring at 101 au.This finding aligns with our results, assuming that the Stokes number (St) does not vary significantly between these two rings.To validate the impact of direct stellar irradiation on dust diffusion, it will be crucial to measure dust settling at the cavity edges through ALMA observations and incorporate dust particles into 3D hydro simulations.Furthermore, in Figure 12, for full disks or transition disks with small cavities, midplane turbulence values increase with radius, a trend that can be readily examined through ALMA observations.
Different molecular lines from the ALMA MAPS Large Program ( Öberg et al. 2021) and other highresolution, high-sensitivity datasets are used to map the temperature structure in protoplanetary disks (Zhang et al. 2021a;Law et al. 2021Law et al. , 2022Law et al. , 2023Law et al. , 2024)).Conversely, these data can also reveal the velocity vec-  2023) identify disk winds and meridional flows attributed to embedded planets.Currently, we anticipate a substantial increase in sample size from the ALMA Large Program exoALMA.Aligning with the theme of the current paper-where temperature structure influences kinematics-we use Figure 17 to illustrate our ability to establish correspondence between the temperature structure and kinematic features.The figure shows temperature contours overlaid by meridional velocity vectors for the 18au-rad, 54au-rad, and 18au-lowdens-rad models.In our fiducial 18au-rad model, the temperature is 10-20 K in the midplane, and the velocity remains low within 200 au.At the transition region from the midplane to the atmosphere, the gas flows inward, reaching 10% of the sound speed.At higher altitudes, the gas flows outward with increasing velocity towards the outer disk.In contrast, for the lowdensity model, since the temperature has no vertical dependence, the gas velocity remains low throughout the disk.The change of the flow structure in response to the temperature structure can be identified in observations.However, the emission surface does not necessarily parallel the streamline.For example, if the emission surface height increases with radius first and then decreases due to a lower optical depth (e.g., Figure 4 in Galloway-Sprietsma et al. 2023), changing directions of the velocity measured along the emission surface can be simply attributed to the multiple crossings of the streamlines, without the need of any radial substructure, such as disk wind or meridional flow.A recent study by Martire et al. (2024) has demonstrated that the difference in inferred rotational velocities from 12 CO and 13 CO can be attributed to vertical thermal stratification (see middle panels of Figure 8), with the midplane (traced by 13 CO) exhibiting a higher rotational velocity than the atmosphere (traced by the more optically thick 12 CO).By accounting for thermal stratification, they can retrieve properties such as disk mass and stellar mass more accurately.In our current paper, we want to highlight that considering thermal structure can lead to further consequences, including spatially varying radial velocity, vertical velocity, and zonal flow.
0.5 1.0 1.5 g (t) / g (t = 0) 0.5 1.0 1.5 g (t) / g (t = 0) 0.5 1.0 1.5 g (t) / g (t = 0) The middle panel of Figure 17 shows the temperature and flow structure of the transition disk 54au-rad.In this case, the temperature varies both vertically and radially.At the cavity edge (Figure 3, 50-100 au), the radial temperature gradient is steeper than that in a smooth disk.A giant clockwise rotating vortex is also present.While the existence and strength of the vortex needs to be tested in 3D simulations, it could result from the radial and vertical variation of the temperature instead of an artifact from the rad-hydro simulation, since a similar feature was also found in a pure-hydro simulation, 54au-bkgT-bkgCool.Future efforts will focus on validating the existence and possible forming mechanism of this vortex through dedicated 3D simulations.We note that the temperature variation across the transition disk cavity has been identified from observations (Leemker et al. 2022) and can shape the variety of tran-sition disk gas substructures (Wölfer et al. 2023), meaning that the correspondence between the temperature structure and the kinematic features in transition disk can also be probed with current observations.It has been well-known that temperature structure strongly shapes disk chemistry, but we also aim to use Figure 17 to demonstrate that the thermal structure can influence disk chemistry by shaping disk kinematics, thereby affecting material transport.While static disk models cover more complete chemical networks, studies considering dynamical effects (e.g., Aikawa & Herbst 1999;Semenov & Wiebe 2011;Furuya et al. 2013;Furuya & Aikawa 2014;Price et al. 2020;Bergner & Ciesla 2021;Van Clepper et al. 2022) by incorporating radial and/or vertical gas/dust mixing often reveal different chemical distributions compared to static models.These dynamic models may provide a better explanation for observations (see reviews by Krijt et al. 2022;Öberg et al. 2023, and references therein).We anticipate that, by accounting for advection and turbulent diffusion due to a specific thermal structure, the chemical distribution can be more accurately predicted.For instance, the layered accretion in our fiducial model 18au-rad, where an outgoing flow is atop an ingoing flow, can transport material inward in the colder layer and outward in the hotter layer, which may also have implications on solid transport in our solar system (Ciesla 2009).The turbulence level measured at both the midplane and the atmosphere in our fiducial model has some similarities to layered accretion models featuring an MRI-active atmosphere and a dead zone in the midplane (Gammie 1996;Simon et al. 2011Simon et al. , 2013;;Bai 2015;Xu et al. 2017;Simon et al. 2018).On the other hand, the high turbulence values in the upper layer are only measured in a few of protoplanetary disks (Flaherty et al. 2020;Paneque-Carreño et al. 2023) (see a recent review by Rosotti 2023).Our objective is to conduct 3D simulations to predict unique channel map features, following previous studies by Hall et al. (2020); Barraza-Alfaro et al. (2021).As demonstrated in Section 3.4 (Figure 13), background temperature and local orbital cooling profiles prove to be effective approximations for rad-hydro simulations.This sets the groundwork for using 3D simulations in our future work to make ALMA and nearinfrared scattered light observational predictions.

CONCLUSIONS
Vertical shear instability (VSI) is a promising candidate to generate turbulence in the outer region of protoplanetary disks.It can be crucial for gas and dust transport in protoplanetary disks.We study VSI using the Athena++ radiation module with stellar irradiation, which self-consistently captured the thermal structure and hydrodynamics.We study disks with different inner cavity sizes, accompanied by pure hydro simulations with various assumptions.We find that temperature structure strongly influences disk kinematics.Our main findings are as follows: 1.The radial optical depth of the star determines the disk's thermal structure.For realistic disk setups (M d = 10 −2 M ⊙ , cavity size < 54 au), the disk can be separated into the cool midplane and super-heated atmosphere, delineated by the τ * = 1 surface (Figure 5).If the disk is optically thin to the stellar irradiation (low mass disk with M d = 10 −4 M ⊙ ), the temperature is almost vertically isothermal.2. The thermal structure determines disk's kinematics.The temperature and cooling time (β) stratification suppresses the classical n = 1 corrugation mode that leads to meridional circulations found in isothermal simulations (Figure 4).Instead, the turbulence becomes more isotropic on a more local scale, in contrast to very large vertical-azimuthal Reynolds stress α Z (∼ 10 −2 ) and weak radialazimuthal Reynolds stress α R (∼ 10 −4 ) found in isothermal simulations (Figure 11).The low mass disk has vertically isothermal profiles, so it closely resembles all features in isothermal simulations.
3. The strongest vertical shear occurs at the transition region between the cool midplane and superheated atmosphere where many vortices form.At this transition region, layered accretion happens with an outgoing flow on top of an ingoing flow (Figure 8).This layered accretion can be perfectly explained by the vertical variation of the stress structure (Figure 7) using Equation 11.
4. Pure hydro simulations with measured temperature structures and estimated orbital cooling profiles can be good approximations for rad-hydro simulations (Figure 13).
5. Zonal flows and gas substructures can develop, and a disk with a larger cavity size has perturbations with a longer length scale and stronger magnitude (Figures 8 and 14).At the cavity edge, the gas has stronger turbulence, which could slow dust settling.Using MCRT simulations, we confirm that transition disks tend to have rings, and the disks with larger cavities tend to have more prominent rings, which are easier to be observed in nearinfrared scattered light images (Figure 16), consistent with the fact that rings in scattered light im-ages are predominantly found in disks with weak Near-IR excess.
We also show the correspondence between the temperature structure and kinematic features (Figure 17).Future work including synthetic observations can be used to predict this correspondence in ALMA observations.The temperature structure can also influence chemistry through shaping the disk kinematics.One can also predict the observational signatures in ALMA chemistry observations by modeling disk chemistry with disk dynamics under such a thermal structure.
The set of equations is Here, ρ represents the gas density, v is the flow velocity, P denotes the gas pressure, I is the unit tensor, E is the total energy, I represents the lab-frame specific intensity of photons emitted by the disk locally, c is the speed of light, and n is the angle in the lab frame.The terms S r (P ) and S r (E) represent the disk's radiation source terms in the momentum and energy equations.They are moments of the source term S I in the radiation transfer equation, given by: where κ a , κ s , and κ P are the Rosseland mean opacity, scattering opacity, and Planck mean opacity, respectively.These opacities are all normalized to the gas.The intensity I in the lab frame is related to the intensity in the co-moving frame I 0 through a Lorentz transformation: where γ ≡ 1/ 1 − v 2 /c 2 is the Lorentz factor, Γ(n, v) ≡ γ (1 − n • v/c), and n ′ is the angle in the co-moving frame given by: The angular-averaged mean intensity in the co-moving frame J 0 is defined as: where Ω 0 represents the angular element in the co-moving frame.The total gas energy density E is given by: where E g represents the gas internal energy.Assuming an ideal gas equation of state (EoS), the internal energy is related to the gas pressure P through the adiabatic index γ g as E g = P/(γ g − 1) for γ g ̸ = 1.The gas temperature T is calculated using T = µP/ (R ideal ρ), where R ideal is the ideal gas constant and µ is the mean molecular weight.We adopted γ g = 1.4 and µ = 2.3 in this paper.We adopted temperature, density and length units to be T 0 , ρ 0 , r 0 respectively.The time unit is given by Ω −1 K,0 = (GM * r −3 0 ) −1/2 .The velocity unit v 0 is then the Keplerian velocity at r 0 .These parameters are used to calculate two key parameters in the radiation module P ≡ a r T 4 0 / (ρ 0 R ideal T 0 /µ) and C ≡ c/a 0 (Jiang et al. 2012).a 0 is the characteristic isothermal sound speed, R ideal T 0 /µ 1/2 .The values for our fiducial model are T 0 = 6.14 × 10 3 K; ρ 0 = 4.28 × 10 −14 g cm −3 ; r 0 = 40 au; hence P = 1.13 × 10 3 ; and C = 6.36 × 10 4 .P represents the ratio between the radiation pressure and gas pressure at these unit quantities, whereas C represents the ratio between speed of light and characteristic sound speed a 0 .The P is larger than unity since we adopt a very large value of T 0 .We adopt this exact value of T 0 since v 0 = a 0 so we do not need to distinguish between the unit velocity and the characteristic sound speed.For typical values of density (∼ ρ 0 ) and temperature (tens of Kelvins) at the midplane, the radiation pressure is much less than the gas pressure.Similarly, the typical ratio between the speed of light and local sound speed is ≫ C, as T 0 is much greater than a typical disk temperature.We set density floor to be 10 −12 and pressure floor to be 10 −15 in code units.We also set a temperature floor to be 0.001 T 0 (6.14 K) and a temperature ceiling to be 0.1 T 0 (614 K) to avoid numerical hotspots.
To account for stellar irradiation, we include a heating source term in the energy equation.This source term is necessary for frequency-integrated radiation transport (e.g., Flock et al. 2017) as the stellar irradiation is at significantly higher temperatures (thousands of Kelvins) compared to the thermal emission from the disk (tens to hundreds of Kelvins).In a future work (Baronett et al., in prep) that uses a multi-group radiation module (Jiang 2022), this source term is not required, and it can better capture the multi-frequency nature of radiation transport both for the stellar irradiation and disk emission.
The stellar irradiation heating flux F * (r) is given by: where T * and R * represent the stellar surface temperature and radius, respectively.Here, σ b denotes the Stefan-Boltzmann constant, which is related to the radiation constant a r as σ b = a r c/4.The radial optical depth for the star at each θ is given by: τ * (r, θ) = where r in represents the inner radius of the computational domain.The first term in the second line of Equation A8 refers to the optical depth within the region interior to the computational domain.These values are not evolved with time but depend on the density and opacity setup of the global disk.Namely, they depend on the inner cavities' radii, gas scale height, and surface density.
To ensure compatibility with MPI (Message Passing Interface), where the ray-tracing needs to cross all the grids in the radial direction and a ray can enter different MeshBlocks located on different CPUs, we adopted the following procedure.First, we calculated the optical depth within each MeshBlock.Then, we declared a user-defined Mesh data array to store all the τ * values at the outer boundaries of each MeshBlock.These values represent the local optical depths integrated from the inner boundary (r mb,in ) to the outer boundary (r mb,out ) of each MeshBlock.For the zeroth column, it stores values of τ * ,bc .In the middle of each timestep, we cumulatively sum the the user defined Mesh data in the radial direction.Then an MPI Allreduce operation is performed to update all the user-defined Mesh data array in all the CPUs so that the inner boundary optical depths of each MeshBlock has their correct global values.Finally, the optical depth τ * within each MeshBlock can be calculated by adding up its current MeshBlock's inner boundary value and its local integrated value.Specifically, within each MeshBlock, the optical depth is given by: where τ * (r mb,in , θ) is the global optical depth at the MeshBlock's inner boundary stored in the user-defined Mesh data array.

B. ENERGY BUDGET
We use the energy equation in Equation A1 and refer to Figure 18 to illustrate the energy budget in our fiducial model, 18au-rad, where each term is time-averaged between t=1000-1200 P in .Assuming a steady state (∂E/∂t=0), the energy flux divergence term on the left-hand side (∇ • [(E + P )v]) should be balanced by three terms on the right-hand side: cooling from the disk radiation (-Sr(E)), heating from the stellar irradiation (-∇ • F * ), and the work done by the stellar gravity (ρa grav • v).We can move the gravity term to the left-hand side and incorporate it into the energy flux divergence, considering it as the flux divergence of gravitational potential energy.
Figure 18 compares the disk cooling (left panel), stellar heating (middle panel), and the energy flux divergence (right panel).In the atmosphere, stellar heating is prominent, and most of the energy is radiated away by disk cooling, as both -∇ • F * and Sr(E) exhibit similar values.However, in the midplane, stellar heating is nearly negligible, given our implementation of single-frequency stellar irradiation (Equation A7), where few photons can penetrate below the τ * = 1 surface.Interestingly, while disk cooling values are relatively small in the midplane, they still surpass stellar heating significantly.This additional disk cooling is offset by the negative energy flux divergence shown in the right panel.
Upon closer examination of each flux component, we find that the negative divergence arises from the advection of energy from the upper and lower atmospheres to the midplane.When we deactivate advection by freezing the hydrodynamics and solely conduct radiative transfer, the midplane temperature experiences a slight decrease by a few percent, indicating that energy advection can increase the temperature of the cold midplane by a few percent.We expect that when we include multi-wavelength stellar irradiation, the midplane can receive more stellar heating so that the influence of advection on the midplane temperature can be weaker.

C. DISRUPTION OF THE INERTIAL WAVE
Linear theory has unveiled the global model of VSI as an overstability, a destabilized inertial wave (Barker & Latter 2015).Recently, Svanberg et al. (2022) (also see Stoll & Kley 2014) used locally and vertically isothermal simulations to investigate the inertial wave patterns of VSI.A significant finding is that inertial waves associated with the corrugation mode could be identified in several radial wave zones separated at Lindblad resonances, each characterized by different frequencies.These wave zones also exhibit slightly different turbulence values.
However, when considering a self-consistent thermal structure that takes into account stellar irradiation, the inertial wave patterns become less distinct.Figure 19 illustrates the time evolution of the vertical velocity at the midplane, while the accompanying frequency and wavelength analyses are presented in Figures 20 and 21.In the isothermal simulation 18au-iso, we observe the classical corrugation model pattern as found in Stoll & Kley (2014) and Svanberg et al. (2022), characterized by alternating peaks and troughs moving radially.These patterns represent group velocity and phase velocities propagating in opposite directions (Figure 3 in Svanberg et al. 2022).In contrast, the remaining simulations do not exhibit this evident wave feature, except for 54au-rad between 30-60 au, within the cavity and at the ring location.This is consistent with our observations in Figure 5, where the corrugation mode is identified in v Z for 54au-rad near the cavity.Similarly, the 18au-lowdens-rad model (not displayed here) exhibits a wave pattern similar to the isothermal simulation.For the other models, we can still observe certain wave patterns, albeit different from the inertial wave patterns identified in the isothermal model.The pronounced velocity in 3r ⊙ -rad is likely an inner boundary effect due to the disk surface's (h/r) discontinuity.A comprehensive understanding of these wave patterns still needs ongoing studies.In summary, the classical VSI inertial wave pattern only manifests when the disk is optically thin to stellar irradiation so that the disk is close to vertically isothermal and has a short cooling time.This can occur when the disk has low dust opacity (18au-rad,lowdens) or when the disk possesses a wide inner cavity (54au-rad).

D. ASYMMETRY ABOVE AND BELOW THE MIDPLANE
While we have demonstrated that the inertial wave pattern associated with the corrugation mode in classical VSI becomes less apparent in rad-hydro simulations or pure hydro simulations with self-consistent thermal structures, we can still observe some elongated stripes in the vertical direction for v Z , as depicted in Figures 5 and 13.In contrast to isothermal simulations, the gas parcels above and below the midplane do not appear to move consistently in the same direction.The vertical velocity in Melon Fuksman et al. ( 2023b) also shows anti-symmetry when the disk midplane is VSI-inactive (f dg =10 −4 therein), whereas corrugation mode only occurs when the VSI is active in the whole domain (f dg =10 −3 therein).In Figure 22, we attempt to quantify these trends by measuring the autocorrelation function (dashed lines) for a horizontal cut within the range of 0.2-0.25 radians in the radial direction and a cross-correlation function (solid lines) between this cut and the one on the other side of the disk (π −θ) at a specific snapshot and then averaging over 200 P in .Their values are normalized by the autocorrelation value with no radial shift (radial shift ratio = 1).Absolute values between 0.1 to 1 are presented in a log-scale, while absolute values between 0 and 0.1 are shown in a linear scale.The cross-correlation of the isothermal model (18au-iso) closely aligns with the autocorrelation function, suggesting that the upper and lower disks move in the same direction, echoing the n=1 corrugation model in Figure 13.In contrast, the remaining models exhibit trends where autocorrelation and cross-correlation functions have the opposite signs during the first few turnovers, indicating that the upper and lower disks are more likely to move in opposite directions.This anti-correlation between upper and lower surfaces resembles that of the n=2 breathing mode found in the initial growing phase of VSI (e.g., Nelson et al. 2013;Barker & Latter 2015).However, this anti-correlation is not as strong as those in the linear growth phase, indicating that more than one modes are operating in this highly non-linear regime.A possible explanation is that the n=1 corrugation mode requires communication between the upper and lower disk.However, in our fiducial radiation models, the communication between two surfaces are disturbed by the temperature and cooling time stratification so that n=2 and other modes take over.Reynolds stresses in R-ϕ, Z-ϕ directions αR, αZ = T R,ϕ /¡P ¿, T Z,ϕ /¡P ¿, turbulence levels .The cross-correlation functions of the vZ between the upper and lower atmosphere shifted in the radial (r) direction with varying shift ratios, averaged between 0.2-0.25 radians above and below the midplane for t = 1000-1200 Pin (t = 500-700 Pin for 54au-rad).The values are normalized by the auto-correlation functions in the upper atmosphere without a shift (shift ratio = 1).If the upper and lower atmosphere has the same vZ, i.e., vZ(r, θ) = vZ(r, π -θ), the correlation should be unity at shift ratio = 1.
Figure1.The temperature-dependent dust opacities adopted for all the radiation-hydro models.The solid line indicates the Planck opacity of the disk, whereas the dashed line indicates the Rosseland mean opacity of the disk.The wavelength-dependent DSHARP opacity(Birnstiel et al. 2018) is convolved at different temperatures to obtain temperature-dependent mean opacities.The stellar temperature is assumed to be 1 T⊙.Its Planck opacity and Rosseland mean opacity are assumed be the same and marked by the star.

Figure 2 .
Figure 2. Temperature comparison between Athena++ and RADMC-3D for the transition disk model (54au-rad) at r = 80 au in the vertical direction.The solid and dashed curves show the temperatures in θ direction for Athena++ and RADMC-3D, respectively.

Figure 4 .
Figure 4.The line integral convolution (LIC) of the velocity field in the meridional plane, (vR, vZ ), color-coded by its magnitude, vmag for fiducial (rcav = 18 au) radiation-hydro (top, 18au-rad) and vertically isothermal (bottom, 18au-iso) models.The flow pattern in the bottom panel is very similar to Figure 7 in Flores-Rivera et al. (2020).

Figure 6 .
Figure6.Time-averaged (t = 1000-1200 Pin and 500-700 Pin for rcav = 54 au, 54au-rad model) gas density, temperature, stellar optical depth in the radial direction (τ * ), disk optical depth in the vertical direction (τ P,d ), and the cooling time (βc) for four radiation-hydro models cut at r = 80 au in the vertical direction.Models for 3r⊙-rad, 18au-rad, 54au-rad, 18au-lowdens-rad are shown in green, orange, purple, and magenta lines, respectively.The horizontal dashed lines in the last panel indicate the critical cooling time, βc.The vertical dashed lines indicate the τ * = 1 surface for the first three models.
Figure 7. From left to right: the turbulence stresses for Z-ϕ and R-ϕ components normalized by the averaged gas pressure, and the root mean square velocity normalized by the local sound speed.The values are calculated between t = 1000-1200 Pin, (at t =500-700 Pin for rcav = 54 au, 54au-rad).Other layouts are similar to Figure 5.

Figure 8 .
Figure 8. From left to right: the time-averaged (t = 1000-1200 Pin and 500-700 Pin for rcav = 54 au, 54au-rad model) values of the radial velocity (vR), azimuthal velocity subtracted by Keplerian velocity (v ϕ -vK), and the vertical shear rate dv ϕ /dZ for four radiation-hydro models in the meridional (R-Z) plane in the same layout as Figure 5.

Figure 10 .
Figure10.The time-averaged accretion rate in the radial direction (integrated in the vertical direction), and αint (vertically integrated αR) for four radiation-hydro models.The M⊙ yr −1 for the low density model is multiplied by 100 to show its value.The color representations are the same as previous figures.

Figure 14 .
Figure14.The gas surface density profiles in the radial direction at t = 1000 Pin or 500 Pin for 54au-rad (more opaque lines) and at the initial condition (more transparent lines).The layout is the same as Figure11.torsat line emission surfaces, as indicated in numerous studies (see a review byPinte et al. 2023).Specifically,Teague et al. (2019);Yu et al. (2021);Galloway- Sprietsma et al. (2023) identify disk winds and meridional flows attributed to embedded planets.Currently, we anticipate a substantial increase in sample size from the ALMA Large Program exoALMA.Aligning with the theme of the current paper-where temperature structure influences kinematics-we use Figure17to illustrate our ability to establish correspondence between the temperature structure and kinematic features.The figure shows temperature contours overlaid by meridional velocity vectors for the 18au-rad, 54au-rad, and 18au-lowdens-rad models.In our fiducial 18au-rad model, the temperature is 10-20 K in the midplane, and the velocity remains low within 200 au.At the transition region from the midplane to the atmosphere, the gas flows inward, reaching 10% of the sound speed.At higher altitudes, the gas flows outward with increasing velocity towards the outer disk.In contrast, for the lowdensity model, since the temperature has no vertical dependence, the gas velocity remains low throughout the disk.The change of the flow structure in response to the

Figure 15 .
Figure 15.The gas surface density normalized by the initial condition as a time evolution from t = 0-1200 Pin. Brown colors mean that the surface density is almost unchanged.Yellow colors indicate the increase of density, whereas blue colors indicate the decrease of density.The layout is the same as Figure 11.

Figure 17 .
Figure 17.Time-averaged temperature and meridional velocity vectors for 18au-rad, 54au-rad, and 18au-lowdens-rad models.10% of the local sound speed is shown on the upper left of each panel.

Figure 19 .
Figure 19.The time evolution of the vertical velocity (vZ) in the midplane from 0-1200 Pin for various models in the same layout as previous figures.

Figure 20 .
Figure20.The Fourier transform of the time evolution of the vertical velocity in the midplane (Figure19).The frequency (ω) is in unit of 1/ΩK,0, where ΩK,0 is the Keplerian frequency at 40 au.R is also normalized by r0 = 40 au.The diagonal solid line indicates ω = Ω = ΩK,0(R/r0) −1.5 , and the dashed line indicates ω = 0.1Ω.The layout is the same as previous figures.

Figure 21 .
Figure21.The wavelength occurrence of the vertical velocity (vZ ) in the midplane measured from every snapshot from t = 0-1200 Pin, with 1 Pin as the interval.The wavelength is measured as the distance of the neighbouring zero crossing points.The brighter colors represent higher counts.The layout is the same as previous figures.
Figure22.The cross-correlation functions of the vZ between the upper and lower atmosphere shifted in the radial (r) direction with varying shift ratios, averaged between 0.2-0.25 radians above and below the midplane for t = 1000-1200 Pin (t = 500-700 Pin for 54au-rad).The values are normalized by the auto-correlation functions in the upper atmosphere without a shift (shift ratio = 1).If the upper and lower atmosphere has the same vZ, i.e., vZ(r, θ) = vZ(r, π -θ), the correlation should be unity at shift ratio = 1.

Table 1 .
Simulation setups and measured values at later times.h/randq are measured from midplane temperature in Figure3.

Table 2 .
Parameters used in the paper stellar optical depth in the radial directionκ P,d , κ R,dPlanck and Rosseland mean opacities (to dust) κP,g, κR,gPlanck and Rosseland mean opacities (to gas) τ P,d , τ R,d disk optical depth in the vertical direction