EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 687 (July 2024) A&A, 687 (2024) L5 Full HTML
Open Access
Issue
A&A
Volume 687, July 2024
Article Number L5
Number of page(s) 6
Section Letters to the Editor
DOI https://doi.org/10.1051/0004-6361/202449323
Published online 26 June 2024

© The Authors 2024

Licence Creative CommonsOpen 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

Thermal relaxation of temperature perturbations in the gas of protoplanetary disks occurs in a two-step process, where collisions thermally accommodate the slowly cooling gas with the quickly cooling dust grains, which then emit or absorb electromagnetic radiation (Malygin et al. 2014, 2017; Woitke 2015; Barranco et al. 2018). The local dust size distribution thus plays a crucial role when it comes to the determination of the thermal relaxation time of the gas. Many thermal and hydrodynamic instabilities depend on this timescale (e.g., convective overstability, Klahr & Hubbard 2014; vertical shear instability, Urpin & Brandenburg 1998; Nelson et al. 2013; zombie vortex instability, Marcus et al. 2015). Turbulence resulting from these instabilities can redistribute the dust grains due to aerodynamic drag, altering the cooling timescale itself. Continuous remodeling of the cooling times in hydrodynamic simulations of protoplanetary disks is thus desirable.

In this Letter, we present axisymmetric simulations of protoplanetary disks that, for the first time, combine the effects of dust grain growth and the emergence of hydrodynamic turbulence on the cooling times. We are specifically interested in the effect on the vertical shear instability (VSI). Various studies have shown its strong dependence on the thermal relaxation time (e.g., Urpin 2003; Lin & Youdin 2015; Manger et al. 2021; Pfeil & Klahr 2021). We note that VSI requires rapid cooling on a timescale ideally shorter than

t c ≲ H g R | β T | γ 2 ( γ − 1 ) | z H g | − 1 Ω K − 1 , $$ \begin{aligned} t_\mathrm{c} \lesssim \frac{H_\mathrm{g} }{R}\frac{|\beta _T|\gamma }{2(\gamma -1)}\left|\frac{z}{H_\mathrm{g} }\right|^{-1}\Omega _\mathrm{K} ^{-1}, \end{aligned} $$(1)

where Hg is the disk’s gas pressure scale height, R is the distance from the central star, βT is the exponent of the radial temperature structure T ∝ RβT, γ is the heat capacity ratio, z is the distance from the disk midplane, and ΩK is the Keplerian angular frequency. Urpin (2003) already showed that the VSI’s growth rate decreases at thermal relaxation times beyond this limit, which has been recently demonstrated in numerical simulations by Klahr et al. (2023). At the same time, VSI creates mostly meridional gas flows, and its ability to vertically corrugate the dust layer has been demonstrated in various simulations (Stoll & Kley 2016; Lin 2019; Flock et al. 2020; Schäfer et al. 2020; Schäfer & Johansen 2022; Fukuhara et al. 2023; Pfeil et al. 2023).

It is, however, unclear whether the VSI evolves fast enough to avoid the effects of dust settling and maintains a dust layer that is thick enough to provide the necessary fast thermal relaxation times. If the grains are large and the cooling times are long, it would be possible that the VSI is not able to develop in the first place or is not able to counteract the settling of the grains. Work by Fukuhara et al. (2023) and Fukuhara & Okuzumi (2024) has demonstrated that an equilibrium state might be possible in which the turbulent mixing of the VSI exactly counteracts sedimentation. Their studies are based on a semi-analytic model of the VSI turbulence and assume constant settling-mixing equilibrium. As radiative cooling, dust dynamics, and dust coagulation are interdependent processes, they would have to be accounted for in a single simulation. Constructing such a fully self-consistent numerical model of the VSI in protoplanetary disks is beyond the scope of this work. If, however, reasonable approximations are made, aspects of this interplay can be studied with much simpler techniques. Here, we present hydrodynamic simulations of protoplanetary disks including four dust fluids of evolving grain sizes in order to investigate the impact of the VSI-induced turbulence on the cooling times. We omitted a full treatment of dust coagulation and only evolved the grain size based on an analytic dust growth model (similar to Birnstiel et al. 2012). Instead of including radiative transfer calculations to model the impact of the evolving dust distribution on the radiative cooling, we relaxed the temperature perturbations on a timescale calculated from the local dust size distribution. Although these methods are only approximations, they enable the first simulations of protoplanetary disks in which the timescale for the thermal accommodation of gas and dust is directly linked to the dynamics of the simulated dust fluids.

