Issue |
A&A
Volume 682, February 2024
|
|
---|---|---|
Article Number | A139 | |
Number of page(s) | 26 | |
Section | Astrophysical processes | |
DOI | https://doi.org/10.1051/0004-6361/202346554 | |
Published online | 16 February 2024 |
Vertical shear instability in two-moment radiation-hydrodynamical simulations of irradiated protoplanetary disks
I. Angular momentum transport and turbulent heating
Max Planck Institute for Astronomy, Königstuhl 17, 69117 Heidelberg, Germany
e-mail: 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 VSI-stability criteria.
Methods. We modeled the gas-dust mixture in our disk models by means of high-resolution axisymmetric radiation-hydrodynamical simulations including stellar irradiation with frequency-dependent opacities. This is the first study of this instability to rely on two-moment radiative transfer methods. Not only are these able to handle transport in both the optically thin and thick limits, but also they can be integrated via implicit-explicit methods, thus avoiding the employment of expensive global matrix solvers. This is done at the cost of 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 dust-to-gas mass ratio of 10% of our nominal value of 10−3 for < 0.25 μm grains), the VSI can operate in surface disk layers while being inactive close to the midplane, resulting in a suppression of the VSI body modes. The VSI reduces the initial vertical shear in bands of approximately uniform 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 VSI-induced angular momentum mixing in the radial direction may be predominantly axisymmetric. The stability regions in our models are well explained by local stability criteria, while the employment of global criteria is still 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 dust-depleted models, respectively, peaking at a few (approximately three) scale heights above the midplane. We conclude that it is unlikely that the VSI can, in general, lead to any 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, large-scale structures produced by different instabilities, such as zonal flows and vortices, are able to concentrate dust grains, possibly triggering the streaming instability and leading to the formation of planetesimals (Johansen et al. 2007; Gerbig et al. 2020; Schäfer et al. 2020). In poorly ionized “dead” zones of protoplanetary disks, where the magnetorotational instability (MRI, e.g., Balbus & Hawley 1991) is expected to be suppressed by nonideal phenomena such as Ohmic and ambipolar diffusion and, depending on the orientation of large-scale magnetic fields, the Hall effect (Bai 2015; Simon et al. 2018), purely HD instabilities have been predicted to take over as the main mechanisms driving turbulence (e.g., Lyra & Umurhan 2019; Pfeil & Klahr 2019). Among these, the vertical shear instability (VSI) stands out as a rather robust phenomenon dominating the production of turbulence in numerous HD (Nelson et al. 2013; Manger & Klahr 2018; Flock et al. 2017b) and even nonideal MHD simulations (e.g., Cui & Bai 2022).
Local and vertically global stability analyses (e.g., Urpin & Brandenburg 1998; Urpin 2003; Lin & Youdin 2015) predict that the VSI can be triggered as long as two conditions are met. On the one hand, the disk’s pressure and density stratification must be baroclinic (∇p × ∇ρ ≠ 0, where p and ρ are the gas pressure and density, respectively), which in vertical and radial hydrostatic equilibrium naturally leads to an azimuthal velocity stratification varying with distance from the midplane, that is, with vertical shear. Baroclinic stratifications naturally result from temperature gradients either in the radial or in the vertical direction, which are typically enforced by the equilibrium between stellar or turbulent heating and radiative cooling. On the other hand, buoyancy forces, which can restore gas perturbations into equilibrium, must be suppressed. Under these conditions, a gas parcel moving in a direction between its corresponding constant specific angular momentum surface and the disk rotation axis is accelerated away from is original position due to its angular momentum excess (e.g., Knobloch & Spruit 1982). Globally, this process can destabilize radially traveling inertial waves (Nelson et al. 2013; Barker & Latter 2015), leading to the production of vertically elongated vertical flows reaching significant fractions of the sound speed (e.g., Flock et al. 2017b).
If the additional buoyancy produced by magnetic tension (e.g., Latter & Papaloizou 2018) and dust back-reaction (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 scale-independent, the temperature perturbations produced by the VSI flows may span several characteristic sizes with different optical depths exceeding the typical VSI wavelengths. On top of this, if larger-scale structures such as vortices are produced by secondary instabilities or independent mechanisms (Marcus et al. 2015; Barranco et al. 2018; Manger & Klahr 2018; Pfeil & Klahr 2021), their optical depths may largely exceed that of the VSI modes, and hence they may cool at a different rate. Thus, the issue of computing accurate cooling timescales in HD simulations is only avoided if these naturally result from an effective solution of the radiative transfer equation, as it is done in radiation hydrodynamics (Rad-HD) codes. Moreover, unlike for isothermal and β-cooling prescriptions, such codes can be used to quantify the turbulent heating produced in VSI-unstable disk regions, which is particularly relevant for planet formation models relying on α-viscous heating prescriptions (e.g., Burn et al. 2022; Voelkel et al. 2022).
In this work we investigated the nonlinear evolution of the VSI in axisymmetric Rad-HD simulations of protoplanetary disks including heating by stellar irradiation with realistic tabulated opacities, focusing on the VSI growth and saturation phases, the produced transport of angular momentum and heating, and the applicability of local stability criteria. Both the baroclinic stratification and the thermal relaxation timescales required to trigger the VSI in our simulations arise from the balance of irradiation heating and radiative cooling. We studied a small disk region between 4 and 7 au around a T Tauri star and consider different degrees of depletion of small dust grains due to coagulation and settling, making the disk either vertically thick or thin and resulting in the formation different stability regions.
We neglected departures from thermal equilibrium between dust and gas species occurring in low-density regions, where the thermal relaxation timescale depends on a combination of the radiative cooling timescale and the collisional timescales of dust and gas species (see, e.g., Malygin et al. 2017; Barranco et al. 2018). This effect is considered in an extension of our stability analysis presented in a second part of this work (Melon Fuksman et al. 2024, Paper II henceforth).
We employed resolutions of up to ∼200 cells per scale height, which are unprecedented in Rad-HD VSI simulations, aiming to investigate the role of secondary instabilities on the VSI saturation in Paper II. The resulting computational overhead is limited by artificially reducing the speed of light, which in this context can be done as long as thermal relaxation timescales are unaffected by this operation, as we verified both analytically and numerically in this work. Also advantageous in terms of computational efficiency is the good parallel scalability of our employed numerical code (Melon Fuksman et al. 2021), which is the first Rad-HD method applied to the study of disk instabilities not to rely on globally implicit solvers.
This article is organized as follows. In 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 large-scale 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 gas-dust mixture in a circumstellar disk and its heating and cooling processes by means of the Rad-HD module (Melon Fuksman et al. 2021) implemented in the open-source PLUTO code (version 4.4, Mignone et al. 2007). In the chosen configuration for this work, the following system of equations is solved:
where ρ, p, and v are the gas density, pressure and velocity, while Er, Fr, 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 Er and Fr are measured in energy density units. The gas energy density is defined as imposing an ideal equation of state to the internal energy density, namely , 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 Fr and Er as , where the Eddington tensor Dij is defined as , with , where n = Fr/||Fr||, f = ||Fr||/Er, 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 Ms located at the center of coordinates, that is, Φ = −MsG/r. We include stellar heating by computing after each HD step the stellar irradiation flux produced by the central star as
where Ts is the star temperature, Rs is the star radius, Bν(T)=(2hν3/c2)/(ehν/kBT − 1) is the Planck radiative intensity at a temperature T and frequency ν, and h and kB 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 , where μ = 2.35 and u are respectively the assumed gas mean molecular weight and the atomic mass unit. The frequency-dependent 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 Rad-HD simulations, we retrieve it from the initial conditions.
Optical depths are computed considering only opacities due to dust absorption, which at optical and near-infrared wavelengths is dominated by small (≲1 μm) grains. This does not mean that there is not significant mass in larger grain sizes, but that larger sizes do not contribute to the total opacity, as specified in the next section. We neglect dust settling of small grains and assumed a perfect coupling between the dust and gas distributions with a constant dust-to-gas mass ratio fdg of small grains, thus computing the dust density of small grains as ρd = fdg ρ. 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.5da with radii in the range a ∈ [5, 250] nm. Opacity values are logarithmically sampled in the frequency range [νmin, νmax]=[1.5 × 1011, 6 × 1015] Hz. We also use these opacities to compute the frequency-integrated radiation-matter interaction terms, which in the gas comoving frame have expressions
where tildes indicate comoving quantities, while ar = 4σSB/c is the radiation constant, σSB is the Stefan-Boltzmann constant, and and are, respectively, the Planck- and Rosseland-averaged absorption and extinction coefficients. The laboratory frame coefficients (G0, G) are obtained via a Lorentz transformation of , which introduces terms of orders v/c and v2/c2 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 , which is the Rosseland-averaged absorption opacity, although the derivations in Appendix B consider a general . Graphs of the frequency-dependent and frequency-averaged opacities used in this work are shown in Melon Fuksman & Klahr (2022).
Our resolution method relies on the reduced speed of light approximation (RSLA), in which the speed of light in the radiation block is reduced to an artificially low value in order to reduce the scale disparity between gas and radiation characteristic speeds. This avoids the employment of the prohibitively small time steps that would be given by the Courant-Friedrichs-Lewy condition using the real speed of light, thus reducing the overall computational cost. The obtained solutions of the Rad-HD equations are in general unchanged for large enough values, as we later verify in our setup.
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, Er, Fr) 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 titer of 0.1 orbits at 1 au, which is enough to prevent artefacts induced by self-shadowing (Melon Fuksman & Klahr 2022) and reach convergence.
We constructed models of protoplanetary disks around a T Tauri star of mass Ms = 0.5 M⊙, effective temperature Ts = 4000 K, and radius Rs = 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 Nr × Nθ = 240 × 400. The total gas mass in that region is thus ∼0.04 M⊙. If a vertically integrated dust-to-gas mass ratio of approximately 0.01 is assumed, the total dust mass is approximately on the upper end of the distribution of estimated dust masses as a function of stellar masses for disks observed with ALMA (Manara et al. 2018).
The only physical parameter that varies among our simulations is the assumed dust-to-gas ratio of small grains, fdg, whereas all other physical quantities are kept constant. We considered two values of fdg: a nominal case with fdg = 10−3 and a dust-depleted case with fdg = 10−4. If we assume a total dust-to-gas mass ratio of 10−2 including all grain sizes, these values correspond respectively to maximum grain sizes of 19 μm and 1.8 mm, in which cases only 10% and 1% of the total dust mass is in the small-size 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 Rad-HD simulations, we employed a higher-resolution 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 (). Dashed, white lines indicate the τP(Ts)=1 surface for stellar photons computed using Planck-averaged opacities, solid white lines show the τP = 1 surfaces for vertically traveling infrared photons, and green dash-dotted lines correspond to the z = ±H surfaces. |
The resulting initial temperature distributions for both fdg 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 fdg = 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 fdg = 10−3 and 10−4, respectively, where R is the cylindrical radius. This corresponds to a disk aspect ratio , where h0 = 0.054 and 0.052 for fdg = 10−3 and 10−4, respectively, while H is the pressure scale height measured at the midplane. The flaring index, , 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 fdg = 10−3 (second from the top) and 10−4 (third from the top) in units of the local squared Keplerian frequency . Bottom panel: squared vertical Brunt-Väisälä frequency for both fdg values in units of (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(Ts)=1, where τP(Ts) is the same optical depth computed with the Planck-averaged opacity at the star temperature. Throughout this work, we shall use the terms “upper” or “surface layers” and “middle layer” to refer to the disk regions above and below that surface, respectively. In the same figure it is shown that this temperature transition occurs a few scale heights above the region where the disk becomes optically thin for vertically traveling reprocessed photons, defined as the vertical τP = 1 surfaces computed with (Fig. 1). In our models, the total vertical optical depth computed in this way from z = 0 to infinity is τP ≈ 10 − 20 for fdg = 10−3 and 0.7 − 2 for fdg = 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),
which is roughly proportional to r−1/2 if opacities vary slowly with temperature and τP(Ts)≪1. This explains the obtained midplane temperature power-law index q ∼ 0.5 (steeper profiles would require the inclusion of an additional heat source, for instance, due to viscosity). Under this surface, the temperature is determined by the transport of reprocessed radiation. As described in Melon Fuksman & Klahr (2022), the fact that we are using a fluid approach to describe this process causes a slight deviation from isothermality right below the temperature transition (z/r ∼ 0.1 − 0.2 in Fig. 2), which barely alters the resulting shear distribution with respect to an isothermal stratification (see Fig. 2).
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:
where , with an analogous expression for the epicyclic frequency , quantifies the vertical variation of the specific angular momentum jz = R2Ω, while Ω = vϕ/R is the orbital angular frequency and . This equation, which is a form of the so-called thermal wind equation in geophysical scenarios, shows that a vertically varying rotation velocity occurs if the disk is baroclinic, namely, such that ∇p × ∇ρ ≠ 0. As evidenced by the last line, baroclinic stratifications naturally occur in vertically isothermal disks as long as ∂RT ≠ 0. For this reason, below the irradiation surface, 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 becomes important, and the obtained vertical shear results from the balance of both terms on the right-hand side of Eq. (5), which have opposite signs. In particular, is positive above the midplane and negative below, and the opposite is true for . 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 and are about 5 times larger than the maximum , 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 α ∝ q2 in Manger et al. (2021). For instant cooling, a local Boussinesq linear analysis predicts that the local growth rate of the fastest-growing mode is (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 Brunt-Väisälä frequency corresponding to vertical buoyant oscillations, (e.g., Rüdiger et al. 2002). This quantity is larger than in isothermal disks at the region of maximal shear and smaller above (see Fig. 2), due to the larger and smaller vertical entropy gradient in these locations, respectively. The combined effect of these quantities on the disk stability is studied in 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. Rad-HD simulations’ setup
Starting from the described initial conditions, we ran Rad-HD simulations in order to study the growth and saturation of the VSI. A list of the presented runs and their corresponding parameters is shown in Table 1. Simulations are labeled as, dgNdgcNc_Nθ where the numbers Ndg and Nc indicate the employed fdg and values, respectively, while Nθ is the resolution in the θ-direction. We employed four different resolutions, Nr × 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 . 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 , respectively. In particular, it is shown there that the radiative cooling time is unaffected for high enough , which is crucial for this work given the dependence on the VSI on that timescale.
The Rad-HD equations are solved using third-order Runge-Kutta time integration for the gas fields and the IMEX method derived from it in Melon Fuksman & Mignone (2019) for the radiation fields. Gas and radiation fluxes are computed using the HLLC solvers in Toro (2009) and Melon Fuksman & Mignone (2019), respectively. We employed in every case third-order WENO spatial reconstruction (Yamaleev & Carpenter 2009) with the modifications introduced in Mignone (2014) for spherical grids. In all cases, Dirichlet boundary conditions are imposed by fixing both radiation and HD fields in the ghost cells to their corresponding values in the initial hydrostatic configurations.
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 highest-resolution simulations are shown in Fig. 3 after 50 and 300 orbits (we define one orbit as the rotation period T5.5 at 5.5 au). For fdg = 10−3, as shown in the velocity snapshots at t = 50 T5.5, the instability first produces “finger” modes (Nelson et al. 2013) characterized by an approximately vertically antisymmetric vz, 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 fdg = 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 ∼102 times lower than for fdg = 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 in our highest-resolution simulations after 300 orbits (also shown after 50 orbits for fdg = 10−3). Bottom row: vertical mass flux in the same snapshots, increased by a factor of 100 for fdg = 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 fdg = 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 fdg = 10−3 and 10−4, respectively. At the middle layer, rms velocities are about 3 times lower than at the upper layers for fdg = 10−3, whereas for fdg = 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 volume-averaged meridional kinetic energy, , shown in Fig. 5 for varying dust content and resolution. For fdg = 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 Flores-Rivera et al. (2020). The saturation of the body modes occurs at t ∼ 125 − 160 orbits in all cases except for Nθ = 256, in which case body modes take up to ∼400 orbits to grow due to the higher numerical diffusion. On the contrary, for fdg = 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 volume-averaged kinetic energy (top) and maximum gas velocity (bottom) in all runs with . |
It can be seen in Fig. 5 that the saturation value of ⟨ϵk⟩ is smaller for increasing resolution for fdg = 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 fdg = 10−4 reaches convergence for Nθ ≥ 512 after ∼30 orbits, when the saturation of the VSI modes occurs above the irradiation surface. However, changes in the gas flow can still be seen for increasing resolution, not only in the wavelength of the body modes but also in the increased growth of parasitic Kelvin-Helmholtz (KH) eddies, visible in Fig. 3 as ripples in between the vertical flows and further studied in Paper II. As discussed in that work, the saturation of the VSI may result from the transfer of kinetic energy from the VSI flows into KH eddies. This is consistent with the earlier VSI saturation for increasing resolution, and may also explain the decreasing saturated value of ⟨ϵk⟩ with resolution for fdg = 10−3 if the Kelvin-Helmholtz instability (KHI) limits the maximum energy of the saturated body modes. This effect may be enhanced by the favored growth of the KHI when numerical diffusion is reduced, and also by the increase with resolution of the number of KH-unstable interfaces between vertical flows, given that larger resolutions lead to smaller VSI wavelengths (Figs. A.2 and C.3). On top of this, some KH eddies are subject to baroclinic amplification, as also analyzed in Paper II. The fast vortices produced by this process are responsible for the maximum velocities in 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, jz = R2Ω. In these snapshots it can be seen that the VSI flows generate regions of approximately uniform jz, where consequently the vertical shear is reduced with respect to its initial state, as . As explored in Paper II, these bands begin forming due to the angular momentum redistribution produced by the motion of VSI-destabilized gas parcels across constant-jz surfaces (see also James & Kahn 1970). The subsequent reduction of the vertical shear causes the baroclinic torque to prevail over the centrifugal torque. Once the vertically moving flows are connected in upper disk regions, angular momentum can be transported in the sense of circulation favored by the baroclinic torque, until uniform-jz 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 highest-resolution 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 quasi-Keplerian initial state (not shown here), since vϕ = jz/R ∼ R−1/2 for a quasi-Keplerian flow and ∼R−1 for constant jz. In each case, the shear in all directions is always maximum in between the uniform-jz zones.
Overall, the increase of the shear in between uniform-jz regions, together with the extra shear produced by KH eddies, result in an increase of the time- and volume-averaged absolute value of the vertical shear with respect to the initial state in all regions where the VSI is active, as shown in Fig. 7. For fdg = 10−4, this occurs everywhere except in the VSI-inactive middle region, where the shear remains unchanged. The growth of this quantity increases with resolution, almost by one order of magnitude in our highest-resolution runs. This effect is likely enhanced by the favored growth of small-scale structures such as KHI eddies as numerical diffusion is reduced (Paper II). Conversely, equivalent averages of the signed ∂zvϕ 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 constant-r 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 time-averaged Reynolds stress produced by the VSI. Since our disk model is axisymmetric, we replaced the ϕ-averages normally employed in this calculation with time averages between 150 and 300 orbits (after the saturation of the body modes for Nθ ≥ 512), taking them as representative of an ensemble average over fluctuations. We thus computed the ϕ-row of the Reynolds stress tensor as
where ⟨ ⋅ ⟩t indicates time average and δ denotes fluctuations with respect to the mean. The components of this tensor can be turned into stress-to-pressure ratios which can be used in effective 1D disk evolution models (e.g., Burn et al. 2022). However, as with other forms of disk turbulence (Klahr & Bodenheimer 2003), the VSI-induced angular momentum transport cannot be modeled by means of an isotropic viscosity coefficient (Shakura & Sunyaev 1973), since the vertical stresses produced by the VSI are orders of magnitude larger than the radial stresses, as shown in Stoll et al. (2017) and next in this section. Despite this, the R-component of Sϕ can still be used to compute the flux of vertically integrated mass and angular momentum, and therefore it can still be parameterized as an α-viscosity model implementable in 1D models. This can be seen by taking the azimuthal average of the angular momentum conservation equation as done, for example, in Balbus & Papaloizou (1999), which in axisymmetric disks can be replaced with the mentioned average over fluctuations. This results in the following expression:
where ⟨ ⋅ ⟩ denotes the chosen averaging procedure, and where we assumed ⟨δvϕδρ⟩=0. The second and third terms in this equation account for the transport of angular momentum due to gas advection and VSI-induced mixing, respectively. In quasi steady-state accretion, the projection of the average mass flux perpendicular to the surfaces of constant specific angular momentum is thus determined by the divergence of RSϕ. Since such surfaces are approximately vertical close to the midplane, 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 ⟨ρvz⟩ and Sϕz are null at z = ±∞, this results in
(see, e.g., Jacquet 2013), where Ω is evaluated at the midplane. In turn, the column density evolves according to
and so the radial transport of vertically averaged mass and angular momentum in steady state depends solely on SϕR.
In Fig. 8 we show pressure-weighted radial averages of SϕR/⟨p⟩t and |Sϕz|/⟨p⟩t. In every case |Sϕz|/⟨p⟩t is between 1 and 2 orders of magnitude larger than SϕR/⟨p⟩t, which evidences that the mixing of angular momentum in the disks occurs predominantly in the vertical direction, as obtained in Stoll et al. (2017). In particular, Sϕz is positive for z > 0 and negative for z < 0, which means that angular momentum is transported away from the midplane. In the 2D profiles of ||Sϕ/⟨p⟩t|| in Fig. 9, we note that this quantity is maximum in between the approximately uniform-jz 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 fdg = 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 (i.e., the ϕR-component of the strain rate tensor of a fluid with viscosity νR), while reproducing the usual parameterization for thin, vertically isothermal disks, we define
or equivalently, (see, e.g., Armitage 2022), where ΩK(R) is the Keplerian angular velocity (the extra factor accounts for the fact that in thin, vertically isothermal disks). We have verified that the graphs of αR defined in this way are practically indistinguishable from those of SϕR/⟨p⟩t,
and so the SϕR/⟨p⟩t averages
shown in Fig. 8 can be interpreted as αR. To obtain an effective 1D viscous evolution model similar to, e.g., Lynden-Bell & Pringle (1974), we can then replace in Eq. (8), which leads to
where ⟨αR⟩p is the pressure-weighted z-average of αR (the corresponding expression in Jacquet (2013) has a multiplying factor of 3 instead of 2 due to their different definition of α). In these expressions we assumed a vertically uniform squared isothermal sound speed , which is a valid approximation given the temperature distribution in the disk middle layer.
Values of ⟨αR⟩p computed via interpolation along constant-R lines are shown in Fig. 10. These are in the range ∼10−6 − 10−5 and ∼10−7 − 10−6 for fdg = 10−3 and fdg = 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. Pressure-weighted z-averages of αR computed for different resolutions and fdg 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 fdg = 10−3 and fdg = 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 . 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 = αzH2Ω given that, unlike in the radial direction, turbulent mixing transports angular momentum and dust particles in the same direction (i.e., away from the midplane), following both the negative angular momentum and dust density gradients. The ϕz-component of the strain rate tensor in cylindrical coordinates is of the form , and so in order to parameterize the eddy viscosity as we define
Even though 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 pressure-weighted α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 fdg = 10−4, αz decreases from values up to 10−2 to almost-zero at the midplane, where the instability is suppressed. On the other hand, for fdg = 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 fdg = 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 so-called 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 large-scale outward transport of dust particles (Takeuchi & Lin 2002; Stoll & Kley 2016), which has been proposed (e.g., Ciesla 2007) to contribute to the puzzling observed presence of high-temperature minerals in outer regions of protoplanetary disks (Olofsson et al. 2009) and our own solar system (Brownlee et al. 2006, Bryson & Brennecka 2021, and references therein) that likely formed in hotter, inner regions.
Fig. 13. Absolute value of the average radial velocity at r = 5.5 au in our highest-resolution simulations. Blue and orange lines correspond respectively to outward and inward velocities, respectively. The pressure scale height is indicated for reference (gray dotted line). |
A meridional circulation can also be seen in the dust-depleted 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 VSI-active disk surface layers, as shown, e.g., in Figs. 8 and 9. Therefore, outward flows can still occur in upper layers even if the VSI is suppressed at the disk middle layer. Disk regions with an unstable surface layer and a VSI-inactive middle layer can be expected for sufficient depletion of small dust (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/T0 = (T − T0)/T0 comparing the temperature distributions T and T0, computed respectively at t = 300 and 4 orbits (the latter being after the initial system relaxation but before VSI modes grow), is shown in Fig. 15 for the maximum-resolution runs.
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 fdg = 10−3, some heating occurs at the disk middle layer due to the radiative dissipation of the mechanical energy of the VSI modes. Thermal equilibrium is reached when radiative cooling balances the VSI heating rate. In our highest-resolution run, this occurs at approximately 250 orbits, i.e., about 100 orbits after the saturation of the body modes, as shown, for instance, by the midplane temperature perturbation profiles in Fig. 16. The temperature perturbations are maximum at the maximum-shear layers in between uniform-jz regions (see Fig. 6), with maximum ΔT/T0 ∼ 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 uniform-jz 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 small-scale amplified eddies (Paper II).
For fdg = 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/T0, 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 fdg = 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/T0 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 fdg = 10−3 results from the fact that the disk is both VSI-active in that region and vertically optically thick.
Fig. 17. Radially averaged relative temperature variation for all employed resolutions and fdg values. Displayed values have been averaged in time between 285 and 300 orbits to get rid of fluctuations. The values for fdg = 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 fdg values. |
By Reynolds-averaging the HD equations, the average local heating power density can be computed as
where νR and νz are the effective eddy viscosities computed in the previous section. We have verified that the term proportional to νz in this expression can be safely neglected, as it can be shown taking Ω ≈ ΩK that the ratio between the z- and R-contributions to Q+ is approximately , which is ≪1 despite αz ≫ αR. Resulting r-averages of the heating rate Q+/(ρϵ) in units of Ω are shown in Fig. 18. For fdg = 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 fdg = 10−4, the heating rate is only nonzero at the VSI-active surface layers, which is not enough to heat the midplane.
Due to the vertical dependence of αR, our obtained heating distributions resemble those in the simulations of MRI-heated disks by Flock et al. (2013), which are also locally minimum at the midplane and peak a few scale heights above it. This observation contradicts the usual assumption resulting from α-viscosity models that turbulent heating is maximum at the midplane, which is commonly used to compute temperature distributions via 1D radiative transfer (e.g., Nakamoto & Nakagawa 1994; Pfeil & Klahr 2019). The vertical temperature distributions resulting from this method peak at the midplane (e.g., Mori et al. 2019), as opposed to our obtained distributions (Fig. 17).
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 Richardson-like criterion in which the VSI is subject to the instability condition
where trel is the local thermal relaxation timescale, while Nz is the vertical Brunt-Väisälä frequency (see Fig. 2). A linear analysis in Klahr (2023) for arbitrary cooling times similar to that in Goldreich & Schubert (1967) shows that the linear growth rate decreases as ∼1/trel for trel < tcrit, where the critical cooling time tcrit is defined as the right-hand side of Eq. (14). On the other hand, the vertically global instability criterion by Lin & Youdin (2015) reads
where H and Γ are respectively the gas pressure scale height and heat capacity ratio. This expression, which corresponds to 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 ∝ Rq, and (II) that β = trel/Ω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 tcrit defined as the right-hand 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 tcrit (top) and the orbital frequency (bottom) at r = 5.5 au for both fdg values. Also shown are the corresponding values for t = tcrit and the locations z = ±ΓH/2 (gray dotted lines) at which the local tcrit coincides with its value in the global stability criterion (Eq. (15)). |
In the simulations presented in this work, the thermal relaxation time trel is entirely determined by radiative transfer, and equals1 the local radiative cooling time, tcool. We estimate this timescale in Appendix B via a local lineal analysis assuming small temperature perturbations with a given wavelength λ on a uniform background, obtaining
with
where |k|=2π/λ, and are respectively the Planck- and Rosseland-averaged mean free paths, and . This expression evidences a dependence of the cooling mechanism on λ: temperature perturbations occurring over much smaller lengthscales than λR and λP (λPλRk2 ≫ 1) get damped at the scale-independent rate corresponding to thermal emission in an optically thin medium. On the other hand, optically thick perturbations (λPλRk2 ≪ 1) cool down through radiative diffusion at a rate , which increases proportionally to k2. The adiabatic limit tcool → ∞ is thus recovered for λRk → 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 first-order perturbation analysis is applicable, but different cooling time prescriptions used in the literature (e.g., Stoll & Kley 2016; Malygin et al. 2017; Dullemond et al. 2022) only differ slightly and should yield equivalent results. In the same appendix, we show that the resulting cooling times are not affected by the RSLA for large enough , in particular for our employed value of .
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λRk2 ≫ 1 everywhere, i.e., induced temperature perturbations of the same wavelength as the velocity flows are optically thin. This is true at all times, and since the optically thin cooling timescale is wavelength-independent, the estimated tcool values are unchanged through time. This validates the employment of β-cooling prescriptions in even higher-resolution axisymmetric simulations or 3D simulations using the same disk models, with the difference that such an approach cannot reproduce the (rather negligible) disk’s thermal evolution (Sect. 5). Further differences can be expected if large-scale optically thick structures are eventually formed (e.g., vortices in 3D).
Resulting slices of tcool/tcrit and tcoolΩ at r = 5.5 au are shown in Fig. 19 for each fdg. 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, tcool decreases up to 2 orders of magnitude despite the surface layers having much smaller optical depths. The reason for this is that the ratio ρϵ/aRT4, 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 fdg = 10−3, the resulting tcool is everywhere below tcrit, which is consistent with the obtained instability in the whole domain.
For fdg = 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 tcool below the irradiation surface forms two layers above the midplane where tcool > tcrit, 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 , 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 (vz ≪ cs) vertical modes with similar amplitude in the midplane region and in the inactive zones above and below. One possible explanation for this is the fact that the local linear VSI growth rate equals the local vertical shear (e.g., Urpin & Brandenburg 1998; Klahr 2023), which is very small close to the midplane. However, under this criterion, and provided the numerical dissipation in our code is low enough, we should still see some growth occurring after several hundreds of orbits, as it can be seen in the local growth rate distribution shown in Fig. 1. Another explanation can be given if the maximum amplitude of the saturated VSI perturbations is limited by the local vertical shear. A hint in this direction is given by the vertically global model by Lin & Youdin (2015), which shows that the VSI growth rates are reduced if the maximum height of the domain, and therefore also the z-averaged vertical shear, are decreased. This hypothesis is also supported by the approximate relation ⟨αr⟩ρ ∝ q2 in the global vertically isothermal model by Manger et al. (2021), which shows that the radial Reynolds stress in the saturated stage increases for larger average vertical shear. Lastly, resolution might also play a role in limiting VSI growth close to the midplane. The linearly unstable directions, strictly confined between the constant-R and constant-jz surfaces (e.g., Knobloch & Spruit 1982; Klahr et al. 2023), are for small |z| only slightly slanted with respect to the vertical direction, and thus large resolutions might be needed to reproduce linear growth in such directions. This can be investigated in high-resolution box simulations of the midplane region.
However useful, this picture has the limitation that VSI modes are nonlocal. Thus, another possibility is that the amplitude of the flows at the middle layer does not only depend on the in situ driving of the VSI but also on the excitation of inertial modes in unstable layers extending to otherwise inactive regions, as it likely occurs in the T states in Fukuhara et al. (2023). This, together with the fact that VSI modes can in general occupy regions of largely different cooling times, makes it difficult to apply local stability criteria in more complex tcool 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 jz is likely enhanced by the imposed axisymmetry of our models. In 3D simulations, we would expect constant-jz zones to be unstable to non-axisymmetric modes (e.g., Papaloizou & Pringle 1984), in which case it remains to be seen if the VSI can still enforce the formation of reduced shear bands. Surprisingly, this is likely the case in a number of published works: in the 3D simulations by Richard et al. (2016), the vertical component of the vorticity, ωz, proportional to ∂Rjz for axisymmetric flows, peaks in between VSI modes, indicating abrupt jumps in jz. 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 jz. More evidence showing abrupt jumps in jz is given by the simulations by Manger & Klahr (2018), in which the azimuthally averaged critical function ℒ, approximately proportional to 1/∂Rjz, has peaks of similar width as the VSI modes. On top of this, the 3D Rad-HD simulations of irradiated disks by Flock et al. (2020) exhibit a formation of approximately uniform-jz bands which is most predominant in the maximum-shear upper layers, as shown in Fig. A.1. It is then possible that non-axisymmetric modes, while preventing the formation of uniform-jz zones, may still permit the creation of regions with reduced ∇jz (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 dust-to-gas ratio fdg = 10−3, where the instability is not suppressed close to the midplane, are slightly lower than the values αr ≈ 10−4 reported in Flock et al. (2020) with H/R = 0.09 − 0.14 and Manger et al. (2021) and Barraza-Alfaro et al. (2021) with H/R = 0.1. On the other hand, Pfeil & Klahr (2021) and Manger et al. (2020) obtained αr ∼ 10−5 respectively in 2D and 3D simulations with H/R ∼ 0.05. Considering the dependence of αr with H/R studied in Manger et al. (2020, 2021), who propose αr ∝ q2(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 high-resolution simulations in Flock et al. (2020), Barraza-Alfaro et al. (2021; see, however, Richard et al. 2016; Manger & Klahr 2018), and thus it is likely that the mixing of angular momentum in 3D is predominantly axisymmetric. We defer such verification to follow-up works on this topic. Still, differences in the spatial distribution of averaged velocities and stresses in 2D and 3D simulations are seen in Stoll & Kley (2014, 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 r-flux along constant-r surfaces, we obtain that the positive and negative contributions cancel each other out almost exactly. More precisely, the θ-averaged radial mass flux takes small positive values between 2 and 3 orders of magnitude smaller than the maximum values shown in Fig. 14. This occurs in contradiction with 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 inward-directed gas flow is forced at the inner radial boundary to turn away from the midplane and merge with the outward-directed meridional flow, while the inverse process occurs at the outer radial boundary. This increases the outward flux, creating an equilibrium configuration in which outgoing and ingoing flows balance each other while no mass escapes the domain. The extent to which this affects the resulting gas structure in our simulations and all mentioned works is unknown. This drawback could be solved by including the inner rim in the domain and allowing gas to escape the radial boundary toward the star, ideally including magnetic effects in well-ionized regions (Flock et al. 2017a).
The ideal MHD simulations of MRI turbulence by Fromang et al. (2011) and Flock et al. (2011) show similar downward-concave parabola patterns as in 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 VSI-unstable 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 smaller-wavelength modes can grow. Given our obtained trend of decreasing VSI strength with decreasing dominant wavelength (see Figs. 8 and 18), we can expect in that region a suppression of the VSI at the midplane due to the reduction of the maximum VSI-unstable wavelengths. It therefore appears unlikely that the VSI can in general lead to any significant heating, since that would simultaneously require (I) large enough vertical optical depths to retain the produced heat and (II) modes with small radial wavelengths producing larger levels of turbulence than predicted by simulations. This hypothesis can be further studied in high-resolution Rad-HD simulations of the inner disk inside of ∼1 au, keeping in mind that close enough to the star the MRI is expected to become predominant and suppress the VSI (Lesur et al. 2023).
The model implemented in Pfeil & Klahr (2019) also predicts that larger αR values of ∼10−4 − 10−3 can lead to an increase of up to a few tens of K outside of 1 au. Such high αR values can only be expected for steeper temperature profiles than in this work (T ∝ R−1, Pfeil & Klahr 2021; Manger et al. 2021) or for large aspect ratios at a few au of up to 0.1 (Manger et al. 2020), both of which are hard to justify in passive irradiated disks. In this context, it is also unclear whether steeper temperature gradients, predicted by models of viscous-dominated heating in optically thick accretion disks (, Pringle 1981), can be maintained solely by VSI-induced turbulence. Moreover, the cited α-disk model in Pfeil & Klahr (2019) assumes maximum heating at the midplane, whereas our maximum heating regions occur in upper regions of lower optical depth, and thus the temperature increase in that work can be regarded as an upper limit in the case of VSI-induced heating. Further explorations are however required to test whether any heating can be produced by the VSI under a plausible choice of parameters.
The fact that the temperature distributions in our nominal disk are at all times minimum close to the r-boundaries (Fig. 16) evidences that the obtained temperature increase is not only limited by vertical diffusion but also by radiative cooling through the radial boundaries. This is a numerical effect produced by the employed Dirichlet boundary conditions, which leads to outflow of radiation as soon as the temperature inside of the domain exceeds that at the boundary. Larger domains could be used in the future to explore whether this effect limits the maximum temperature perturbations.
Determining the equilibrium temperatures produced by the VSI is also made difficult by the lack of convergence of the temperature distributions with resolution due to the reduction of the VSI strength, as shown in Fig. 17. However, the average temperatures obtained with our two highest resolutions are rather similar, with differences of ∼10%. Whether such resolutions are enough to accurately predict the temperature increase produced by the VSI can be explored in the future via even higher-resolution simulations. However, given our obtained trend of decreasing heating with increasing resolution, this effect should not alter the conclusion that the heating produced by the VSI is most likely negligible.
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 tcrit 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 . This is in fact verified in our disk models, since tcool is smaller and larger than for fdg = 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 ΔLs/2 of the stability region predicted by the local application of Eq. (15) is at least ∼H. Given the functional form of tcool(z) in that work, ΔLs/2 must be larger than z = ΓH/2 to verify tcool ≫ tcrit at z = ±ΓH/2. Thus, the midplane stability for ΔLs/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 tcool 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 dust-gas thermal equilibration time (see, e.g., Barranco et al. 2018). This becomes important at the disk upper layers, where the dust and gas densities are low enough that trel is increased due to their infrequent collisions, with a stabilizing effect in those regions (Pfeil & Klahr 2021). On the other hand, trel 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 two-moment radiation-hydrodynamical simulations with realistic temperature stratifications and cooling times, both of which result from the balance of stellar irradiation and radiative cooling, which are treated in our code using realistic tabulated dust opacities. We assumed typical parameters for disks around T Tauri stars considering two different values of the dust-to-gas mass ratio of small grains, namely a nominal case with fdg = 10−3 and a dust-depleted case with fdg = 10−4, making our studied disk region between 4 and 7 au vertically optically thick and thin, respectively. Our main findings can be summarized as follows:
-
1.
The temperature gradient in our disk models leads to a localized increase of the vertical shear in the temperature transition region between the midplane and the atmosphere. The radial and vertical components of the temperature gradient tend to increase and decrease the overall vertical shear, respectively. At its peak, the resulting vertical shear reaches up to 2 times its value in equivalent vertically isothermal models, which is translated into similarly larger local VSI growth rates than predicted by linear theory for vertically isothermal disks. The stabilizing buoyancy frequency is also maximum in that region, while above it decreases below its vertically isothermal value.
-
2.
In our nominal case, we observe two distinct growth stages corresponding to the growth and saturation of finger and body modes occupying the entire domain. Instead, in the dust-depleted case, the VSI is suppressed at the disk middle layer, and body modes never grow. This behavior results in a single growth phase of VSI modes localized at the disk surface layers.
-
3.
In every case, the VSI reduces the vertical shear inside of bands of approximately uniform specific angular momentum. Conversely, shear is maximal at the interfaces between these bands. Even though this process is likely enhanced by the enforced axisymmetry, bands of reduced shear can also be seen in 3D simulations, suggesting that their formation mechanism may in some cases predominate over their destruction by non-axisymmetric modes.
-
4.
The time- and volume-averaged absolute value of the vertical shear increases with resolution as a result of the growth of small-scale velocity fluctuations, exceeding the mean shear value by almost one order of magnitude in our highest-resolution runs. This is not translated into an increase of the VSI strength, which is instead enforced by the disk’s baroclinicity and limited by the VSI saturation mechanism.
-
5.
The Reynolds stress is in every case maximum at the maximum shear regions above the irradiation surfaces and in between the approximately uniform-jz bands, with z components between 1 and 2 orders of magnitude larger than the R components. The pressure-weighted αr values of ∼10−5 obtained in the nominal case coincide with reported values for 3D vertically isothermal simulations. This suggests that, even in 3D, the radial mixing of angular momentum may be predominantly axisymmetric.
-
6.
The anisotropy of the Reynolds stress tensor causes differences between its radial components in cylindrical and spherical coordinates, with αR < αr by factors of ∼4 and ∼2 in the nominal and dust-depleted cases, respectively. Determining whether this is a general trend for different disk parameters requires further testing. While αR should be used to compute the vertically integrated mass flux, αr is usually reported in the literature.
-
7.
As in previous works (e.g., Stoll & Kley 2016), the anisotropic angular momentum transport produced by the VSI produces a meridional circulation with inflow at the midplane and outflow in upper layers. This pattern only occurs at the disk’s surface layer in the dust-depleted case. This circulation does not lead to accretion in our simulations due to our imposed no-outflow and no-inflow boundary conditions, which likely increases the resulting outflow velocities with respect to their values if mass could escape the domain.
-
8.
The turbulent heating produced by the VSI in our models causes negligible temperature increases of 0.1% and 0.01% in the nominal and dust-depleted models, respectively. In the nominal case, the energy dissipated by the VSI modes heats the disk middle layer with a heating rate distribution peaking a few (∼3) scale heights above the midplane, contrarily to the heating distributions in α-viscosity models, which peak at the midplane. We conclude that the VSI is in general unlikely to produce any significant temperature increase, as that would either require it to produce larger heating rates than expected in dense disk regions inside of ∼1 au, where largely optically thick radial wavelengths cease to be VSI-unstable, or unrealistically large radial temperature gradients or aspect ratios leading to larger average αR values (10−4 − 10−3) than can be obtained in models of passive disks. Further explorations of the parameter space via high-resolution Rad-HD simulations of inner disk regions should either deny or confirm this hypothesis.
-
9.
The Richardson-like local criterion by Urpin (2003) explains the instability regions observed in our simulations, showing that, despite VSI modes are vertically global, a local analysis can provide good insight into the stability of specific disk regions. On the other hand, the global criterion in Lin & Youdin (2015) derived for vertically isothermal disks with uniform cooling time can still be applied to the vertically isothermal disk middle layer, which justifies its general application to assess the stability close to the midplane as long as VSI modes are radially optically thin.
-
10.
We verified via analytical computations and numerical testing that the RSLA is safely applicable to this problem, since for high enough 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 tcool if this value is similar or larger than tcrit. This also shows that IMEX schemes for Rad-HD, with generally good scalability properties in comparison with globally implicit methods, can be applied to this problem without introducing unphysical 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 low-density disk surface layers, in particular in outer disk regions currently observable with ALMA (e.g., Andrews et al. 2018). These effects are included in Paper II to study the stability of such regions.
In reality trel = Γtcool if one includes the effect of PdV work in the relaxation process (see, e.g., Goldreich & Schubert 1967; Lesur et al. 2023; Klahr 2023), but this widely used approximation is accurate enough for the purpose of a comparison with tcrit, especially considering that the suppression of the VSI for trel > tcrit as a function of trel 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/16-1/2. We thank our collaboration partners on this project in Kiel, Sebastian Wolf and Anton Krieger, under contract WO 857/17-1/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 Max-Planck Computing and Data Facility in Garching (Germany). We thank the anonymous referee, whose comments helped us improve the quality of this work.
References
- Al-Mohy, 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 e-prints [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]
- Barraza-Alfaro, 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]
- Flores-Rivera, 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., Murray-Clay, 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]
- Lynden-Bell, 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 Rad-HD disk simulation by Flock et al. (2020) after 350 inner orbits of evolution. Approximately uniform-jz 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 distribution taken from the initial condition of each Rad-HD 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, Er, Fr)=(p0 + δp, Er0 + δEr, Fr0 + δFr), where (p0, Er0, Fr0) are the hydrostatic unperturbed values and (δp, δEr, δFr) are time-dependent perturbations. Inserting this expression into Equation (1) and keeping only linear terms in the perturbations, we obtain the following system of equations:
where the quantity η = 4 aRT4/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 , while in this work we neglect scattering and assume . We then simplify this expression by considering one-dimensional radiation transport along a single direction n0 = Fr0/||Fr0||. In this way the radiation pressure tensor takes the simple form , 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 , 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:
where the matrix A is defined as
while and are respectively the Planck and Rosseland mean free paths and , which is always positive for f > 0 and vanishes for f = 0.
Given an initial temperature perturbation, which in this analysis takes the form (δp0, 0, 0), the time-dependent solution of Eq. (B.2) can be obtained via matrix exponentiation as eAt(δp0, 0, 0)⊺. The matrix A has three different eigenvalues, {sc, s+, s−}, and therefore all perturbations evolve as a linear combination of {esct, es+t, es−t}. We derive expressions for these eigenvalues in Sections B.1–B.3 and show in Section B.4 that δp ∼ esct, where |sc|≪|s±|. The cooling time is therefore given by , while the other two eigenvalues correspond to the faster adjustment of the radiation fields. These conclusions remain unaltered when for high enough . We study the applicability of the RSLA in Section B.5, where we use this analysis to obtain constraints to the values of 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
where I3 is the rank-3 identity matrix and
In principle we can obtain all three roots of P by applying a cubic formula, but instead we can get accurate analytic expressions for these roots by considering the limits of small and large |s|.
We first consider the limit
which also guarantees , since our employed opacities verify for all temperatures (Melon Fuksman & Klahr 2022; this is also true in general for opacity laws of the form as long as n > 0.81). In that case, the characteristic equation P(s)=0 loses its dependence on and becomes
whose solution is
This expression can be simplified by noting that the term iξ′λRk is much smaller than 1 for large λ (λRk ≲ 1) and much smaller than the 𝒪(λRλPk2) term for small λ (λRk ≫ 1), and thus it can be neglected. On the other hand, the factor ξ − ξ′f only affects the value of sc if λRλPk2 is not much larger than 1, which can only occur if f ≪ 1, in which case ξ′f ≈ f2 ≪ 1 and ξ ≈ 1/3. Therefore we can safely replace ξ − ξ′f by 1/3, which leads to the simpler expression
Using the k distributions in our simulations (Section 6), we estimated the relative error of this approximation to be of at most 0.01%.
If this analysis is repeated in 3D, the same expressions are obtained with the difference that λPλRk2 is replaced by a term of the same order of magnitude that tends to λPλRk2/3 wherever λPλRk2 ≲ 1. Thus, the particular form of the 𝒪(λPλRk2) terms for anisotropic transport does not affect the resulting cooling rate, which always tends to sthin for λPλRk2 ≫ 1 and to Equation (B.9) otherwise. Consequently, Equation (B.9) still holds when considering multidimensional radiation transport. Analogously, we would have obtained the same result using any other closure than M1 or even solving the full frequency-integrated radiative transfer equation, since always tends to Erδij/3 for λPλRk2 ≲ 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 tcool = |sc|−1. Specifically,
which is the constant-opacity limit of Equation (16). The wavelength-dependent behavior of tcool can be understood by noting that the value of k regulates how much radiation energy is transported away through advection before it can reheat the dust-gas mixture via absorption. This can be seen by comparing the relative size of the different terms in Eq. (B.2). For large k, Er quickly evolves into a state such that δEr ≈ ηδp/(1 + 𝒪(λPλRk2)) ≪ ηδp, which can be seen by zeroing the time derivatives of the radiation fields. This means that δEr can be dropped from the cooling source term of δp, which is proportional to (ηδp − δEr). In this way, radiation transport limits the maximum value of δEr enough to cause temperature perturbations to be damped solely due to optically thin cooling, which is scale-independent. This is not verified for smaller k values if λPλRk2/3 is not much larger than 1, in which case . 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
where and . This expression only depends on bP far from the midplane, where Er ≠ aRT4. On the other hand, the rightmost term proportional to bR only appears in the diffusion regime, in which case it introduces an oscillatory component. Since in that case , it is straightforward to prove that the relative variation of sc due to that term is at most 4bRλ/L, where L is the variation lengthscale of Er. Given the initial assumption that λ ≪ L, that term can be neglected in all regimes and the resulting cooling time is given by Equation (16).
B.2. Consistency check and constraints on
We now test the consistency of the small |sc| condition assumed in Eq. (B.6). This limit can be better understood by considering the components of δEr and that evolve proportionally to esct + ikx, which we denote as and . Those components must also satisfy Eq. (B.2), which leads to
meaning that if sc satisfies Eq. (B.6), the time derivatives of the radiation fields can be dropped. This means that, in that case, and evolve quasi-statically, as their readjustment time is much shorter than |sc|−1, which explains why disappears from the equations. In particular, if f ≪ 1, this leads to the usual expression for the radiative flux in the diffusion limit, , whose perturbation in this analysis becomes .
Using the obtained expression for sc (Eq. (B.9)), we can see that Eq. (B.6) is satisfied for all k as long as , which leads to the equivalent condition
For , 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 sc is still valid in that region, as we explain next. Far enough from the midplane, the timescale |sthin|−1 becomes shorter than the light-crossing time of a distance λR, typically of a few hundred au in those regions. In such locations, λRλPk2 → ∞ (otherwise λ ∼ λR/P, which would render the uniform background assumption inapplicable), in which case s = sthin provided a large-k condition is satisfied independently from Eq. (B.13). To see this, we rearrange the characteristic equation as
with
from which we can see that s = sthin is obtained as long as |φ(sthin)| ≪ 1. We can then combine both conditions, which leads to the conclusion that Eq. (B.9) is valid provided
that is, as long as either the cooling timescale is much shorter than the light-crossing time of a mean free path or the gas-dust mixture is largely optically thin along a perturbation wavelength. For the wavelengths in this work (see Section 6), this relation holds in the entire domain for all considered values in Section B.5, except for in the case fdg = 10−3, in which the left-hand side of Eq. (B.16) is of order 1 close to the midplane. We further study the effect of reducing on the cooling timescale in Section B.5.
B.3. Remaining eigenvalues
Having verified the validity of Eq. (B.9), it only remains to show that δp ∼ δp0esct, 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|≫|sc|, which transforms the characteristic equation into
Since we only want to obtain the order of magnitude of the two solutions to this equation, s±, it is enough to compute their values in the optically thick and thin limits. If λRk ≪ 1, which also guarantees λPk ≪ 1, the two solutions of these equations are
while for λPk ≫ 1, and therefore also λRk ≫ 1, we have
where the argument of the square root is always positive if f < 1. Therefore, |s±| is at least of order . We can then check the consistency of the assumption |s±|≫|sc|, for which it is enough to verify that |s±|≫|sthin|, since |sthin|≥|sc|. For λRk ≪ 1, this is equivalent to Eq. (B.13), which as already mentioned is verified in that regime for large enough . The same is true for λRk ≫ 1, in which case the resulting condition is , which in that limit holds everywhere for .
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 {esct, es+t, es−t}. The leading orders of δp = (eAt(δp0, 0, 0)⊺)1 can be evaluated in a rather straightforward way by applying the Cayley-Hamilton theorem, which allows us to compute eAt as a linear combination of {I, A, A2} with time-dependent coefficients, where I is the 3 × 3 identity matrix. This results in expressions of the form
if λRk ≲ 1, and
if λPk ≫ 1 , where 𝒪i(x) denotes a term of order x. For , the leading term in these equations is proportional to esct, whereas the other terms are much smaller, both close to the midplane, where η ≪ 1, as well as away from the midplane, where λPk is always much larger than 1. This proves that . On the other hand, the relative size of these terms changes for , which together with the changes in the eigenvalues can affect the overall relaxation rate for sufficiently small . We study this effect in the next section.
B.5. Dependence on the reduced speed of light
In the previous sections we showed that tcool can be different from its physical value when , either due to changes in the value of sc 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 by numerically computing δp(t) in our domain, defining as the e-folding time of δp, i.e., . This definition coincides with the real cooling time when . We achieve this by numerically computing eAt for different times in the interval [0, 2 tcool(c)] using the SciPy (Virtanen et al. 2020) expm matrix exponentiation routine, which follows the method in Al-Mohy & Higham (2010). The A matrix is in turn computed at each location using the k distributions obtained in Appendix D for Nθ = 512.
We repeated this computation for different values in the interval [10−5, 100]. The obtained distribution of the relative error , with , is shown in Fig. B.1 for all considered and fdg values. For the two highest values, we also show in Fig. B.2 the signed relative variation , which shows that the cooling timescale does not always decrease with , since the e-folding time depends on the contribution of the different terms in Eqs. (B.20) and (B.21). The computed relative errors grow for decreasing and stay on the order of 1% and below for all parameters except for , 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 . 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 as a function of for fdg = 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 for and both considered fdg values. |
An analytical constraint for the minimum value of that does not alter tcool can be obtained by noting that, in the midplane region where the largest errors occur for , the most stringent condition for the value of among those derived in Sections B.2-B.4 is given by Eq. (B.16), which guarantees that sc is not affected by the reduction of . This condition is verified except for , in which case the left-hand 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 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 tcool and Ω−1 (or any other dynamical time) to guarantee the independence of the growth of VSI modes on , since if tcool > tcrit, the linear growth rate and the strength of the saturated instability depend on the precise value of tcool (Lin & Youdin 2015; Manger & Klahr 2018; Klahr 2023). It is therefore necessary that tcool is unchanged by the reduction of , which is approximately the case if , as already shown. It remains to be evaluated if is much larger than the maximum gas velocity, which in our runs is on average at most 7% of the isothermal sound speed. Since 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 values at a resolution of 480 × 512. Solid and dashed lines correspond to fdg = 10−3 and 10−4, respectively. |
Fig. C.2. Same as Fig. 8 but comparing the Reynolds stress components for different 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 fdg values, measured at z/r = 0.1 and 0.28 (top and bottom panels, respectively). |
To validate the nominal choice of 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 values spanning from 10−5 to 10−3, with an additional run with for fdg = 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 fdg = 10−4, all curves match almost perfectly except for the random oscillations in |v|max, as expected from the analysis in Appendix B. For fdg = 10−3, we do notice that, both during the growth of the finger modes and the body modes, the runs with 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 values overlap except during the growth of the body modes, where, however, no identifiable trend as a function of is followed. A lower kinetic energy with 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 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 is high enough not to introduce unphysical effects in any stage of the VSI. On top of this, it even seems that 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 VSI-suppressed region for fdg = 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 T5.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 with .
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 ν = αcsH. 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λRk2 for all resolutions and fdg values, indicating the optical depth of the VSI modes (λPλRk2 ≫ 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 cut-off imposed on the VSI modes by the employment of a finite-sized 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 third-order 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λPk2 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λPk2 values are 𝒪(10), and therefore
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 (). Dashed, white lines indicate the τP(Ts)=1 surface for stellar photons computed using Planck-averaged opacities, solid white lines show the τP = 1 surfaces for vertically traveling infrared photons, and green dash-dotted lines correspond to the z = ±H surfaces. |
|
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 fdg = 10−3 (second from the top) and 10−4 (third from the top) in units of the local squared Keplerian frequency . Bottom panel: squared vertical Brunt-Väisälä frequency for both fdg values in units of (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 in our highest-resolution simulations after 300 orbits (also shown after 50 orbits for fdg = 10−3). Bottom row: vertical mass flux in the same snapshots, increased by a factor of 100 for fdg = 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 volume-averaged kinetic energy (top) and maximum gas velocity (bottom) in all runs with . |
|
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 highest-resolution runs. Right panels: corresponding specific angular momentum distributions after 300 orbits. |
|
In the text |
Fig. 7. Same as Fig. 4, but showing instead constant-r 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. Pressure-weighted z-averages of αR computed for different resolutions and fdg 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 fdg. |
|
In the text |
Fig. 13. Absolute value of the average radial velocity at r = 5.5 au in our highest-resolution simulations. Blue and orange lines correspond respectively to outward and inward velocities, respectively. The pressure scale height is indicated for reference (gray dotted line). |
|
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 fdg values. Displayed values have been averaged in time between 285 and 300 orbits to get rid of fluctuations. The values for fdg = 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 fdg values. |
|
In the text |
Fig. 19. Radiative cooling times normalized by the critical cooling time tcrit (top) and the orbital frequency (bottom) at r = 5.5 au for both fdg values. Also shown are the corresponding values for t = tcrit and the locations z = ±ΓH/2 (gray dotted lines) at which the local tcrit coincides with its value in the global stability criterion (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 as a function of for fdg = 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 for and both considered fdg values. |
|
In the text |
Fig. C.1. Average kinetic energy (top) and maximum velocity (bottom) for different values at a resolution of 480 × 512. Solid and dashed lines correspond to fdg = 10−3 and 10−4, respectively. |
|
In the text |
Fig. C.2. Same as Fig. 8 but comparing the Reynolds stress components for different 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 fdg 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λRk2 for all resolutions and fdg values, indicating the optical depth of the VSI modes (λPλRk2 ≫ 1 indicates optically thin cooling). |
|
In the text |
Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.