Dust Coagulation Reconciles Protoplanetary Disk Observations with the Vertical Shear Instability. I. Dust Coagulation and the VSI Dead Zone

Protoplanetary disks exhibit a vertical gradient in angular momentum, rendering them susceptible to the vertical shear instability (VSI). The most important condition for the onset of this mechanism is a short timescale of thermal relaxation (≲0.1 orbital timescales). Simulations of fully VSI active disks are characterized by turbulent, vertically extended dust layers. This is in contradiction with recent observations of the outer regions of some protoplanetary disks, which appear highly settled. In this work, we demonstrate that the process of dust coagulation can diminish the cooling rate of the gas in the outer disk and extinct the VSI activity. Our findings indicate that the turbulence strength is especially susceptible to variations in the fragmentation velocity of the grains. A small fragmentation velocity of ≈100 cm s−1 results in a fully turbulent simulation, whereas a value of ≈400 cm s−1 results in a laminar outer disk, being consistent with observations. We show that VSI turbulence remains relatively unaffected by variations in the maximum particle size in the inner disk regions. However, we find that dust coagulation can significantly suppress the occurrence of VSI turbulence at larger distances from the central star.


INTRODUCTION
Around 1 % of the mass of protoplanetary disks is initially composed of solids (Lodders 2003;Magg et al. 2022).Despite its small contribution to the overall mass budget, this dust is the building material for planetesimals and planets and an essential observable for infrared and radio observations.It can have a considerable influence on the gas dynamics within the disk via drag forces (Weidenschilling 1980;Youdin & Goodman 2005) and is the main source of opacity.Therefore, cooling and heating are mostly determined by the solids for the bulk of the disk (Semenov et al. 2003;Woitke 2015;Malygin et al. 2017).Many linear instabilities of the gas flow depend on the local rate of thermal relaxation (Klahr & Bodenheimer 2003;Petersen et al. 2007a,b;Klahr & Hubbard 2014;Lin & Youdin 2015;Marcus et al. 2015;Lyra & Umurhan 2019) or the ionization state of the gas (Balbus & Hawley 1991;Blaes & Balbus 1994), and are therefore sensitive to the assumed dust size distribution (Barranco et al. 2018;Fukuhara et al. 2021;Kawasaki & Machida 2023).
In this work, we are specifically interested in the evolution of the vertical shear instability (VSI, Urpin & Brandenburg 1998), which requires a short thermal relaxation time of the gas (Lin & Youdin 2015;Manger et al. 2021;Fukuhara et al. 2021).VSI was studied in much detail in isothermal and adiabatic disk models at various rates of β cooling (e.g., Nelson et al. 2013) and in models with radiative transfer (e.g., Stoll & Kley 2016;Stoll et al. 2017;Flores-Rivera et al. 2020).Due to the numerical obstacles of incorporating dust evolution models in hydrodynamic simulations (Drążkowska et al. 2014;Gonzalez et al. 2017;Drążkowska et al. 2019;Lombart et al. 2022), most previous studies consider a static dust population, perfectly coupled to the gas.These studies often aim for a detailed analysis of the instability mechanism itself (e.g., Nelson et al. 2013;Manger et al. 2021;Svanberg et al. 2022).They showed the VSI's ability to cause large-scale vortex formation (Richard et al. 2016;Manger & Klahr 2018;Pfeil & Klahr 2021) and strong corrugations in the dust layer (Stoll & Kley 2016;Flores-Rivera et al. 2020).Simulations assuming perfectly coupled dust or isothermal conditions cannot, however, model the conditions in real protoplanetary disks, for which observations show an evolved dust population (Pérez et al. 2012;Tazzari et al. 2016;Huang et al. 2018;Ohashi & Kataoka 2019;Sierra et al. 2021), substructures (ALMA-Partnership et al. 2015;Andrews et al. 2018;Dong et al. 2018), and planets (Keppler et al. 2018).In this work, we intend to go one step further by considering an evolved-yet static-dust population in two-dimensional simulations of smooth protoplanetary disks.
Our work is motivated by the results of Dullemond et al. (2022), which show that VSI turbulence in an isothermal disk model is not consistent with observations of thin dust layers in protoplanetary disks.In Pfeil & Klahr (2021), we have explored the impact of a more realistic cooling time prescription on the strength of VSI turbulence.For this, we assumed the presence of a static, µm-sized dust population in the inner parts of a protoplanetary disk (at ∼ 5 au).For these setups, we found that the collisional decoupling of the gas and dust particles inhibits thermal relaxation in the disk atmosphere and thus reduces VSI turbulence.The respective collisional coupling time scale depends on the size distribution and is, thus, sensitive to the fragmentation velocity and other dust properties.Fukuhara et al. (2021) further studied this effect in models with a more detailed prescription of the dust size distribution.They found that coagulation can indeed inhibit the VSI by depleting the number of small grains that provide radiative cooling.In their most recent study, Fukuhara et al. (2023) attempted to simulate this in a more selfconsistent way, by taking into account the effect of the VSI on the diffusivity and the cooling times.Since they could not afford to dynamically evolve the dust population within their hydrodynamic simulations, they relied on analytic prescriptions for the cooling time for a static dust size distribution.
In this work, we study the effect of a more realistic steady-state dust distribution for varying coagulation parameters using DustPy (Stammler & Birnstiel 2022) and PLUTO (Mignone et al. 2007).We deduce thermal relaxation times from dust coagulation models in Section 3 which are then implemented in hydrodynamic simulations, from which we study the VSI activity in Section 4. This makes it possible to study the influence of dust coagulation and the coagulation parameters on VSI turbulence.These steps are schematically displayed in Figure 1.In the next step, we introduce passive dust fluids to our simulations in Section 4.1 to study the effect of the emerging VSI turbulence on the thickness of the dust layer.To make our results comparable to observations, we create synthetic intensity maps with RADMC-3D (Dullemond et al. 2012) in Section 5.

Cooling Requirements for the Vertical Shear Instability
Vertical shear, in the geophysical context also known as thermal wind (Holton & Hakim 2012), is a consequence of the radial temperature gradient in the vertically stratified protoplanetary disks.The temperature gradient itself is maintained by stellar irradiation.Consequently, fluid parcels can be displaced upward into a region of lower specific angular momentum experiencing an outward acceleration.A perturbation along such a trajectory violates Rayleigh's stability criterion and leads to a continued acceleration of the fluid parcel.This mechanism is called the Vertical Shear Instability (Urpin & Brandenburg 1998) and results in vertically elongated and radially narrow flow patterns.However, as the gas parcels enter the lower-density regions of the disk atmosphere, they are subjected to buoyancy forces, which, in a stably stratified atmosphere, would lead to an oscillation around the disk midplane.The characteristic frequency of this oscillation is the Brunt-Väisälä frequency where z is the distance from the disk midplane, ρ g is the gas density, P is the pressure S is the gas entropy, and C P is the gas' specific heat capacity at constant pressure.Thermal relaxation counteracts the restoring force of this oscillation by adjusting a gas parcel's specific entropy to the background.In order for the vertical shear to overcome buoyancy and trigger the VSI, thermal relaxation must be fast.Lin & Youdin (2015) have shown that vertically global VSI grows the fastest if the cooling timescale fulfills where R is the distance to the central star, β T is the power law exponent of the temperature profile, H g is the pressure scale height, Ω K is the local Keplerian frequency, and γ = C P /C V is the gas' heat capacity ratio.Equation 2 was derived under the assumption of a vertically constant thermal relaxation time.As we specifically consider the height dependence of thermal relaxation, we will use the local definition of a critical cooling time (Urpin 2003;Klahr et al. 2023) for local VSI modes Planck mean opacity (averaged over the size distribution) Settling-mixing equilibrium Equations ( 15) and ( 16) dsharp_opac opacity model  P (, ) Dust size distribution left side of Figure 1  thin NLTE (, ) Equations ( 5)-( 8) Fit with Equation (B17)