2. Methods

2.1. Thermal relaxation times

In protoplanetary disks, thermal relaxation of temperature perturbations in the gas is mostly achieved via thermal accommodation with the dust. We thus approximated

t thin NLTE ≈ t coll gas = γ γ − 1 1 n s v ¯ gas σ s , $$ \begin{aligned} t_\mathrm{thin} ^{\text{NLTE}}\approx t_\mathrm{coll} ^\mathrm{gas} =\frac{\gamma }{\gamma -1}\frac{1}{n_\mathrm{s} \bar{v}_\mathrm{gas} \sigma _\mathrm{s} }\, , \end{aligned} $$(2)

where t coll gas $ t_\mathrm{coll}^\mathrm{gas} $ is the thermal accommodation timescale of the gas with the dust due to collisions (Probstein 1969; Burke & Hollenbach 1983; Barranco et al. 2018). The thermal accommodation timescale therefore depends on the mean molecular gas velocity v ¯ gas $ \bar{v}_{\mathrm{gas}} $ and the Sauter mean radius of the dust size distribution (Sauter 1926):

a s = ∫ n ( a ) a 3 d a ∫ n ( a ) a 2 d a , $$ \begin{aligned} a_\mathrm{s} =\frac{\int n(a)a^3\, \mathrm{d} a}{\int n(a)a^2\, \mathrm{d} a}\, , \end{aligned} $$(3)

where a refers to the grain size and n(a) is the number density distribution of the particles. From this, the number density n s = ρ d / ( 4 3 π ρ m a s 3 ) $ n_{\mathrm{s}}=\rho_{\mathrm{d}}/(\nicefrac{4}{3}\pi \rho_{\mathrm{m}} a_{\mathrm{s}}^3) $ and the collisional cross-section σ s =π a s 2 $ \sigma_\mathrm{s}=\pi a_\mathrm{s}^2 $ follow. Here, ρm is the material density of the dust and ρd is the dust volume density. In most cases, the size distribution of dust in a protoplanetary disk can be approximated as a truncated power law by following

n ( a ) = n tot ( p + 1 ) a max p + 1 − a min p + 1 a p , $$ \begin{aligned} n(a) = \frac{n_\mathrm{tot} (p+1) }{a_\mathrm{max} ^{p+1} - a_\mathrm{min} ^{p+1}}a^{p}\, , \end{aligned} $$(4)

where ntot is the total dust number density and amin and amax are the minimum and maximum grain sizes that truncate the distribution. Typical values of p are in the range of −4 to −2, depending on the physical details of the grain collisions.

Considering a discretized version of the size distribution with bins of size Δ a = a i + 1 2 − a i − 1 2 $ \Delta a=a_{i+\nicefrac{1}{2}} - a_{i-\nicefrac{1}{2}} $, we can write the same size distribution as

n ( a ) = ∑ i = 1 N n i ( p + 1 ) a i + 1 2 p + 1 − a i − 1 2 p + 1 a i p Θ ( a i + 1 2 − a ) Θ ( a − a i − 1 2 ) , $$ \begin{aligned} n(a) = \sum _{i=1}^N \frac{n_i\, (p+1) }{a_{i+\nicefrac {1}{2}}^{p+1} - a_{i-\nicefrac {1}{2}}^{p+1}}a_i^{p}\Theta (a_{i+\nicefrac {1}{2}}-a)\Theta (a - a_{i-\nicefrac {1}{2}}) \, , \end{aligned} $$(5)

where each bin i contains a total number density of ni and the size-grid cell interfaces are a i − 1 2 > a min $ a_{i-\nicefrac{1}{2}} > a_{\mathrm{min}} $ and a i + 1 2 < a max $ a_{i+\nicefrac{1}{2}} < a_{\mathrm{max}} $. The term Θ denotes the Heaviside step function and truncates every bin at a i + 1 2 $ a_{i+\nicefrac{1}{2}} $ and a i − 1 2 $ a_{i-\nicefrac{1}{2}} $. This assumes that the size distribution is a continuous power law with exponent p in each bin. With these definitions, we can rewrite the numerator of Eq. (3) as

∫ n ( a ) a 3 d a = p + 1 p + 4 ∑ i = 1 N n i a i + 1 2 p + 4 − a i − 1 2 p + 4 a i + 1 2 p + 1 − a i − 1 2 p + 1 ∑ i = 1 N n i χ i . $$ \begin{aligned} \int n(a)a^3\, \mathrm{d} a = \frac{p+1}{p+4} \sum _{i=1}^N n_i \frac{a_{i+\nicefrac {1}{2}}^{p+4} - a_{i-\nicefrac {1}{2}}^{p+4}}{a_{i+\nicefrac {1}{2}}^{p+1} - a_{i-\nicefrac {1}{2}}^{p+1}} \nonumber \sum _{i=1}^N n_i \chi _i \,. \end{aligned} $$

Likewise, we can write the denominator of Eq. (3) as

∫ n ( a ) a 2 d a = p + 1 p + 3 ∑ i = 1 N n i a i + 1 2 p + 3 − a i − 1 2 p + 3 a i + 1 2 p + 1 − a i − 1 2 p + 1 ∑ i = 1 N n i ξ i . $$ \begin{aligned} \int n(a)a^2\, \mathrm{d} a = \frac{p+1}{p+3} \sum _{i=1}^N n_i \frac{a_{i+\nicefrac {1}{2}}^{p+3} - a_{i-\nicefrac {1}{2}}^{p+3}}{a_{i+\nicefrac {1}{2}}^{p+1} - a_{i-\nicefrac {1}{2}}^{p+1}}\nonumber \sum _{i=1}^N n_i \xi _i\, . \end{aligned} $$

If the bin interfaces are fixed in time and if we assume constant p in every bin at all times, ξi and χi are constants. With these definitions, the Sauter mean of the entire size distribution can be written as

a s = ∑ i = 1 N n i χ i ∑ i = 1 N n i ξ i . $$ \begin{aligned} a_\mathrm{s} =\left.\sum \limits _{i=1}^N n_i \chi _i \right. \sum \limits _{i=1}^N n_i \xi _i \, . \end{aligned} $$(6)

The calculation of the Sauter mean radius in a hydrodynamic simulation of dust and gas is therefore reduced to the calculation of the dust number density ni for the given dust size bins. Knowledge of the Sauter mean allows us to calculate t thin NLTE $ t_\mathrm{thin}^{\text{NLTE}} $ from Eq. (2). If gas and dust are treated through radiative transfer, as in Muley et al. (2023), t thin NLTE $ t_\mathrm{thin}^{\text{NLTE}} $ can be used to dynamically relax the gas and dust temperatures with respect to each other. In a simpler use case, t thin NLTE $ t_\mathrm{thin}^{\text{NLTE}} $ can be used as a dynamically evolving cooling timescale of the gas via

T ( n + 1 ) = T eq + ( T ( n ) − T eq ) exp ( − Δ t t thin NLTE ) , $$ \begin{aligned} T^{(n+1)} = T_\mathrm{eq} + (T^{(n)}-T_\mathrm{eq} )\exp \left(-\frac{\Delta t}{t_\mathrm{thin} ^{\text{NLTE}}}\right)\, , \end{aligned} $$(7)

where Teq is the equilibrium temperature (set by stellar irradiation) and Δt is the current simulation time step. We demonstrate the latter case in hydrodynamic simulations of protoplanetary disks with a focus on the VSI activity in the next section.

2.2. Hydrodynamic simulations

We set up axisymmetric simulations in the r − ϑ plane of a protoplanetary disk in spherical coordinates. Our disk’s initial hydrostatic structure follows the standard accretion disk with mass Mdisk by Lynden-Bell & Pringle (1974), with a truncated power law in column density with exponent βΣ and characteristic radius Rc:

Σ ( R ) = ( 2 + β Σ ) M disk 2 π R c 2 ( R R c ) β Σ exp [ − ( R R c ) 2 + β Σ ] . $$ \begin{aligned} \Sigma (R)=(2+\beta _{\Sigma })\frac{M_{\text{disk}}}{2\pi R_\mathrm{c} ^2} \left(\frac{R}{R_\mathrm{c} }\right)^{\beta _{\Sigma }}\exp {\left[-\left(\frac{R}{R_\mathrm{c} }\right)^{2+\beta _{\Sigma }}\right]} \, . \end{aligned} $$(8)

From this, the vertical disk structure in z follows via

ρ = ρ mid exp [ ( H g R ) − 2 ( R R 2 + z 2 − 1 ) ] , $$ \begin{aligned} \rho = \rho _\mathrm{mid} \exp {\left[\left(\frac{H_\mathrm{g} }{R}\right)^{-2}\left(\frac{R}{\sqrt{R^2+z^2}} - 1\right)\right]} \, , \end{aligned} $$(9)