PLUTO hydrodynamic simulations
Cooling time (Equation 18):  thin NLTE (, ) Four dust fluids: Resolution: Discrete dust size distribution Equation ( 20)  In fact, numerical studies like Manger et al. (2021) investigated the dependency of the VSI turbulence on a vertically constant thermal relaxation time and found VSI not to develop for cooling times beyond the critical value for global modes.This may be due to numerical resolution, as Lin & Youdin (2015) show that VSI exists for all cooling times, yet at reduced efficiency.Urpin (2003) derived growth rates in this regime, which show a decay proportional to t −1 c .This behavior was recently confirmed in high-resolution 1 studies of the VSI and other thermal instabilities in disks by Klahr et al. (2023).It is still subject to investigation how longer growth times will translate into turbulence levels for the non-linear regime, especially in terms of angular momentum transport, diffusion, and gas rms velocities.The saturation behavior of VSI and other thermal baroclinic instabilities especially for longer cooling times at sufficient resolution is still being investigated (Latter & Papaloizou 2018;Cui & Latter 2022;Klahr et al. 2023).

Optically Thin Thermal Relaxation
Thermal relaxation of the gas in a protoplanetary disk is mostly achieved via thermal coupling with the dust in a two-stage process.At low temperatures, the emission timescale of the gas molecules is long, which means that cooling is only possible via thermal accommodation with the strongly emitting dust particles through collisions.
1 PLUTO-4.2simulation with 256 cells per gas scale height, WENO reconstruction, and RK3 time integration (Klahr et al. 2023).Barranco et al. (2018), derived the thermal relaxation times for the non-LTE case between dust grains and the gas based on the calculation of cooling rates (see Appendix A for a recap of the derivations).For a given dust size distribution n(a), the Sauter mean radius is an instructive parameter in this context, defined as (Sauter 1926) where the size integral is executed over the entire size distribution.Corresponding to the Sauter mean, we define a respective number density n S = ρ d / 4 /3 πρ m a 3 S and a collisional cross-section σ s = πa 2 S , where ρ m = 1.67 g cm −3 is the interior density of the dust grains.With these definitions, we write the thermal accommodation timescale for the gas molecules and the dust grains (Probstein 1969;Burke & Hollenbach 1983) as where vg = c s 8 /π is the average gas molecule velocity of a Maxwell-Boltzmann distribution with the isothermal speed of sound c s .Similarly, a timescale for the thermal relaxation of the dust component can be derived, which reads with the dust-to-gas density ratio ρ d/ρ g = ε and the specific heat capacity of the dust particles C d .As a typical value we pick C d = 800 J kg −1 K −1 , as used by Barranco et al. (2018) (see Wasson 1974;Piqueux et al. 2021;Biele et al. 2022).If the collisional coupling is efficient, i.e., temperature perturbations in the gas are transferred to the dust, the thermal equilibrium of the grains will be restored by the emission of radiation.This happens on the black body timescale, depending on the dust density distribution ρ d (a) in units of [g/cm 4 ] and the respective Planck mean opacity distribution κ P (a, T ), in units of [cm 2 /g] with the Stefan-Boltzmann constant σ SB .The total thermal relaxation time of the dust gas mixture can then be calculated following Equation ( 19) from Barranco et al. ( 2018) . In practice, this means the slowest channel of energy transfer acts as a bottleneck and the longest timescale of thermal relaxation determines the cooling time scale of the gas.If the dust's emissivity is low, energy cannot be emitted effectively by the grains, and temperature perturbations cannot decay, no matter how well the grains and molecules are coupled (t NLTE thin ≈ t rad d ).This situation is unlikely to occur in protoplanetary disks because of the large dust opacities.Another case is the collisional decoupling of dust grains and gas molecules.At low densities and in regions where small grains are depleted, heat cannot be transferred between the main carriers of thermal energy (the gas molecules) and the emitters (the dust grains).The high emissivity of the grains does not matter in such a case, since temperature perturbations stay locked in the poorly emitting gas (t NLTE thin ≈ t coll g ).Muley et al. (2023) introduced a three-temperature radiation transport scheme, which treats dust and gas temperatures separately, yet coupled via collisions.They also find that in most cases the collisional time scale is the most relevant to determine thermal relaxation.
In this case, the cooling time is proportional to the square root of the maximum particle size.This can be shown by assuming the size distribution to be a truncated power law with maximum particle size a max , minimum size a min , and power law exponent p = −3.5.Then a s = √ a max a min and thus t coll g ∝ (n S σ S ) −1 ∝ √ a max .Sticking collisions between grains typically in-crease the maximum particle size until a fragmentationcoagulation equilibrium is reached.In this case, a max ≈ a frac ∝ v 2 frag holds (Birnstiel et al. 2012), and we deduce that the collisional timescale is directly proportional to the fragmentation velocity in this case.Laboratory experiments aim to determine the actual value of v frag which is dependent on the composition and porosity of grains (Blum 2000;Wurm et al. 2001;Blum et al. 2006;Musiolik & Wurm 2019).Typical values lie within a range of 100-1000 cm s −1 .
An additional uncertainty arises from the unknown relative grain velocities, which depend on the strength of turbulence, differential drift, and settling.Especially the strength of turbulence in protoplanetary disks is highly uncertain and also a subject of this article.The simplest assumption for the turbulent transfer of energy across length scales is the Kolmogorov cascade.For the resulting energy spectrum, relative grain velocities can be approximated as δv ≈ √ 3αStc s (Ormel & Cuzzi 2007), with the Stokes number St (see Equation 14).This is the underlying assumption for the derivation of a frag .In this turbulence prescription, which is Based on the assumption of a mixing length model (Prandtl 1925), turbulent stresses result in an effective viscosity where c s is the local sound speed and H is the pressure scale height of the disk (Shakura & Sunyaev 1973).From this, turbulent rms velocities can be related to α by assuming a turbulent correlation time of With this, a frag ∝ α −1 , implying t coll g ∝ α − 1 /2 .Low α therefore corresponds to longer cooling times, as a consequence of the presence of larger particles.Additionally, lower levels of turbulence correspond to smaller dust scale heights, leading to a depletion of the upper layers and an additional dampening of the VSI in these regions.Fukuhara et al. (2021) investigated the effect of varying maximum particle sizes throughout a protoplanetary disk and found that the presence of VSI depends on particle sizes via the cooling time dependency.
In the following sections, we investigate this effect through the use of more realistic dust coagulation models and subsequent hydrodynamic simulations.We aim to determine the implications for the interpretation of observational data and the respective feedback onto the dust layer by turbulent mixing through the VSI.