where we approximated the midplane dust density as ρ mid ≈ Σ ( r ) 2 π H g $ \rho_{\mathrm{mid}}\approx \nicefrac{\Sigma(r)}{\sqrt{2\pi} H_{\mathrm{g}}} $. The angular frequency in hydrostatic equilibrium follows accordingly as

Ω 2 ( R , z ) Ω K 2 = ( H g R ) 2 ( β T + β ρ − ( β Σ + 2 ) ( R R c ) β Σ + 2 ) − β T R R 2 + z 2 + β T + 1 , $$ \begin{aligned}&\frac{\Omega ^2(R,z)}{\Omega _\mathrm{K} ^2} = \left(\frac{H_\mathrm{g} }{R}\right)^2\left(\beta _T+\beta _\rho -(\beta _\Sigma + 2)\left(\frac{R}{R_\mathrm{c} }\right)^{\beta _\Sigma +2}\right) \nonumber \\&\qquad \qquad \qquad - \frac{\beta _T R}{\sqrt{R^2+z^2}} + \beta _T + 1 , \end{aligned} $$(10)

where βρ is the power law exponent of the radial midplane density profile ρ mid ∝ R β ρ exp [ − ( R R c ) 2 + β Σ ] $ \rho_{\mathrm{mid}}\propto R^{\beta_\rho}\exp\left[-\left(\nicefrac{R}{R_{\mathrm{c}}}\right)^{2+\beta_\Sigma}\right] $.

For the dust, we initialized the simulation with a Mathis et al. (1977) (MRN) distribution with four dust fluids, according to Eq. (5). We conducted simulations with three different initial dust grain sizes aini = 1 μm, 10 μm and 100 μm. As the grains undergo coagulation, we let the maximum grain size evolve in time following the monodisperse growth rate a ˙ max = a max ε 0 Ω K ( 1 − a max a lim ) $ \dot{a}_{\mathrm{max}}=a_{\mathrm{max}} \varepsilon_0\Omega_\text{K} (1-\nicefrac{a_{\mathrm{max}}}{a_{\mathrm{lim}}}) $ (Kornet et al. 2001), where ε 0 = Σ d , 0 Σ g , 0 $ \varepsilon_0=\nicefrac{\Sigma_{\mathrm{d,0}}}{\Sigma_{\mathrm{g,0}}} $ is the initial dust-to-gas column density ratio. Here, we have multiplied the original growth rate by the factor ( 1 − a max a lim ) $ (1-\nicefrac{a_{\mathrm{max}}}{a_{\mathrm{lim}}}) $ to bound the growth to the fragmentation limit. The equation then has an analytic solution for a given initial grain size aini:

a max ( t ) = a lim a ini e t ε 0 Ω K a lim + a ini ( e t ε 0 Ω K − 1 ) , $$ \begin{aligned} a_{\mathrm{max} } (t) = \frac{a_\mathrm{lim} a_\mathrm{ini} e^{t\, \varepsilon _0\, \Omega _\text{K}}}{a_\mathrm{lim} +a_\mathrm{ini} \left(e^{t\,\varepsilon _0\, \Omega _\text{K}} - 1\right)} \, , \end{aligned} $$(11)

where the limiting particle size alim = min(afrag, afrag-drift) is determined by the fragmentation velocity vfrag = 500 cm s−1, the radial pressure gradient, and the level of turbulence, characterized by the α parameter (see Birnstiel et al. 2012, for definitions). As our simulations were initialized without turbulence, we set this parameter to a very low value of 10−5. The particle growth was thus limited by the drift-fragmentation limit (∼1 cm) instead of the turbulent fragmentation limit. We did not change the mass content of the four dust-size bins during the growth. Only the bin interfaces and the respective mass-averaged particle size of each bin changed due to growth. Growth therefore also influences the thermal accommodation times. As mass is shifted to larger sizes, fewer small grains are present, and thus dust-gas collision rates decrease with time.

The four passive dust fluids were advected using the same technique as in Pfeil et al. (2023) and thus evolved under the influence of settling and radial drift in the terminal velocity approximation. The respective velocity components were calculated from the dust populations’ mass-averaged particle size, where an MRN size distribution was assumed within each bin at all times. The bin boundaries were defined to be equidistant in log space. The size grid was then defined between amin as the lower boundary of the first cell and amax as the upper boundary of the last cell.

The initial dust distribution of each population followed the same vertical structure as the gas, but instead of the gas scale height with the dust scale height,

H d = H g δ ini δ ini + St i , $$ \begin{aligned} H_\mathrm{d} = H_\mathrm{g} \sqrt{\frac{\delta _\mathrm{ini} }{\delta _\mathrm{ini} +\mathrm{St} _i}}\, , \end{aligned} $$(12)

where St i = π a i ρ m 2 Σ g $ \mathrm{St}_i=\nicefrac{\pi a_i \rho_{\mathrm{m}}}{2 \Sigma_{\mathrm{g}}} $ is the Stokes number of the respective dust fluid and δini = α is the initially assumed dust diffusivity. Simulations initialized with different particle sizes also have dust layers of different initial heights and cooling times. This can be seen in Fig. 1, where the top panels show the initial distribution of dust in the inner and outer disk regions as dashed lines.

thumbnail Fig. 1.

Dust-to-gas ratios and thermal relaxation times radially averaged over the 30–55 au inner region and the 115–150 au outer region of our simulations. As the VSI develops, it keeps the dust particle at high altitudes and thus retains the necessary cooling times in the inner disk. The outer regions, however, are mostly sedimenting, as the VSI is not active there, except for the case with the smallest initial particle size.

The cooling times increased with height since the dust density decreases with distance from the midplane. We show the cooling times in units of the local critical cooling time of the VSI. Values above one correspond to regions that are not susceptible to linear instability due to the vertical shear.

We ran the simulations for 1000 orbital timescales at 50 au, that is, 353 553 yr. We present the simulation parameters in Table 1.

Table 1.

Simulation parameters of the three presented runs.

3. Results

We present the radially averaged vertical dust distributions and the respective cooling times at the beginning and at the end of the simulation in Fig. 1, where we distinguish between an inner disk region (30–55 au) and an outer disk region (115–150 au). In Pfeil et al. (2023), we conducted simulations with static cooling times and showed that the VSI preferentially develops in the inner disks if the dust has already grown. We observed a qualitatively similar result here. The ongoing dust growth in the new simulations, however, modifies the results and imposes stricter conditions on whether the VSI can emerge. The turbulence is in every case stronger in the inner disk regions, where it maintains a relatively thick dust layer with an almost vertically constant dust-to-gas ratio ε and a stable vertical distribution of cooling times. The outer regions show more settling due to the longer growth timescale of the VSI. Continued settling during the simulation is, however, of minor importance for the cooling times compared to the coagulation process because thermal relaxation is mostly done by small particles that have very long settling times.

Whether the VSI develops at all strongly depends on the initial maximum particle size and thus the initial vertical distribution of small dust grains. If the particles are small at the beginning of the simulation, as in the case with aini = 1 μm, thermal accommodation times are short almost everywhere in the disk, and the VSI develops quickly. The resulting turbulence keeps the dust layer vertically extended and thus maintains the necessary cooling times self-sufficiently. This can also be seen in the first row of Fig. 2, which shows the dust distribution and the respective cooling times at the end of the simulation. Regions inside of ∼100 au show the typical filamentary VSI pattern in the dust densities and thus also in the cooling times. The outer regions develop slower, and the VSI is not present after 350 000 yr of evolution.

thumbnail Fig. 2.

Snapshots after ∼350 000 yr of our protoplanetary disk simulations with thermal accommodation times calculated from the dust distribution for three different initial particle sizes. Simulations initialized with larger particles and thus more settled dust layers are mostly not VSI active in the outer disk regions. If the VSI, however, can develop in the first place, it can stabilize the dust layer and sustain itself, even if grain growth is considered.

The situation is different in the simulation with initially large particles with aini = 100 μm. The particles are strongly settled at the beginning of the simulation, meaning that the cooling times in the upper disk atmosphere are accordingly long. The VSI begins to develop in the inner disk at a lower intensity than in the small-particle case. Nonetheless, a sufficiently thick layer of small dust is maintained, thus keeping the cooling times low enough in the inner regions.

The inner disk regions become VSI turbulent in all three cases, albeit on different timescales and at different intensities. The outer areas evolve differently. For small initial particle sizes, we observed that the VSI begins to develop at the end of the simulation. For larger particles, however, the outer disk areas are completely VSI inactive due to the already long cooling times at the beginning of the simulation. Regions beyond 70 au are not VSI active for the largest initial particle size of aini = 100 μm. For ten times smaller particles, VSI is active up to 90 au at the end of the simulation.