DustPy COAGULATION MODELS
In the previous sections, we discussed the importance of thermal relaxation for the VSI.We have also highlighted that the cooling times are highly sensitive to the present dust population, most importantly, the maximum particle size.
In this section, we present a series of dust coagulation simulations, conducted with DustPy, that further illustrate the impact of dust coagulation on the cooling times.We use the output of these simulations to calculate cooling time distributions for our subsequent hydrodynamic simulations with the PLUTO code.
For our disk model we employ the standard Lynden-Bell & Pringle (1974) profile for a solar-mass star and a 0.05 M ⊙ disk with dust-to-gas ratio (metallicity) Z = 0.01 (see Table 1) We set the radial column density gradient to β Σ = −0.85,and the cutoff radius to R c = 100 au.Our radial temperature profile is determined by passive stellar irradiation and assumed to be constant in the vertical direction (see Chiang & Goldreich 1997;D'Alessio et al. 1998;Dullemond et al. 2018) where L * is stellar luminosity, and φ = 0.02 is the flaring angle.Gas evolution and dust drift alter the dust size distribution in protoplanetary disks.The overall effect of these transport phenomena on the shape of the distribution is, however, most relevant in the final stages of disk evolution, when the growth front has reached the outer disk edge and the mass budget is quickly decreasing (i.e., when the dust accretion rate is no longer radially constant, Birnstiel & Andrews 2014).At what point in time after disk formation this becomes relevant is dependent on the disk's size, its radial structure, the dust-to-gas ratio, the strength of turbulence, the fragmentation velocity, etc.In this study, we are interested in the effect of dust coagulation on the cooling times and, through the cooling times, on the VSI.In the inner parts of the disk, a steady state distribution, determined by fragmentation and coagulation, will be reached and approximately maintained as long as the outer disk edge is not yet moving inward.We have therefore decided to completely disregard any transport effects (except the vertical settling-mixing equilibrium).We are thus calculating a steady state dust distribution for each parameter set that is only determined by fragmentation and coagulation.The output of our models is, therefore, time-independent once the equilibrium size distribution is reached at each radius.In that way, we avoid selecting an arbitrary simulation snapshot.Note that this is still an idealized assumption.In reality, radial drift and gas evolution could slightly alter the radial structure and the size distributions at similar timescales.Typically, drift-limited size distributions are slightly steeper than in the fragmentation limit (Birnstiel et al. 2011).In recent studies, the VSI itself was also shown to alter the radial disk structure (Manger et al. 2021).
Our DustPy models are run for 10 5 yr, after which coagulation-fragmentation equilibrium is reached at every radial grid cell.
We conduct simulations for three different fragmentation velocities v frag = 100, 200 and 400 cm s −1 and for a turbulence parameter α = 10 −3 .Additionally we probe two different turbulent diffusivities with α = 10 −4 and 10 −2 , at v frag = 100 cm s −1 .At this point we do not further specify the origin of the diffusivity α, making it a free parameter for the coagulation models.We show the resulting dust size distribution at 50 au and 100 au on the left-hand-side of Figure 2 and some key particle properties are shown in Table 1.We can see that the particles grow to larger sizes at smaller distances to the central star, in accordance with analytic estimates of the fragmentation-limited particle size (Birnstiel et al. 2012).The respective size distributions can be approximated with power laws with exponents −p ≈ 3.6 − 3.7.These values lie within the typical range for fragmentation-limited size distributions derived by Birnstiel et al. (2011).

Thermal Relaxation Times Derived from Dust Coagulation Simulations
We derive the vertical structure from these, vertically integrated, DustPy models by assuming vertical hydrostatic equilibrium for the gas and vertical settlingmixing equilibrium for the dust.Gas densities thus follow with gas scale height H g and ρ mid = Σ(R) / √ 2πH 2 g .We assume an ideal equation of state P = ρc 2 s .The vertical dust distribution is determined by the diffusion parameter δ and the Stokes number of the individual size bins on the size distribution, which is defined as Volume dust densities for each size are then derived by calculating the dust scale height The resulting temperature and density structure is used to calculate the Planck mean opacities of the dust.We use the DSHARP opacity model by Birnstiel et al. (2018) as implemented in the dsharp_opac python package with the standard DSHARP particle properties.Thermal relaxation times of the gas can then be calculated from the disk structure and opacities via equations ( 5)-( 8).For the given parameters in our simulations, we find that the thermal relaxation time is limited by the collision timescale outside of ∼ 10 au.At smaller radii, the disk might become optically thick, meaning the relaxation time of temperature perturbations depends on the respective length scale.We are therefore only modeling the parts of the disk around 50 au, where thermal relaxation operates in the optically thin regime.Figure 2 shows the size distributions and the vertical profile of the thermal relaxation times for the respective coagulation and turbulence parameters at 50 au and 100 au.We find that the cooling times increase with height above the midplane.The reason for this is that cooling is achieved via collisions between dust particles and gas molecules, which become rarer at lower densities.This also means that models with larger particles have longer thermal relaxation times because of the reduced number densities of dust particles and the stronger settling.Higher fragmentation velocities are counteracting the VSI.Likewise, models with weaker turbulence parameter α can also be expected to have less VSI activity, as demonstrated by our numerical simulations.

PLUTO SIMULATIONS BASED ON COAGULATION MODELS
We set up axisymmetric PLUTO simulations with the same radial structure as our DustPy models to study the evolution of VSI with the respective model's cooling times.Pressure forces act in the outward direction of the disk and therefore decrease the equilibrium rotation frequency of the gas, especially at the steep outer edge of the disk.We define our hydrostatic initial rotation profile accordingly as where β ρ is the power law exponent of the midplane gas density ρ ∝ R βρ and β T is the power law exponent of the radial temperature profile T ∝ R β T .Thermal relaxation is realized as in Pfeil & Klahr (2021), by analytically relaxing the gas pressure toward the equilibrium profile (determined by stellar irradiation).Density is kept constant in this cooling step, which makes a relaxation in pressure equal to a relaxation in temperature for an ideal equation of state.
where (n) denotes the number of the current simulation timestep of length ∆t.The equilibrium temperature T eq is defined by stellar irradiation (Equation 12).Cooling times, presented in the previous section, are derived from DustPy simulations (see Figure 2) and subsequently fitted as a function of local gas density and temperature for each simulation (for a detailed description of the fits, see Appendix B).
Fitting the spatial distributions of the thermal relaxation times as functions of density and temperature also introduces uncertainties in the cooling times for PLUTO.For all models except one, these errors lie within 25 % with respect to the real distribution of cooling times.For the case of the most settled particles (v frag = 100 cm s −1 , α = 10 −4 ), however, the fitting function seems to diverge further from the real distribution and the fit deviates up to 58 % from the cooling times close to the midplane.This is likely due to the difference between this particular highly settled model and the other less settled cases.Since the cooling times vary over several orders of magnitude throughout the simulation domain and between the models, we deem this uncertainty acceptable-also because the overall distribution of cooling times is still well reproduced (this can be seen in the matching contours in Figure 12).It is worth noting, however, that in this work, we only study the overall trends of VSI turbulence with the coagulation parameters and do not aim to exactly reproduce specific systems or observations.
Cooling time fit tNLTEΩK as used in PLUTO at 50 au at 100 au Dust size distributions at 50 au (solid lines) and 100 au (dashed lines) of our DustPy models (left side).On the right-hand side, we show the respective vertical cooling time profiles, assuming vertical settling-mixing equilibrium for the given α and the critical VSI cooling time.Models with larger particles also exhibit longer cooling times due to collisional decoupling between dust and gas.We also show the height-dependent cooling time for local VSI modes as purple lines (see Equation 3).
The resulting analytic cooling time prescriptions are used within our PLUTO simulations to calculate t NLTE thin from the local disk structure.Since cooling is dominated by the small grains, which predominantly move along with the gas, minor disturbances in the gas densities, as caused by the VSI, can also influence the cooling times in this model.We emphasize that this is a minor effect in our simulation, and does not have an impact on the resulting turbulence.It should be noted, that our cooling time prescription, which is derived from dust coagulation models, is static throughout the simulation.
Although our coagulation models assumed a certain turbulent diffusivity delta to calculate relative particle velocities, we set up our hydrodynamic simulations to be inviscid.This is because we want to study the onset of the VSI and the resulting turbulence strength.Applying the same diffusivities as for the coagulation models (δ = 1 × 10 −4 − 1 × 10 −2 ) as viscosity in PLUTO would likely stop the VSI from emerging in the first place (Barker & Latter 2015).Note that setting up viscous simulations would also not be fully self-consistent since relative particle velocities in DustPy are inferred from perfectly isotropic turbulence and the resulting Kolmogorov cascade, which would not be the case for the developing VSI turbulence in our simulations.
We carry out the calculations for 500 orbital periods at 50 au (= 176 777 yr).Simulation domains are set up in spherical coordinates and extend from 25-150 au in the radial direction, and over ±3 pressure scale heights from the midplane of the disk in the polar direction.We resolve one scale height at 50 au with 85 cells and employ logarithmic griding in the radial direction to preserve the cells' aspect ratios, resulting in a 2011 r × 513 ϑ grid.Pe-riodic boundary conditions are set up in the azimuthal direction with only one grid cell, making our simulations axisymmetric.Radial and polar boundaries are setup up as reflective for the orthogonal velocity components and as zero-gradient for the respective tangential velocity components.Pressure and density in the boundary cells are kept at the initial condition.
In Figure 3, we show the vertical velocities in our simulations at the end of the simulation time.It is evident, that the spatially varying cooling times set constraints on where the VSI can be active and where vertical motions are suppressed by buoyancy.As a comparison, we also show an isothermal simulation (i.e., ideal VSI), in which the resulting turbulence is present in the entire simulation domain and at higher turbulent Mach numbers.For the case of v fr = 400 cm s −1 and α = 10 −3 , we find the disk to be completely quiescent outside of ∼ 80 au, due to the long cooling times.In this case, dust would settle into a very thin layer in the outer disk, which we will further investigate in the next sections.Similarly, the disk regions outside of ∼ 100 au show only very little VSI activity for the coagulation model with v fr = 100 cm s −1 and α = 10 −4 .
To characterize the development and strength of the VSI turbulence, we measure the Favre-averaged (i.e.density-weighted) turbulent Mach numbers over the whole simulation domains, where the average in a direction x (polar, radial, or both) is defined as The isothermal run shows a snapshot after only 200 orbits.White contours mark the position at which the critical cooling time for the VSI is reached (Equation 3), i.e., VSI is theoretically possible within the white lobes.
where v r and v ϑ represent the radial and polar velocity components.Since our simulations are set up hydrostatically, these components measure turbulent fluctuations caused by the VSI.While velocities in our isothermal simulation saturate after ∼ 100 orbits at ⟨M⟩ ≈ 4 × 10 −2 , all other, non-ideal simulations, reach lower Mach numbers and have longer growth time scales (see Figure 4).The vertical profile of the Mach numbers shows the typical vertical increase and a sharp upper cutoff, similar to the results in Pfeil & Klahr (2021).The collisional decoupling of dust particles and gas molecules is the reason for this behavior.Figure 4 also shows the three Mach numbers corresponding to the diffusivities chosen to calculate turbulent relative velocities between particles in our coagulation model (α= 1×10 −4 , 1×10 −3 and 1 × 10 −2 ).As can be seen, the three lines do not exactly correspond to the measured Mach numbers of our simulations.This is, however, also not to be expected, since the direct conversion between Mach numbers and particle collision speed (see Equation 10) assumes a perfect Kolmogorov spectrum and, thus, isotropic turbulence which is not given for the VSI.The calculation of collision speeds would furthermore depend on the correlation time spectrum which was not taken into account here.We conclude that the level of VSI turbulence is highly dependent on the physical details of the dust coagulation process.If dust grains can grow up to the fragmentation limit-which is to be expected in most parts of protoplanetary disks in the early evolutionary stages-we can expect weak collisional coupling between dust grains and gas molecules in the optically thin, outer regions, leading to inefficient cooling and only weak VSI turbulence.The magnitude of the impact of dust coagulation on the hydrodynamic turbulence depends mostly on the maximum size of the grains, where larger grains correspond to less cooling and, thus, stronger damping of VSI.

Dust Dynamics in the PLUTO Simulations
In the previous section, we have shown that the VSI activity in protoplanetary disks is highly sensitive to the properties of the present dust grain population, especially the largest grain size.However, we cannot directly infer the VSI's feedback on the dust population.Dullemond et al. (2022) have clearly shown that the ideal VSI is inconsistent with the observed thickness of protoplanetary disks in millimeter-wave observations with ALMA (Villenave et al. 2020(Villenave et al. , 2022)).Our simulations In models with larger particles, cooling times are generally longer, which results in lower growth rates and lower Mach number turbulence.The vertical profiles on the right-hand side change accordingly.Cooling times in models with larger maximum particle size increase more rapidly with height above the midplane, which also cuts off the VSI turbulence.Isothermal models typically have vertically increasing turbulent velocities.The three dashed horizontal lines show the Mach numbers corresponding to the three α values that we assumed for our coagulation models (see Equation 10).Note that the conversion between turbulent Mach numbers and diffusivities assumes a perfect Kolmogorov turbulence spectrum (see discussion in Section 6), which is likely not given for the anisotropic VSI turbulence.10).Note that such a conversion assumes a perfect Kolmogorov turbulence spectrum (see discussion in Section 6).
show that the level of turbulent vertical velocities can vary by orders of magnitude across the disk, depending on the details of the dust size distribution.
In this section, we explore how these different levels of turbulence impact the thickness of the dust layer.For this, we restart the simulations after the VSI has reached a saturated level of turbulence.We add four dust fluids, resembling a power law size distribution n(a) ∝ a p , and thus Σ d (a) ∝ n(a)m(a) ∝ a p+3 .Normalizing to the total dust column density (column dust-to-gas ratio Z = 0.01) and integrating the distribution over the size bin i with boundaries a i and a i+1 , we get The maximum grain sizes a max and exponents p are derived from the underlying DustPy models (measured at a distance of 50 au as the size including 99.9 % of the dust mass, see Table 1).Similar to the DustPy simulations, the minimum grain size is set to 0.1 µm, which is a typical size assumed for monomers in protoplanetary disks (Tazaki & Dominik 2022) and which is constant throughout the simulations.We divide the power law size distribution into four sections, equally spaced in logarithmic size space between a min and a max .The initial vertical dust distribution is determined by the midplane Stokes numbers and the level of turbulence assumed in the respective DustPy runs, follwing Equation 16.Dust is allowed to flow in from the outer boundary of the simulation domain with the initial vertical distribution.As to the time of this work, the PLUTO code has no builtin dust fluids.Therefore, we make use of the available gas tracer fluids.To model radial dust drift and vertical settling we modify the tracer fluxes according to the respective grain sizes' relative velocity to the gas, which is given by the prescriptions of Nakagawa et al. (1986) (terminal velocity approximation).Each dust fluid is advected with the gas velocity plus the drift correction of the mass-averaged size of the respective size bin.In Appendix C we present tests of this method that verify its accuracy.
We continue the previous, gas-only, VSI simulations with dust for another 150 orbits (measured at 50 au).Figure 6 depicts the distribution of dust-to-gas ratios in our simulations after 150 orbits.In our model with α = 10 −3 and v fr = 400 cm s −1 , we have the largest particles of ≈ 0.14 cm radius, while the smallest particles are present in the model with α = 10 −2 and v fr = 100 cm s −1 , with a maximum size of ≈ 15 µm (see Table 1).As a comparison, we initialize the isothermal simulation with the largest grains, to get an estimate of the effect of ideal VSI on a grown dust population (as in Dullemond et al. 2022).The effect of the different levels of VSI turbulence, depending on the coagulation parameters and the respective thermal relaxation times becomes visible in the dust-to-gas ratios, where the simulations with larger particles, longer cooling times, and less VSI turbulence have more settled dust layers.Especially the outer disk regions are affected by this, as can be seen in the cases with v fr > 100 cm s −1 and α < 10 −2 .
We can furthermore see, that the isothermal simulation provides a good approximation for the models with the smallest particles.This is to be expected because the models with the smallest particles also have the shortest cooling times, making the VSI modes almost isothermal.
To visualize the clear distinction between the inner VSI active region and the outer VSI inactive regions, we plot the time and radially averaged total dust-to-gas ratios in Figure 7.For the models with fully VSI active disks, we find flat top, or double-peaked dust distributions throughout the entire disks.In contrast, models with larger grains and inactive outer disks, show flattopped, or double-peaked profiles in the inner disk regions and highly settled outer regions.
A perfect flat-top distribution would indicate spatially homogeneous diffusion and could easily be fitted by an analytic expression (see Equation 21, Fromang & Nelson 2009).The double hump, on the other hand, cannot be a feature of iso-tropic turbulence and reflects the action of the quasi-periodic VSI motions.
At this point, we can only speculate what the feedback of these dust distributions onto the VSI would be.Lin & Youdin (2017) and Lin (2019) studied the influence that dust back-reaction could have on the VSI and found that this process generally damps the VSI turbulence.For the highly settled cases, with midplane dust-to-gas ratios near unity, one would have to include hydrodynamic back-reaction, as in the work by Schäfer et al. (2020); Schäfer & Johansen (2022).In these scenarios, the presence of VSI would probably be further inhibited by the hydrodynamic feedback of the dust onto the gas.Cooling times would also increase significantly in these regions.The areas above the midplane would be in the collision-limited regime, whereas the midplane could become optically thick (see Section 5).

RADIATIVE TRANSFER POST PROCESSING
We have shown the impact of the dust grain sizes on the strength of the VSI and the morphology of the dust layer in the previous section.Now, we want to determine the visual appearance of the simulated disks in synthetic millimeter-wave observations.Our goal is a qualitative comparison of our results with ALMA observations of edge-on or almost edge-on protoplanetary disks.Specifically, the works of Villenave et al. (2020Villenave et al. ( , 2022Villenave et al. ( , 2023) ) have shown that many protoplanetary disks appear settled in λ = 1.25 mm images obtained with ALMA.Oph 163131 is the most prominent example with a very thin dust disk of height H d,100 au ≈ 0.5 au.Villenave et al. (2022) obtained this result by modeling the appearance of one of the disk's gaps.
For our approach, we create radiation intensity maps of edge-on disks (i = 90 • ) from the dust distributions of the last snapshot of our hydrodynamic simulations with RADMC-3D.For comparison, we also simulate the intensities arising from steady-state dust distributions under the assumption of a fixed diffusivity.In this settlingmixing equilibrium, the vertical dust distribution can be written (Fromang & Nelson 2009).Opacities are calculated for each of the four populations using the standard DSHARP particle properties with the dsharp_opac python package (Birnstiel et al. 2018).We consider a photon package to be fully extinct after being scattered over a length of five optical depths.Our models are axisymmetric and we treat the anisotropic scattering angle for 60 angular sample points.Before running the ray  total dust-to-gas ratio Figure 6.Total dust to gas ratios in VSI simulations restarted after 425 orbits with four passive dust fluids.Each simulation is started with a dust distribution similar to the one derived from the respective DustPy simulations.Snapshots are taken after 150 orbits of evolution.White contours mark the position at which the critical cooling time for the VSI is reached (Equation 3), i.e., VSI is theoretically possible within the white lobes.
tracing algorithm, we use the mctherm task to calculate the dust temperatures from a thermal Monte Carlo simulation.For this, we use 10 7 photon packages.
To mimic the effect of a finite beam size in ALMA observations, we convolve our images with a circular Gaussian beam, which for DSHARP observations had a typical FWHM of 35 mas.We place our disk at a distance of 100 pc to the observer.
We show the resulting images for the VSI simulation with v fr = 100 cm s −1 in Figure 8, v fr = 200 cm s −1 in Figure 9, and for v fr = 400 cm s −1 in Figure 10.The right-hand side of each figure depicts three minor axis cuts through the intensity map at the locations of the vertical lines in the images.The images within each figure are created from disk models with identical particle sizes.As a result of optical depth effects, we find that the models with v fr = 200 cm s −1 (Figure 9) v fr = 400 cm s −1 (Figure 10), have a double-peaked intensity profile in the optically thick regions, marked by the hatched areas in each image.Above the midplane, these models have optical surfaces closer to the central star.Therefore, we observe the hotter inner regions above the midplane and the cooler outer regions in the disk midplane, as illustrated in Figure 11.Doublepeaked profiles have already been observed in synthetic images of a VSI active disk in Blanco et al. (2021).Their work is based on the simulation presented in Flock et al. (2020) and also treats radiative transfer through radia- Radially and time-averaged dust-to-gas ratios in the inner and outer parts of our simulations.The inner regions are VSI active, forming plateau-like dust distributions in all simulations with a cutoff at the edges of the VSI active zones.The outer disk regions appear much more settled in that cases of v fr = 200 cm s −1 and v fr = 400 cm s −1 , in which the outer regions are quiescent.
tive diffusion in combination with ray-tracing from the central star for up to 10 µm dust particles.
The disk model with the smaller particles (v fr = 100 cm s −1 ), is subject to the strongest VSI and the strongest vertical mixing (row a of Figure 8).Therefore, the disk midplane is not as strongly enriched and remains optically thin outside of ∼ 45 au.We are therefore not observing any double-peaked minor axis intensity profiles in these cases.The minor cut intensity profiles in the inner disk match best with the analytic profile with δ = 10 −4 or δ = 10 −3 (rows b and c in Figure 8).In the outer disk, they show almost no settling, since the VSI is still active under the given conditions (more comparable with large diffusivities as in row c in Figure 8).Similar to the conclusions of Dullemond et al. (2022), we confirm that such a disk structure is not consistent with observations of highly settled edge-on disks like Oph 163131.
In our disk model with v fr = 200 cm s −1 , we find a vertically extended and optically thick inner disk with the typical two intensity peaks.However, we can already see the effect of the radially increasing cooling times in the regions beyond 100 au.While the inner, VSI active regions appear the be more consistent with the analytic models of high diffusivity (row d in Figure 9), we can see that the outer regions are most consistent with a low diffusivity of ∼ 10 −4 (row b in Figure 9).This would still not be in agreement with observations of Oph 163131, which find the disk to be highly settled at r ≈ 80 au.
Ramping up the fragmentation threshold further, as in our model with v fr = 400 cm s −1 , results in a highly settled outer dust layer outside of ∼ 60 au, as can be seen in Figure 10.The minor axis cuts illustrate the transition from an optically thick vertical structure in the inner regions to a mostly optically thin profile in the outer regions, which occurs at the outer edge of the VSI active region.For the inner regions, we find a good agreement between the VSI simulation and the analytic model with δ = 10 −3 (row c in Figure 10).Similar to Dullemond et al. (2022), we find that the VSI can still lift up large particles in these inner regions.In contrast, the outer regions are strongly settled, and more consistent with the analytic profile with δ = 10 −4 (row b in Figure 10).At this level of settling, it is unlikely that the outer regions of the VSI simulation could still be distinguished from a fully settled disk (δ = 0), due to the applied beam smearing.Note that in this simulation, we allow dust to flow into the simulation domain with a vertical distribution equal to the initial condition (which assumes δ = 10 −3 ).Any remaining vertical extent of the dust layer in the outer disk therefore likely exists as a result of the boundary condition.

Other Modes of Thermal Relaxation
We assumed the dust to be the only source of cooling in the outer regions of protoplanetary disks.However, molecules like CO, H 2 O, CO 2 , etc, with electromagnetic dipole moments, might also contribute to the cooling of the gas through line emission when gas and dust are thermally decoupled (Woitke et al. 2009;Malygin et al. 2014).In this case, thermal energy must also be transferred from the bulk constituent of the disk, H 2 , to the emitting species via collisions.Cooling the VSI modes could, thus, again become a matter of collision timescales at the very low densities of the outer disk.At low temperatures, emission lines may also become extremely inefficient at cooling the gas at the required rate.Freeze-out of emitting molecules might also reduce the rate of thermal relaxation that can be achieved by emission line cooling.How much material can freeze out and thus be stopped from cooling the H 2 , depends also on the availability of small grains.Cooling of the disk via gas emission lines is, thus, also dependent on the details of the dust population.Future studies should aim to incorporate some treatment of gas cooling via emission lines.Models for this exist (Woitke 2015), but are very complex and currently not feasible for implementation in a hydrodynamic simulation.Furthermore, we have omitted the optically thick regions of protoplanetary disks (R < 10 au) in our simulations.Optically thick in this context does not refer to the bulk optical depth of the disk (τ ∼ Σκ), as discussed in the previous section, but to the optical depth of individual VSI flow structures, which in the inner disk measure only a fraction of the disk scale height in the radial direction(denoted as l in the following).We attempted to simulate these regions in Pfeil & Klahr (2021) by assuming a characteristic diffusion length scale.However, self-consistent modeling requires some treatment of radiative transfer, as in Stoll & Kley (2016) or Flores-Rivera et al. (2020).Our findings nonetheless allow us to make predictions about the effect of dust coagulation on the cooling times in these regions, based on the results obtained here.If radiative diffusion becomes the dominant channel for thermal relaxation, we can write the respective cooling time as where κ R is the Rosseland mean opacity, which is mostly determined by the small grains of density ρ small (Lin & Youdin 2015;Dullemond et al. 2022).If coagulation is increasing the maximum particle size, the density of small particles will be reduced, therefore reducing the diffusion timescale.At the same time, the size distribution-averaged opacity will also be reduced.Therefore, dust coagulation would effectively reduce the diffusion time scale and thus be beneficial for the VSI in the inner disk regions.

Implication for the Vertical Shear Instability
We have shown that the vertical shear instability is highly sensitive to the underlying dust size distribution, which determines the timescale of thermal relaxation.Manger et al. (2021) and Klahr et al. (2023) have shown that the VSI growth rate almost instantaneously drops to almost zero once the critical cooling time threshold is reached.This is also what we observe as a sudden cutoff in the VSI activity at large disk radii.Therefore, the VSI active zones in protoplanetary disks are not extending throughout the entire outer disk.Our simulations predict a VSI dead zone at large radii, which is caused by the reduced efficiency of cooling.
Our simulations omit a treatment of dust backreaction onto the gas.Schäfer et al. (2020) have shown, that if the dust can settle into a thin layer in the disk midplane before the VSI starts to grow, dust feedback can counteract the VSI.Since dust coagulation, settling, and the onset of the VSI, occur on comparable timescales, it is not trivial to predict the outcome of such a situation without a realistic disk simulation that treats all of the aforementioned effects simultaneously.Our results show that if some dust settling and coagulation can occur before the onset of the VSI, the effect of the reduced cooling time would reduce the VSI activity and therefore probably enhance the dampening effect of the dust's dynamic back-reaction onto the gas.6.3.The Need for a Self-consistent Three-dimensional Model and the Limitations of Our Approach Simulations that aim to study the VSI under realistic conditions cannot ignore the implications of an evolved dust population, as presented in our and previous studies (see Fukuhara et al. 2021Fukuhara et al. , 2023)).Measurements of the spectral index in protoplanetary disks (Tazzari et al. 2016;Pérez et al. 2012;Huang et al. 2018;Sierra et al. 2021) and polarization observations (Ohashi & Kataoka 2019) imply that dust coagulation is occurring and that grains in the outer disk can reach sizes of between 0.1-1 mm, similar to the outcome of the DustPy models that our VSI simulations are based on.Note, however, that our studies are no self-consistent representations of protoplanetary disks.The dust size distributions used to calculate the cooling time in our setups are static.In a real disk, they would evolve together with the VSI.Settling and stirring of the dust layer would impact the cooling times.It is unclear if this would lead to some sort of equilibrium situation in which the dust stirring by the VSI can maintain a thick enough dust layer to support the necessary cooling times.Continuous coagulation of grains would counteract the turbulent mixing further.Fukuhara et al. (2023) presented an approach to study this equilibrium by using analytic, yet physically motivated, cooling time profiles.They iterated between VSI simulations and calculations of the resulting steady-state dust distribution from the measured turbulent diffusivity.In that way, they were able to reach a convergent state in which the VSI turbulence creates the necessary diffusivity to maintain the underlying cooling times.Their studies did, however, not consider the effect of the changing diffusivity on the grain size itself through coagulation and fragmentation.This poses an additional uncertainty in their and our studies.We can already see that the measured Mach numbers in our simulations do not always correspond to the α values used in the underlying coagulation models (see Figure 4).Note that M is only part of the generation of turbulent collision speeds (Ormel & Cuzzi 2007).The turbulent spectrum in correlation time space is also required to calculate the acceleration that can be imposed on various particle sizes.Collision speeds can only be obtained from the large scale rms velocity U (L) and the associated length scale L = √ αH, for an ideal Kolmogorov turbulence cascade which causes isotropic turbulent diffusivities (Youdin & Lithwick 2007;Binkert 2023).
If any source of additional turbulence would be present that causes the turbulent diffusivities used in our coagulation models, this would also have an effect on the developing VSI.Even small viscosities of α = 1 × 10 −4 − 1 × 10 −3 are enough to hinder the evolution of the VSI (Barker & Latter 2015).Future studies should try to apply a more realistic, self-consistent prescription of diffusivities in the coagulation model.
In our cooling time calculations, we have also neglected the effects of radial drift.Drift-limited size distributions are characterized by smaller maximum particle sizes and are more top-heavy than fragmentationlimited distributions.This results in longer thermal accommodation timescales and would further inhibit the VSI turbulence.
The effect of the drag force onto the gas was also not considered in our simulations.Schäfer et al. (2020) and Schäfer & Johansen (2022) have shown that backreaction can indeed inhibit the VSI turbulence close to the disk midplane if the dust has time to sediment before the VSI is saturated.Future studies should therefore aim to incorporate more realistic dust dynamics.
In our two-dimensional simulations with dust, we have observed flat top or double-peaked dust-to-gas ratio distributions.This reflects the periodic and non-isotropic nature of the VSI-driven turbulence, which is not accounted for in the coagulation simulations.However, as our simulations are two-dimensional, the prominence of these features might be artificially enhanced, as the φ−dimension is missing as a degree of freedom.Threedimensional simulations (Manger & Klahr 2018;Flock et al. 2020;Manger et al. 2021;Pfeil & Klahr 2021) are needed for the study of the non-linear saturation and fully developed turbulent state of VSI-driven turbulence, before deriving the turbulence properties as diffusivity, correlation times, and energy spectra.
The main conclusions of our study and Fukuhara et al. (2021), however, remain unchanged by all these considerations.Dust coagulation and dynamics are essential components in studies of cooling-time-sensitive instabilities like the VSI.
This highlights the need for a more-self consistent numerical approach.Cooling times have to be constantly recalculated throughout a simulation from the present dust size distributions in order to study such systems.In the inner, optically thick parts of the disk, radiative transfer models have to be employed to study the effect of coagulation on diffusive radiative cooling.

SUMMARY AND CONCLUSIONS
In this work, we studied the effect of evolved dust size distributions on the VSI activity in protoplanetary disks.We conducted hydrodynamic simulations based on five different dust coagulation models for different fragmentation velocities and assumed turbulence strengths, which resulted in maximum particle sizes between ∼ 10 µm and ∼ 0.1 cm.Based on these dust size distributions, we calculated the cooling times for our subsequent hydrodynamic simulations.Our results show a strong dampening effect of dust coagulation on the VSI, as predicted by previous studies (Lin & Youdin 2015;Fukuhara et al. 2021;Pfeil & Klahr 2021;Dullemond et al. 2022;Fukuhara et al. 2023).The reason for this is the collisional decoupling between dust particles and gas molecules that is enhanced if dust coagulation is increasing the maximum particle size.Reduced collision rates inhibit the thermal accommodation of dust and gas and therefore reduce the cooling rate of the gas.
The effect can be strong enough to hinder the development of the VSI, leading to a highly settled dust layer even for moderate fragmentation velocities of v fr ≳ 200 cm s −1 .At the same time, the inner regionsin which the gas and dust components remain well coupled-can maintain some level of VSI turbulence.This finding is consistent with recent observations of highly settled dust layers in protoplanetary disks (Villenave et al. 2020(Villenave et al. , 2022)).Our simulations also show that even a low level of VSI can still significantly alter the vertical distribution of dust, which we can observe in the inner disk regions of our simulations with the largest particles.Synthetic millimeter-images of these VSI active regions are mostly consistent with analytic models that assume large diffusivities of δ ∼ 10 −3 − 10 −2 .At the same time, outer disk regions can appear completely settled in our simulations.We thus report the existence of a VSI dead zone in the outer regions of protoplanetary disks.The existence of the VSI dead zone in the outer regions of protoplanetary disks reconciles recent millimeter-wave observations with models of hydrodynamic turbulence.
Future studies of VSI-active disks should aim to incorporate a more self-consistent treatment of dust coagulation and dynamics.Additionally, cooling via gas emission lines has to be considered to gain a better understanding of the impact of thermal relaxation on the VSI in protoplanetary disks.For this, thermochemical modeling is required to track the amounts and the evolution of relevant species, which in fact also depends on the dust coagulation process.Modeling the optically thick parts of protoplanetary disks and the impact of stellar irradiation furthermore requires radiative transfer modeling.
After applying our methodology to smooth disks in this article, we will extend our studies to disks with substructure in Part II.Specifically, Oph 163131 (Villenave et al. 2020;Wolff et al. 2021;Villenave et al. 2022) and HD 163296 (Dullemond et al. 2018;Rosotti et al. 2020;Doi & Kataoka 2021) have been extensively surveyed with a focus on the dust diffusivities and provide good conditions for comparison with simulations.This means that our two-dimensional cooling time distribution can be described as a function of ρ g and T with a total of eight parameters.For the actual fitting procedure, we use the scipy routine curve_fit.Residuals between the actual cooling time maps and our fitting functions are displayed in Figure 12.With the exception of our model with v fr = 100 cm s −1 and α = 10 −4 , all fits have maximum deviations from the data of < 30 % in the entirety of the simulation domain.The respective model with higher deviations corresponds to a highly settled particle layer with α = 10 −4 , making it distinct from the other models with α = 10 −3 .Thus, some deviation had to be expected for this case.Note that we are generally interested in trends of the VSI activity with the given coagulation parameters, which are well captured by our fits.Uncertainties of 10-60 % in the cooling times thus do not influence the conclusions of our work.

C. DUST ADVECTION AND DIFFUSION IN PLUTO
As a consequence of the sub-Keplerian azimuthal gas velocity, particles in aerodynamic force equilibrium with the gas also have sub-Keplerian terminal velocities.Radial pressure forces, which contribute the to gas' radial force balance, do not significantly act on the dust grains.Therefore, the grains embedded in the gas cannot stay on circular orbital trajectories and spiral inward at a given terminal drift speed.Nakagawa et al. (1986) derived this radial drift velocity as where in our simulations P is the gas pressure, and therefore subject to fluctuations arising from the VSI.We use the zeroth order approximation for small dust-to-gas ratios on the right-hand side.This approximation is robust in the VSI-active regions, where the dust-to-gas ratios are generally smaller than 0.05 in our simulations.A similar derivation can be made for the vertical velocity component.Pressure forces keep the gas on elevated trajectories around the central star (acting against the vertical component of the gravitational force).Again, these forces have a negligible effect on the grains.An expression, equal to Equation C20, can be found for the settling velocity of the grains The PLUTO code already allows for the treatment of passive tracer fluids, which are simply advected with the gas following ∂(ρ g ε) ∂t + ⃗ ∇ • (ε ρ g ⃗ v g ) = 0. (C22) In our case, the advected quantity ε represents a local dust-to-gas ratio ε = ρ d /ρ g .In the short friction time, terminal velocity approximation, the respective dust flow is modified to simulate a dust fluid that is aerodynamically coupled to the gas, i.e., undergoes radial and azimuthal drift, and vertical sedimentation, with terminal velocities given by Equation C20 and Equation C21. to the tracer flux, which is valid for small dust-to-gas ratios (Youdin & Goodman 2005).The tracer flux is then calculated with the upstream dust density based on the above velocity at the respective cell interface ⃗ F drift,interface = ρ d,upstream ⃗ v d−g,interface . (C23) C.1.Test Case for Dust Drift Youdin & Shu (2002), presented an analytic description for the time evolution of a dust fluid with a fixed grain size in a protoplanetary disk due to radial drift.Their prescription was further developed in Birnstiel & Andrews (2014), and their general solution to the advection equation is given by Σ(r, t) = Σ(r 0 , t) v(r 0 )r 0 v(r)r , (C24) with the velocity v(r) and the original location of the characteristic r 0 , defined via This results in a steady state dust distribution, given by where the width of the dust density distribution is related to the pressure bump's width via We conducted a simulation of this setup embedded in a protoplanetary disk with a power law background gas density and four dust fluids with Stokes numbers 10 −3 , 5 × 10 −3 , 10 −2 and 5 × 10 −2 .The test domain spans from 10 au to 30 au with 200 grid cells, and contains a pressure bump of 2 au width at 20 au distance to a solar mass star.Within the pressure bump, where the analytical solution applies, the resulting dust profiles are in excellent agreement with the analytic solutions.As a second test, we set up a two-dimensional, axisymmetric simulation in spherical coordinates.We set up a vertically isothermal disk in hydrostatic equilibrium.Similar to the radial pressure bump, an analytic solution for the settling-mixing equilibrium can be derived which reads where ϑ = tan(z/R).We find a good agreement between the steady state profiles and the theoretically predicted steady state, as can be seen in Figure 15.

Figure 1 .
Figure 1.Workflow from our dust coagulation models with DustPy to the hydrodynamic simulations with PLUTO to the radiative transfer modeling with RADMC-3D.Used methods and tools are shown as dark blue boxes.Input parameters and intermediate results are shown in light blue.The results of our work are schematically displayed as orange boxes.The details of our methodology are laid out in Section 3 (DustPy and cooling times), Section 4 (PLUTO simulations and results), and Section 5 (radiative transfer and synthetic observations with RADMC-3D).

Figure 3 .
Figure 3. Vertical velocities in units of the local speed of sound in our six PLUTO runs after 500 orbital time scales at 50 au.The isothermal run shows a snapshot after only 200 orbits.White contours mark the position at which the critical cooling time for the VSI is reached (Equation3), i.e., VSI is theoretically possible within the white lobes.

Figure 5
depicts the radial dependence of the Mach numbers in our simulations.The lowest turbulence levels of ⟨M⟩ ≈ 8 × 10 −3 are reached in our simulations based on the DustPy model with α = 10 −3 and v fr = 400 cm s −1 , i.e., in the model with the largest par-ticles (a max (50 au) ≈ 0.14 cm).For this simulation, we observe a decrease in turbulence outside of 40 au.At 60 au, turbulent Mach numbers have already decreased by a factor 10 compared to the inner regions.Also our models with v fr = 200 cm s −1 and α = 10 −3 and the model with v fr = 100 cm s −1 and α = 10 −4 show a radially decreasing level of turbulence in the outer disk.

vFigure 4 .
Figure 4. Time evolution of vertical shear instability simulations based on the different dust models.Turbulent Mach numbers are shown as a function of time (radially and vertically Favre-averaged) and as a function of height above the midplane (timeaveraged and radially Favre-averaged).In models with larger particles, cooling times are generally longer, which results in lower growth rates and lower Mach number turbulence.The vertical profiles on the right-hand side change accordingly.Cooling times in models with larger maximum particle size increase more rapidly with height above the midplane, which also cuts off the VSI turbulence.Isothermal models typically have vertically increasing turbulent velocities.The three dashed horizontal lines show the Mach numbers corresponding to the three α values that we assumed for our coagulation models (see Equation10).Note that the conversion between turbulent Mach numbers and diffusivities assumes a perfect Kolmogorov turbulence spectrum (see discussion in Section 6), which is likely not given for the anisotropic VSI turbulence.

Figure 5 .
Figure5.Radial dependency of the turbulent Mach numbers in a polar and time Favre average over 200 orbits in our VSi simulations.VSI simulations based on DustPy models with larger particles have lower levels of turbulence.For our model with the largest particles v fr = 400 cm s −1 and α = 10 −3 , the outer disk, beyond 80 au is completely quiescent.The three dashed horizontal lines show the Mach numbers corresponding to the three α values that we assumed for our coagulation models (see Equation10).Note that such a conversion assumes a perfect Kolmogorov turbulence spectrum (see discussion in Section 6).
× 10 −5 1.98 × 10 −5Note.This value was used as a maximum particle size in our PLUTO simulations with dust.We also show the respective Sauter mean radii and Stokes numbers.
Figure7.Radially and time-averaged dust-to-gas ratios in the inner and outer parts of our simulations.The inner regions are VSI active, forming plateau-like dust distributions in all simulations with a cutoff at the edges of the VSI active zones.The outer disk regions appear much more settled in that cases of v fr = 200 cm s −1 and v fr = 400 cm s −1 , in which the outer regions are quiescent.

Figure 8 .
Figure 8. Upper row a: RADMC-3D intensity maps of our VSI simulation with v fr = 100 cm s −1 and α = 10 −3 , seen edge-on.Rows b, c, and d show intensity maps calculated from analytic dust distribution that assume different diffusivities δ.The grain sizes are identical in all simulations.We convolve the images with a typical ALMA beam with FWHM of 35 mas for a distance of 100 pc shown as a grey circle.Hatched areas mark regions that have optical depth τ ≥ 1. Horizontal hatches correspond to areas for which the τ = 1 surface lies on the far side of the disk.Diagonally hatched regions mark τ = 1 surfaces that lie on the observer's side of the disk.The panels on the right-hand side show minor axis cuts through the images along the vertical lines in the intensity maps.Purple lines in all plots are the minor axis cuts from the VSI simulation (panel a).

Figure 9 .
Figure 9. Upper row a: RADMC-3D intensity maps of our VSI simulation with v fr = 200 cm s −1 and α = 10 −3 , seen edge-on.Rows b, c, and d show intensity maps calculated from analytic dust distribution that assume different diffusivities δ.The grain sizes are identical in all simulations.We convolve the images with a typical ALMA beam with FWHM of 35 mas for a distance of 100 pc shown as a grey circle.Hatched areas mark regions that have optical depth τ ≥ 1. Horizontal hatches correspond to areas for which the τ = 1 surface lies on the far side of the disk.Diagonally hatched regions mark τ = 1 surfaces that lie on the observer's side of the disk.The panels on the right-hand side show minor axis cuts through the images along the vertical lines in the intensity maps.Purple lines in all plots are the minor axis cuts from the VSI simulation (panel a).

Figure 10 .
Figure10.Upper row a: RADMC-3D intensity maps of our VSI simulation with v fr = 400 cm s −1 and α = 10 −3 , seen edge-on.Rows b, c, and d show intensity maps calculated from analytic dust distribution that assume different diffusivities δ.The grain sizes are identical in all simulations.We convolve the images with a typical ALMA beam with FWHM of 35 mas for a distance of 100 pc shown as a grey circle.Hatched areas mark regions that have optical depth τ ≥ 1. Horizontal hatches correspond to areas for which the τ = 1 surface lies on the far side of the disk.Diagonally hatched regions mark τ = 1 surfaces that lie on the observer's side of the disk.The panels on the right-hand side show minor axis cuts through the images along the vertical lines in the intensity maps.Purple lines in all plots are the minor axis cuts from the VSI simulation (panel a).

Figure 11 .
Figure11.Origin of the double-peaked intensity profiles in Figure9and Figure10.The τ = 1 surfaces for the layers above the midplane lie closer to the central star due to the lower densities.The respectively higher temperatures lead to higher intensities above the midplane.Here shown are the τ = 1 surfaces for row c in Figure9.

,Figure 14 .
Figure 12.Cooling time maps in our 5 VSI PLUTO simulations on the left-hand side, and the residuals between the fit used in our simulations, and the actual cooling time distributions derived from DustPy simulations.

Table 1 .
Dust coagulation parameters of our five DustPy simulations and the respective maximum particle size measured at 50 au in the DustPy simulation.