To visualize the impact of the VSI on the cooling times, we have plotted the height at which the critical cooling time for the VSI is reached (Eq. (1)) as a function of time and radius in Fig. 3. The beginning of all three simulations is characterized by the growth of the particles and the subsequent settling of the larger grains in the first ∼40 000 yr. During this time, the surfaces of the VSI-susceptible zone move toward the midplane of the disk. This is mostly due to the coagulation process, which transforms small, not sedimented grains into larger and quickly sedimenting ones. Cooling times are mostly determined by the small grains that remain in the atmosphere and only sediment slowly. The critical surface thus only moves toward the midplane because mass is transferred from small grains to big grains. During this initial particle growth, the first VSI modes develop in the inner regions of the disk and soon begin to corrugate the dust layer. The emerging VSI turbulence stabilizes the extent of the VSI-susceptible zone and even moves small grains up into the disk atmosphere, thus extending the susceptible region with time in the inner disk. In the simulation with the smallest initial particle size, the critical height reaches a stable value of ∼1.5 H in the regions close to the star. If the simulations are initialized with larger, and thus strongly sedimented particles, we observed that the initial vertical extent of the susceptible zone is already small. For aini = 10 μm, the critical height is stabilized at ∼1.2 H; for aini = 100 μm, it is at ∼0.8 H.

thumbnail Fig. 3.

Evolution of the vertical extend of the VSI susceptible region as a function of time and radius. In the inner disk, VSI quickly develops and stabilizes the settling dust layer in less than 50 000 yr. In the outer disk, VSI has not yet developed in our simulations. The thin lines show comparison simulations in which the VSI is artificially suppressed.

In the outer regions, VSI either develops very slowly (in the case of the smallest initial particle size of 1 μm) or not at all during the runtime of our simulations. This is a result of the limited simulation time and could also be a resolution issue, as the fastest growing modes become smaller with longer cooling times. The trend is nonetheless clear: If a hydrodynamic simulation is initialized with large (100 μm), settled particles (here for δini = 10−5), VSI only develops slowly, or possibly not at all, during the lifetime of a protoplanetary disk. If the particles are initially small, or vertically dispersed, VSI can develop quickly and even extend the susceptible region by diffusing small particles to the upper disk layers.

4. Discussion

Our studies allow for a first insight into the interplay of VSI and its influence on the cooling times, and vice versa. There are nonetheless various limitations to the presented approach. We only modeled the size distribution with four dust fluids and assumed each size bin to have an MRN particle size distribution. At the beginning of the simulation, the size distribution is still continuous. When the dust fluids begin to drift and sediment, this changes, and the complete size distribution is no longer an MRN distribution. In the future, the coagulation and fragmentation process should be considered when the dust populations evolve. In that way, a meaningful size distribution could be maintained during the simulations. We also set up the initial dust structure in settling-mixing equilibrium at all heights. This is strictly speaking only valid for Stokes numbers ≪1 because the terminal velocity approximation is used in the derivations (Dubrulle et al. 1995). Thus, at large distances from the midplane, our initial conditions might not be realistic. We modeled the dust as a passive fluid. Therefore, backreaction, which was shown to be able to quench the VSI activity in the disk midplane, especially if the dust settles, could not be accounted for (Schäfer et al. 2020). The settling, however, occurs mostly once the dust has grown. We thus do not expect it to be able to hinder the VSI growth in the inner disk.

We have not tested other initial size distribution power laws than the MRN distribution. If initially more small dust were present, the conditions could be more favorable for the VSI. The same is true if the fragmentation velocity is smaller, which counteracts the production of large grains.

Furthermore, we have modeled the dust size as a vertically global quantity. Simulations of coagulation that consider a vertically varying dust size and coagulation at all heights have shown that dust growth can start in the upper layers, where the relative sedimentation velocities are high. This sedimentation-driven coagulation depletes the upper layers of grains, which coagulate, fragment, and sediment continuously (Zsom et al. 2011). Thus, a full radial-vertical treatment of coagulation must be applied in future studies.

Other sources of cooling have also not been accounted for in our study. Especially in the upper disk atmosphere, cooling through molecular line emission could be the dominant channel of thermal relaxation (Woitke 2015). These considerations, however, require thermo-chemical modeling of the different molecular species and are currently beyond the scope of the presented work. A future self-consistent study of VSI with realistic thermal relaxation also has to involve the treatment of radiative transfer, ideally in a three-temperature framework, as recently presented by Muley et al. (2023). This is especially important when the effects of stellar irradiation and thermal relaxation in the optically thick regime should be included.

One of the caveats in the study of the VSI in protoplanetary disks remains the poorly constrained initial conditions of the simulations. Specifically, the state of the dust size distribution and the initial vertical extent of the dust layer are important for the onset of the VSI because they determine the cooling times. If the dust at the beginning of the disk’s lifetime is very small (∼1 μm), we would expect the VSI to develop within the first 350 000 yr of disk evolution, even in the outer regions. Small dust would be expected if the grains resemble interstellar dust or if an intense source of turbulence during the disk formation process has caused strong fragmentation. The initial phases might be gravitoturbulent, which could cause such high initial levels of turbulence (Gammie 2001; Johnson & Gammie 2003; Hirose & Shi 2017; Zier & Springel 2023). Some observations and studies, however, hint toward early grain growth (Galametz et al. 2019; Bate 2022). Although these studies do not suggest the complete removal of sub-μm-sized grains, a reduction in the respective number density would already lead to less favorable conditions for the VSI, as shown by our simulations.

5. Conclusions

We have for the first time conducted hydrodynamic simulations of protoplanetary disks with thermal accommodation times that are consistent with the simulated dust densities and the present grain size. This made it possible to assess the influence of different initial conditions and dust grain growth on the developing VSI. We find that the initial dust distribution has the biggest influence on the resulting spatial distribution of the turbulent gas flows. Coagulation also increases the cooling times and thus counteracts the VSI. This effect is, however, not a significant hurdle for the development of hydrodynamic turbulence in the inner regions of protoplanetary disks (inside of 70 au), even if the dust is already large (100 μm) and vertically settled at the beginning of the simulation. The VSI can develop during the initial growth phase and even extend the susceptible region by diffusing small grains vertically. In the outer regions of protoplanetary disks (beyond 70 au), the situation is more difficult for the VSI. If the dust is assumed to be large at the beginning of the simulation (100 μm) and is assumed to be in settling-mixing equilibrium (with a low initial diffusivity of δ = 10−5), we find that the potentially VSI-susceptible region is constrained to a small area around the midplane and will probably not develop VSI-induced turbulence during the disk’s lifetime.

We note that although our three simulations only differ in their initial dust grain size and reach identical maximum particle sizes, the resulting spatial distribution of dust and the level of VSI turbulence are vastly different between the three runs. Whether a specific disk stratification will develop VSI turbulence is thus highly dependent on the initial dust size distribution and the initial dust scale height.

Acknowledgments

The authors thank the anonymous referee for their helpful comments. T.P., H.K., and T.B. acknowledge the support of the German Science Foundation (DFG) priority program SPP 1992 “Exploring the Diversity of Extrasolar Planets” under grant Nos. BI 1816/7-2 and KL 1469/16-1/2. T.B. acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme under grant agreement No 714769 and funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under grants 361140270, 325594231, and Germany’s Excellence Strategy – EXC-2094 – 390783311. Computations were conducted on the computational facilities of the University of Munich (LMU) and on the c2pap cluster at the Leibniz Rechenzentrum under project pn36ta.

References

  1. Barranco, J. A., Pei, S., & Marcus, P. S. 2018, ApJ, 869, 127 [NASA ADS] [CrossRef] [Google Scholar]
  2. Bate, M. R. 2022, MNRAS, 514, 2145 [CrossRef] [Google Scholar]
  3. Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, 148 [Google Scholar]
  4. Burke, J. R., & Hollenbach, D. J. 1983, ApJ, 265, 223 [NASA ADS] [CrossRef] [Google Scholar]
  5. Dubrulle, B., Morfill, G., & Sterzik, M. 1995, Icarus, 114, 237 [Google Scholar]
  6. Flock, M., Turner, N. J., Nelson, R. P., et al. 2020, ApJ, 897, 155 [Google Scholar]
  7. Fukuhara, Y., & Okuzumi, S. 2024, PASJ, in press, https://doi.org/10.1093/pasj/psae042 [Google Scholar]
  8. Fukuhara, Y., Okuzumi, S., & Ono, T. 2023, PASJ, 75, 233 [NASA ADS] [CrossRef] [Google Scholar]
  9. Galametz, M., Maury, A. J., Valdivia, V., et al. 2019, A&A, 632 [Google Scholar]
  10. Gammie, C. F. 2001, ApJ, 553, 174 [Google Scholar]
  11. Hirose, S., & Shi, J. M. 2017, MNRAS, 469, 561 [NASA ADS] [CrossRef] [Google Scholar]
  12. Johnson, B. M., & Gammie, C. F. 2003, ApJ, 597, 131 [Google Scholar]
  13. Klahr, H., & Hubbard, A. 2014, ApJ, 788, 21 [Google Scholar]
  14. Klahr, H., Baehr, H., & Melon Fuksman, J. D. 2023, ApJ, submitted [arXiv:2305.08165] [Google Scholar]
  15. Kornet, K., Stepinski, T. F., & Rózyczka, M. 2001, A&A, 378, 180 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Lin, M. K. 2019, MNRAS, 485, 5221 [CrossRef] [Google Scholar]
  17. Lin, M. K., & Youdin, A. N. 2015, ApJ, 811 [Google Scholar]
  18. Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603 [Google Scholar]
  19. Malygin, M. G., Kuiper, R., Klahr, H., Dullemond, C. P., & Henning, T. 2014, A&A, 568, 91 [Google Scholar]
  20. Malygin, M. G., Klahr, H., Semenov, D., Henning, T., & Dullemond, C. P. 2017, A&A, 605, 30 [Google Scholar]
  21. Manger, N., Pfeil, T., & Klahr, H. 2021, MNRAS, 508, 5402 [NASA ADS] [CrossRef] [Google Scholar]
  22. Marcus, P. S., Pei, S., Jiang, C. H., et al. 2015, ApJ, 808, 87 [NASA ADS] [CrossRef] [Google Scholar]
  23. Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425 [Google Scholar]
  24. Muley, D., Fuksman, J. D. M., & Klahr, H. 2023, A&A, 678, A162 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610 [Google Scholar]
  26. Pfeil, T., & Klahr, H. 2021, ApJ, 915, 130 [NASA ADS] [CrossRef] [Google Scholar]
  27. Pfeil, T., Birnstiel, T., & Klahr, H. 2023, ApJ, 959, 121 [NASA ADS] [CrossRef] [Google Scholar]
  28. Probstein, R. 1969, in Probl. Hydrodyn. Contin. Mech. Contrib. Honor Sixtieth Birthd. Acad., eds. L. I. Sedov, & M. Lavrent’ev (Philadelphia: Society for Industrial and Applied Mathematics), 568 [Google Scholar]
  29. Sauter, J. 1926, VDI-Verlag, Forschungsarbeiten auf dem Gebiete des Ingenieurwesens, 279 [Google Scholar]
  30. Schäfer, U., & Johansen, A. 2022, A&A, 666, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Schäfer, U., Johansen, A., & Banerjee, R. 2020, A&A, 635 [Google Scholar]
  32. Stoll, M. H., & Kley, W. 2016, A&A, 594, 57 [Google Scholar]
  33. Urpin, V. 2003, A&A, 404, 397 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Urpin, V., & Brandenburg, A. 1998, MNRAS, 294, 399 [Google Scholar]
  35. Woitke, P. 2015, EPJ Web Conf., 102, 00011 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. Zier, O., & Springel, V. 2023, MNRAS, 520, 3097 [NASA ADS] [CrossRef] [Google Scholar]
  37. Zsom, A., Ormel, C. W., Dullemond, C. P., & Henning, T. 2011, A&A, 534, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

All Tables

Table 1.

Simulation parameters of the three presented runs.

In the text

All Figures

thumbnail Fig. 1.

Dust-to-gas ratios and thermal relaxation times radially averaged over the 30–55 au inner region and the 115–150 au outer region of our simulations. As the VSI develops, it keeps the dust particle at high altitudes and thus retains the necessary cooling times in the inner disk. The outer regions, however, are mostly sedimenting, as the VSI is not active there, except for the case with the smallest initial particle size.
In the text
thumbnail Fig. 2.

Snapshots after ∼350 000 yr of our protoplanetary disk simulations with thermal accommodation times calculated from the dust distribution for three different initial particle sizes. Simulations initialized with larger particles and thus more settled dust layers are mostly not VSI active in the outer disk regions. If the VSI, however, can develop in the first place, it can stabilize the dust layer and sustain itself, even if grain growth is considered.
In the text
thumbnail Fig. 3.

Evolution of the vertical extend of the VSI susceptible region as a function of time and radius. In the inner disk, VSI quickly develops and stabilizes the settling dust layer in less than 50 000 yr. In the outer disk, VSI has not yet developed in our simulations. The thin lines show comparison simulations in which the VSI is artificially suppressed.
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.