Razor-thin dust layers in protoplanetary disks: Limits on the vertical shear instability

Context: Recent observations with the Atacama Large Millimeter Array (ALMA) have shown that the large dust aggregates observed at millimeter wavelengths settle to the midplane into a remarkably thin layer. Aims: We intend to find out if the geometric thinness of these layers is evidence against the vertical shear instability (VSI) operating in these disks. Methods: We performed hydrodynamic simulations of a protoplanetary disk with a locally isothermal equation of state, and let the VSI fully develop. We sprinkled dust particles and followed their motion as they got stirred up by the VSI. We determined for which grain size the layer becomes geometrically thin enough to be consistent with ALMA observations. We then verified if, with these grain sizes, it is still possible to generate a moderately optically thick layer at millimeter wavelengths, as observations appear to indicate. Results: We found that even very large dust aggregates with Stokes numbers close to unity get stirred up to relatively large heights above the midplane by the VSI, which is in conflict with the observed geometric thinness. For grains so large that the Stokes number exceeds unity, the layer can be made to remain thin, but we show that it is hard to make dust layers optically thick at ALMA wavelengths (e.g., tau(1.3mm)>=1) with such large dust aggregates. Conclusions: We conclude that protoplanetary disks with geometrically thin midplane dust layers cannot be VSI unstable, at least not down to the disk midplane. Explanations for the inhibition of the VSI include a reduced dust-to-gas ratio of the small dust grains that are responsible for the radiative cooling of the disk. A reduction of small grains by a factor of between 10 and 100 is sufficient to quench the VSI. Such a reduction is plausible in dust growth models, and still consistent with observations at optical and infrared wavelengths.


Introduction
According to canonical theory, the evolution of protoplanetary disks is thought to be driven by a combination of viscous evolution and photoevaporation (e.g., Clarke et al. 2001).The viscosity in such disks is thought to be caused by turbulence produced by the magnetorotational instability (MRI).While the MRI is inhibited in very dense regions of the disk (the so-called dead zones), it may still be operational in the hot inner regions (r ≪ 1 au), in the irradiated surface layers, and in the weakly ionized outer regions (r ≫ 1 au) of the disk (e.g., Dzyurkevich et al. 2013).
However, in recent years the turbulent viscous disk theory for the outer regions of protoplanetary disks has been called into question.Using radiative transfer modeling of the Atacama Large Millimeter Array (ALMA) image of HL Tau, Pinte et al. (2016) infer that the turbulence in that disk must be weak (α ≲ 10 −3 ).Direct measurements of the turbulent line width in CO 2-1 with ALMA show mostly upper limits to the turbulent velocities of ≲10% of the local sound speed (e.g., Flaherty et al. 2020).While these velocity upper limits are still consistent with turbulent α values of up to 10 −2 , and marginally consistent with MRI-turbulent disks (Flock et al. 2017), these and other measurements have stimulated the exploration of the possible consequences of the absence of MRI turbulence in protoplanetary disks and, consequently, the possibility that these disks may be much less turbulent than previously thought.
The implications of very low turbulent α values for protoplanetary disks are numerous.For instance, Bae et al. (2017) show that in low-α disks a single planet can produce multiple rings.Indeed, Zhang et al. (2018) demonstrate that a protoplanetary disk model with very low α and a single embedded planet can reproduce the observed many-ringed structure of the disk around AS 209 remarkably well (see, however Ziampras et al. 2020).Low turbulent velocities also have strong implications for dust growth, gap formation, planet migration, and many other A105, page 1 of 17 Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0),which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.This article is published in open access under the Subscribe-to-Open model.Subscribe to A&A to support open access publication.A&A 668, A105 (2022) things.The reason why, until not long ago, turbulent α values of ≲10 −3 were not seriously considered by scientists in the field is that most stars with protoplanetary disks are observed to undergo substantial gas accretion.This requires α values in excess of about 10 −3 (e.g.Hartmann et al. 1998).However, wind-driven accretion may provide a solution to this dilemma (e.g., Ferreira & Pelletier 1995;Tabone et al. 2021;Martel & Lesur 2022).
In the absence of the MRI, a protoplanetary disk can be prone to other, non-magnetohydrodynamic instabilities that cause turbulence or turbulence-like velocity fluctuations (e.g., Pfeil & Klahr 2019).Of particular importance in the outer regions of the protoplanetary disk is the vertical shear instability (VSI, e.g.Nelson et al. 2013;Stoll & Kley 2014, 2016;Flores-Rivera et al. 2020).This instability, when fully developed, produces upward and downward vertical streams of gas that slowly oscillate.These oscillations form a radially propagating wave (Svanberg et al. 2022).When viewed as a kind of turbulence, it is highly anisotropic, with turbulent "eddies" being radially narrow, but vertically extended sheets of gas moving either up or down.This "turbulence" is only weakly effective as a replacement of MRI turbulence for the radial transport of angular momentum, with values on the order of α VSI,radial ∼ 10 −4 (Stoll & Kley 2014).However, due to the strong upward and downward motions of the gas, crossing the midplane, with velocities in the range of 5-20% of the isothermal sound speed, the effect of the VSI on the dust population of the disk is very pronounced (Stoll & Kley 2016;Flock et al. 2017;Lin 2019;Lehmann & Lin 2022).Even large dust particles can be stirred up to high elevations above the midplane.
From the observational side, however, there is now increasing evidence that many, if not most, protoplanetary disks contain a layer of large dust aggregates at the midplane that contains a substantial amount of dust mass and is geometrically extremely flat, i.e., having a very small scale height.The first evidence came from the ALMA image of the disk around HL Tau, where radiative transfer modeling put an upper limit on the vertical scale height of the dust layer of 1 au at a radius of 100 au (Pinte et al. 2016).The high resolution ALMA images of the DSHARP campaign (Andrews et al. 2018) also suggest very flat geometries of the dust layers seen at λ = 1.3 mm wavelengths.Detailed radiative transfer analysis of the DSHARP observations of HD 163296 shows that the dust in the inner ring of that source (at r ≃ 67 au) appears to be vertically extended almost to the gas pressure scale height, but the outer ring (at r ≃ 100 au) appears to be less than 10% of the gas pressure scale height, i.e., highly settled (Doi & Kataoka 2021).
To get better constraints on the vertical extent (geometric thickness) of the midplane dust layers of protoplanetary disks, dedicated observing campaigns with ALMA for nearly-edge-on disks are required.The first such campaign already yielded indications of strong settling of large grains (Villenave et al. 2020).But when the disk of Oph 163131 was reobserved with ALMA by Villenave et al. (2022), the vertical scale height of the dust layer could be constrained to be less than 0.5 au at a radius of 100 au, which is about 7% of the gas pressure scale height at that radius.From these observations, and under some assumptions of the grain size, these authors derive an upper limit of α ≲ 10 −5 on the turbulence.
However, if the VSI is operational in these outer disk regions, one might expect that the dust layer should be much more vertically extended, due to the high efficiency of the vertical dust stirring of the VSI.The purpose of this paper is to quantify this.
The grain sizes are not perfectly known, nor is the gas disk density.We address the question whether the grains could be so large that they remain in a thin layer in spite of the VSI.And if not, what could be the reason that the VSI is not operational in this disk.
The paper is structured as follows: We start with an analysis of the stirring-up of particles in Sect. 2. In Sect.3, we explain why St ≫ 1 particles are not a probable explanation for the thin dust layers.We propose a natural explanation for the absence of VSI in protoplanetary disks in Sect.4, and we finish with a discussion and conclusions.

Stirring up of large dust aggregates by the VSI
In this paper we wish to find out if the presence of the VSI in the outer regions of a protoplanetary disk, such as the one around Oph 163131, would inevitably lead to the big-grain dust layer observed with ALMA to be more geometrically extended than observed.This would then be clear evidence that the VSI does not operate in that disk.
The effect of the VSI on dust particles in the disk was studied by several papers (e.g., Lorén-Aguilar & Bate 2015;Stoll & Kley 2016;Flock et al. 2017;Lehmann & Lin 2022).The models we present in this section are not fundamentally different from those earlier papers.However, we explore the parameters and compare the results to the observational constraints.

Conveyor-belt estimate of vertical mixing efficiency of the VSI
A simple estimate of the height above the midplane that a dust aggregate can be lifted by the VSI would be the following.
Assume that the VSI consists of long-lived vertical upward and downward moving slabs of gas, acting as vertical conveyor belts for the dust.As a dust aggregate gets dragged upward, the vertical component of gravity increases linearly with z, leading to a vertical settling of the aggregate with respect to the upward moving gas.The maximum height that the dust aggregate can reach is the height z at which the vertical settling speed equals minus the vertical gas speed.The settling speed of a particle with Stokes number St ≪ 1 at a height z above the midplane is By setting v sett + v z,VSI = 0, with v z,VSI the typical vertical gas velocities of the VSI, we obtain the maximum elevation above the midplane that an aggregate can obtain: where h p is the pressure scale height of the gas (see Appendix A) and c s is the isothermal sound speed of the gas.For typical vertical velocities of the VSI of |v z,VSI | ∼ 0.1 c s , a dust aggregate of St = 0.1 can be stirred up, according to this simple estimate, to about one gas pressure scale height.If Eq. ( 2) gives values larger than h p , the estimate is no longer accurate, and we limit it to h p for convenience.
In practice the mean elevation ⟨z 2 ⟩ of the dust aggregate will be smaller than this value, because the VSI motions are not stationary (see Fig. 3).But this estimate does explain why the VSI can stir even large dust aggregates (St ≃ 1) very far away from the midplane.
The conveyor-belt estimate can be compared to the more traditional settling-mixing equilibrium (see Appendix H).

Particle motion model
A more accurate estimation of how dust aggregates are stirred up from the midplane by the VSI is to compute their detailed motion within a hydrodynamic model of the VSI.We employ the PLUTO code (Mignone et al. 2007) for this.The setup of the disk follows the fiducial model of Appendix A.
We assume the disk to be perfectly locally isothermal and inviscid, which we expect to maximize the VSI activity.Given that the VSI establishes itself primarily in the radial and vertical coordinates, we model it in 2D using spherical coordinates r and θ.The radial coordinate r has 882 grid points logarithmically spaced between 0.2 r 0 and 5 r 0 , where r 0 = 100 au is the reference radius.The vertical coordinate θ (where θ = π/2 is the equatorial plane) has 160 grid points linearly spaced between π/2 − 0.3 and π/2 + 0.3, which corresponds to 20 cells per scale height at r = r 0 .This is enough to resolve the large-scale structure of the VSI (Manger et al. 2020).At r = r 0 the range in θ corresponds to ±4 h p , dropping to ±2.7 h p at r = 5 r 0 .The temperature is fixed in time, and depends on radius r as T ∝ r q with q = −1/2.It is chosen such that at r = r 0 the gas pressure scale height is h p (r 0 ) = 0.0732 r 0 .
We compute the gas dynamics without accounting for dust dynamics.Once the VSI is fully developed, after about 300 orbits at r = r 0 , we extract 200 time snapshots, 0.1 orbits apart in time.The vertical gas motions in the first of these snapshots is shown in Fig. 1 for the entire radial and vertical range of the model.The same is again plotted in Fig. 2 in natural (linear cylindrical) coordinates, which gives a better view of the proportions.In Fig. 3 the vertical gas velocity at the location r = r 0 and z = 0 is shown as a function of time, to show how the VSI motions oscillate with a period of a few local orbits.
For these 20 orbits we now follow, as a post-processing step, the motion of N = 2000 large dust particles that have been randomly placed between 0.75 r 0 and 1.6 r 0 radially, and between −0.001 r 0 and +0.001 r 0 vertically.The particle velocities are initialized as being equal to the local Kepler velocity.The particles all have the same St 0 , meaning that if they were placed at r = r 0 and z = 0, they would have Stokes number St = St 0 .At each time step, for each particle, we recompute St based on the local conditions, consistent with keeping the grain size constant.The equations of motion of the particles include the force of gravity as well as the friction with the gas.We implement the numerical integration of these equations in a Python program.The particles do not have dynamical feedback onto the gas, allowing the gas hydrodynamics to be precomputed, and the dust particle dynamics to be computed in post-processing mode.
A105, page 3 of 17 Given that the particles, in spite their comparatively large size (of the order of ∼millimeter), are much smaller than the mean free path of the gas molecules, the friction force is the simple Epstein drag law.At each time step the local gas temperature, density and velocity are linearly interpolated in time and space from the precalculated 200 snapshots from the hydrodynamic simulation, using the RegularGridInterpolator function of the SciPy library.
In Fig. 4, the results of the model are shown for St 0 = 0.01, St 0 = 0.1, St 0 = 1, and St 0 = 10, after 2.5 orbits at radius r 0 .This short time is enough to achieve approximately the typical heights above and below the midplane that the particles acquire, and this does not change much in time after that.The conveyor-belt estimate of the height z max above the midplane that the particles are stirred (Eq.( 2)) appears to be a reasonably good estimate, as can be seen by comparing the vertical locations of the particles with the light-blue dashed lines in the figure.
It is evident that the St 0 = 0.01 particles are stirred up to one gas pressure scale height.This is entirely due to the VSI, and not due to any α-diffusion, which is not included in the model.For St 0 = 0.1, which is typically the highest Stokes number expected in dust coagulation models in the outer disk regions (Drazkowska et al. 2021), the dust particles are still stirred up to a substantial fraction of the gas pressure scale height.Even for St 0 = 1 the particles get up to a height z/R ∼ 0.005, which is about 7% of the gas pressure scale height, which is marginally consistent with the upper limit obtained for Oph 163131.If we go to St 0 = 10, the particles remain close to the midplane, producing a thin layer well within the vertical geometric thinness of the observed dust layer of Oph 163131.However, as is shown in Sect.3, it is unlikely that the particles in this dust layer have St 0 ≳ 1.
To see this more quantitatively, we compute the root-meansquare of the z/r ratio of the particles ⟨(z/r)2 ⟩ as a function of time since insertion.This is a measure of the vertical extent of the dust layer.The results are shown in Fig. 5-left.As can be seen, for St 0 ≥ 0.1 the particles quickly reach a steady-state vertical extent, which is smaller for larger values of St 0 .
For St 0 = 1 the radial drift of the particles becomes so fast that after 17 orbits the first particles leave the grid of the model at the inner edge.The calculation is then halted.
In general it should be noted that the rapid radial drift of large dust particles is a long-standing problem in the interpretation of millimeter wave observations of protoplanetary disks (Birnstiel et al. 2009).One explanation could be that the gas disks are so massive, that even millimeter particles in the outer regions of protoplanetary disks do not drift at excessive speeds (Powell et al. 2017).Another explanation is that radial drift of these particles is inhibited by dust trapping in one or more local pressure maxima (Pinilla et al. 2012).This would imply that the dust we observe with ALMA in the outer regions of protoplanetary disks (r ≳ 20 au) is either trapped in vortices or in rings.Both features are indeed observed with ALMA in numerous disks (e.g., van der Marel et al. 2013;Dong et al. 2018;ALMA-Partnership 2015;Huang et al. 2018) and thus lend support to this picture.This means that the very flat dust midplane layers, if they consist of a series of concentric rings, could very well be made up of large dust particles, without experiencing the strong radial drift that the particles in our model undergo.
In principle this means that our models should be repeated for the case of disks with radial pressure bumps.However, since the origin of these pressure bumps is not yet clear, this would introduce a series of new and unconstrained model parameters.Also, a variety of additional phenomena could occur in these traps (e.g., Carrera et al. 2021;Lehmann & Lin 2022).So for this paper we limit ourselves to disks without pressure traps.

Dynamics of dust modeled as a fluid
The dust motion can also be modeled directly within the hydrodynamics model.For this we employ the Fargo3D code1 (Benitez-Llambay & Masset 2016), which has dust dynamics built in (Krapp & Benitez-Llambay 2020).Like before, we assume a locally isothermal equation of state, maximizing the VSI activity.The dust is treated as a pressureless fluid, which feels friction with the gas.In the standard setup of Fargo3D, the gas feels the opposite force from the dust.However, to make the comparison with the results of Sect.2.2, we switch this feedback off.It is known that for high metallicity Z the VSI can be hampered simply by the mass of the dust (Schäfer et al. 2020;Lehmann & Lin 2022), which would be one possible explanation for the razor-thin dust disks seen in ALMA.But in this section we assume that this effect is not taking place.The background viscosity is set to α = 10 −6 .
The results for the case of St 0 = 0.1 are shown in Fig. 6.The dust has been allowed to settle from the very beginning of the simulation over the entire modeling time.Throughout this time frame, the pattern of the dust remains corrugated.It is very comparable to the results of Sect.2.2.The main difference is that in the dust-fluid approach the vertical width of the corrugated dust "layer" is thicker than in the particle approach of Sect.2.2.This is due to numerical diffussivity.
Next, we put the result of this model into the RADMC-3D 2 radiative transfer code and compute the images at an inclination of 84 • , at a wavelength of λ = 1300 µm.The big grains were assumed to have a radius of 100 µm, and we used the corresponding opacity for them (see Appendix B).The results are shown in Fig. 7.The corrugated geometry of the dust "layer" is clearly seen.To stress the effect this has on identifying any potential radial gaps in the dust layer, we artificially added a gap between 87 au and 98 au and a slight reduction of the density between 3.5 au and 60 au according to the model of Villenave et al. (2022) for Oph 163131.This was done a-posteriori: the big-grain dust density from the hydrodynamic model of Fargo3D was multiplied by a radial function that reduces the density by a factor of 0.1 between 87 and 98 au and by 0.5 inward of 60 au.After that, it was inserted into the RADMC-3D code.As seen in Fig. 7, these features are not recognizable due to the strong vertical waves.These images are not convolved with the ALMA beam, as they merely serve as an illustration.For the case of Oph 163131, Villenave et al. (2022) show that the spatial resolution of ALMA easily suffices to rule out that the dust layer is as strongly corrugated as in Fig. 7.
To highlight the difference between the VSI model and an equivalent model with a flat midplane layer, we show in Fig. 8 the comparison between these cases, for several inclinations.Again, no beam convolution is applied.Although the disk around Oph 163131 is used as a basis for these models, they are not meant to directly fit Oph 163131, but instead to illustrate the typical protoplanetary disk case (hence the different inclinations shown).It is clearly seen that at high inclinations, the VSI models look very different from the flat models, at scales easily resolvable with ALMA for objects at typical distances of about 100 pc.Also shown is the case where the big dust grains are vertically smeared out in a Gaussian layer with a vertical thickness half that of the gas (h big = 0.5 h p ).This mimicks the case when the vertical dust transport by the VSI would be treated as a vertical turbulent mixing instead of an actual advective transport.This case looks also substantially different from the VSI case.But it will depend on the distance of the object and the ALMA baselines whether they can be distinguished.The differences become less clear at lower inclinations, because the models only differ in vertical direction and are the same in radial direction.
A caveat of these synthetic images of the VSI-stirred dust disk is that we have inserted only a single grain size.If we assume that the large grains follow a size distribution of a certain width, then the strong wiggles seen in the image get smeared out.The degree of smearing-out depends on the width of the size distribution.But it would not affect the conclusions of this paper.
This simulation confirms the results of Sect.2.2 that even particles with a rather high Stokes number of St 0 = 0.1 are stirred up to a substantial fraction of the gas pressure scale height, easily measureable with ALMA, and clearly in contrast with for instance the ALMA observations of Oph 163131.
As in the models of Sect.2.2, the corrugated pattern of the dust layer is not static.It follows the time-dependent variations of the VSI velocity profile, where upward gas motions turn into downward motions and vice versa over time scales of a few local orbits.The model also confirms that the dust is not "vertically mixed" as in the simple vertical mixing-settling model of Dubrulle et al. (1995), Dullemond & Dominik (2004b) and Fromang & Nelson (2009).Instead, the corrugated structure of the dust is maintained, and the vertical extent is better described by the "conveyor-belt model" of Sect.2.1.

Conclusion of this section
In this section we have shown that with a VSI operating in the disk, even particles with a Stokes number close to unity get stirred up to high elevations above the midplane, of the order of the gas pressure scale height.This is in conflict with ALMA observations of several protoplanetary disks, most strikingly the disk around Oph 163131 (Villenave et al. 2022).
However, for St 0 ≫ 1 the midplane dust layer indeed becomes very geometrically thin, even in a disk in which the VSI is operating.So we need to rule out that these particles could have St 0 ≫ 1, which is the topic of Sect. 3.
Once we ruled it out, we had to investigate how the VSI could be suppressed in the outer regions of protoplanetary disks.This is explored in Sect. 4.

The case against the midplane dust layer consisting of St ≫ 1 particles
As Sect. 2 showed, the geometric thinness of the midplane dust aggretate layers in protoplanetary disks can most easily be explained by particles that have St ≳ 1, since they remain in a thin layer in spite of a possible VSI operating in the background.However, the dust rings seen in ALMA observations at λ = 1.3 mm tend to have optical depths larger than about 0.3 at that wavelength (Dullemond et al. 2018).For Oph 163131 in particular, Villenave et al. (2022) find with their radiative transfer modeling that the midplane dust layer is partially optically thick.We demonstrate in this section that having both St ≳ 1 and τ 1.3 mm ≳ 0.3 requires a vertically integrated dust-to-gas ratio of at least Z ≳ 0.08, but likely Z ≳ 0.16.This value increases linearly with increasing St and τ 1.3 mm .To arrive at this, we start with the dust opacity model described in Appendix B. In Fig. 9 the λ = 1.3 mm opacity as a function of grain size for this dust model is shown.Clearly the opacity κ 1.3 mm is a strong function of the grain size a.It has a maximum value of 44.6 cm 2 g −1 at a grain size of a = 0.28 mm.For a → 0 the asymptotic value is κ 1.3mm → 1.65 cm 2 g −1 .For a → ∞ we can express κ λ in terms of the geometric opacity κ geom defined as the geometric cross section πa 2 divided by the grain mass: where ρ s is the mean material density of the dust aggregate, and Q λ is the ratio of the opacity to the geometric opacity (van de Hulst 1957de Hulst , 1981)).For a ≳ 1 mm, Eq. ( 3) provides a good fit to the real opacity for a constant Q λ = 1.6.We note, however, that this equation can be used for all values of a, in which case Q λ will depend on a and drop well below unity for a ≪ λ/2π.If we define the surface density of the big dust grains as Σ big , then the optical depth of the big dust grain layer becomes where a is the radius of the big dust grains.
Next, let us compute, for the big dust grains, the Stokes number St.For the outer regions of the protoplanetary disk we can assume that we are firmly in the Epstein friction regime, so that we can write (Birnstiel et al. 2010): where Σ g is the surface density of the gas.We can combine Eqs. ( 4), ( 5) and eliminate ρ s a, to obtain A105, page 7 of 17 A&A 668, A105 (2022) In the dust opacity model of Appendix B, for λ = 1.3 mm, the maximum value of Q λ is 3.0, which is only reached in a narrow range of grain radii (see Fig. 9).For most values of a, Q λ ≲ 1.6.This means that if both St big ≳ 1 and τ λ,big ≳ 0.3, then Eq. ( 6) shows that Z big ≳ 0.08 • • • 0.16, i.e., the "metallicity" must be extremely high.The question is then: is this a realistic scenario?Can the geometrically thin, but optically marginally thick (τ ≳ 0.3) dust rings seen in many protoplanetary disks (and most strikingly seen in Oph 163131) be rings of dust particles with St big ≳ 1 and Z big ≳ 0.16, or even St big ≫ 1 and Z big ≫ 1, with a dynamics similar to the rings of Saturn?This scenario is completely different from the standard picture of dust dynamics in protoplanetary disks, which assume St big ≪ 1 and Z big ≪ 1.
Although we cannot rule it out, we consider this scenario unlikely.The conditions derived in this section are the minimal conditions required.To be more comfortably within the limits, one would need Z big /St big ≫ 0.16, leading, for St big ≳ 1, to very large values of Z big .
Using this constraint, we are then forced to consider mechanisms quenching the VSI entirely in order to understand the geometrical thinness of the dust rings.One way would be to load so much mass worth of dust in this midplane layer, that the gas is no longer able to lift it up from the midplane.As was shown by Lin (2019), Schäfer et al. (2020) and Lehmann & Lin (2022), the VSI is suppressed if the vertically integrated dust-to-gas ratio ("metallicity") of the big grains exceeds about Z ≳ 0.02 • • • 0.05.Another way is to increase the cooling time scale, which is a natural consequence of grain growth (Fukuhara et al. 2021).

Importance of the cooling efficiency for the VSI
The VSI operates in disks which have, to good approximation, a locally isothermal equation of state.That is, at any given position (r cyl , z, ϕ) in the disk the temperature of the gas T g (r cyl , z, ϕ) is fixed and does not vary in time.The justification for this is that in the outer regions of protoplanetary disks, the thermal household is determined by a balance between irradiation from the central star and thermal radiative cooling by the dust.The radiative cooling time scale t rad cool for any perturbation of this equilibrium is short compared to the orbital time scale.
However, the cooling time is not completely negligibly small compared to the orbital time scale.As we shall show, only a moderate amount of dust coagulation is enough to increase the cooling time (or more accurately, the relaxation time; see Appendix D) beyond the limit where the VSI is stopped.
It was shown by Lin & Youdin (2015) that if the thermal relaxation is not fast enough, the vertical entropy gradient acts as a strongly stabilizing force against the VSI.They derive the following upper limit on the radiative relaxation time t relax : where q is the powerlaw index of the midplane temperature profile of the disk T mid ∝ r q , and γ is the usual adiabatic index of the gas, which for the outer disk regions is γ = 5/3 because the rotational and vibrational modes of H 2 are not excited at those temperatures.For our fiducial disk model (Appendix A) we have q = −1/2, and at r = r 0 = 100 au we have h p /r = 0.0732.So we obtain t relax < 0.055/Ω K as the upper limit on the thermal relaxation time for VSI to be operational.
Where in the protoplanetary disk this condition is met, and where not, was, among other things, explored by Pfeil & Klahr (2019).They found that the VSI is typically operational for radii r cyl ≳ 10 au, which are the regions of protoplanetary disks that have been resolved with ALMA, and where these geometrically thin dust layers are detected.(Fukuhara et al. 2021) explore how dust evolution can change this, and they found that the coagulation of dust grains can increase the relaxation time scale and act against the VSI.A similar conclusion for the Zombie Vortex Instability was found by Barranco et al. (2018).

Gas cooling via small dust grain emission
In this section we revisit the question of the cooling efficiency, and estimate t relax in a simplified, yet robust way, including realistic dust opacities and the effect of dust depletion due to coagulation.We mimick the effect of coagulation by a simple conversion factor X ∈ [0, 1] that says that a fraction X of the small grains has been converted into big grains that are not participating in the radiative cooling of the gas (these are probably the grains we observe with ALMA), while only a fraction (1 − X) of the small grains remain to radiatively cool the gas.In essence, we make the simplifying assumption that the dust consists of only two components: small submicron dust grains that are wellmixed with the gas, and are solely responsible for the radiative cooling, and big millimeter-size grains that tend to settle to the midplane unless they are stirred up by the gas.
In the outer regions of a protoplanetary disk, the gas near the disk midplane is cold: T mid ≲ 70 K.This means that the gas has very few emission lines, and no continuum, by which it can radiatively cool: typically only the rotational transitions of CO and its isotopologs, and maybe a few more complex molecules.Effectively this means that the gas is unable to radiatively cool by itself.It can only cool by transmitting its thermal energy to the available dust grains in the gas, which then can radiate away this energy.
In the midplane regions of the disk, the thermal coupling of gas and dust through collisions of gas molecules with the dust particles, is relatively efficient, though not perfect.The gas-dust thermal coupling time scale is estimated in Appendix F, but first we assume that the gas and dust thermally equilibrate fast enough that we can set T g = T small , i.e., the gas temperature equals the small-grain dust temperature.
We ignore the effect of the large-dust-aggregates midplane layer, and focus only on the gas and the small dust grains floating in the gas.We assume that these small dust grains are well-mixed with the gas in vertical direction, so that the dust-to-gas ratio for these small dust grains is vertically constant.
Under these conditions, the fastest cooling happens in the optically thin regime, because the dust opacity is independent of the dust density (the amount of dust per unit volume of the disk).
The rate of thermal emission of the small dust grains per unit volume of the disk is: where ρ d is the volume mass density of the small dust grains, κ abs ν,small is their absorption opacity as a function of frequency ν, and B ν (T small ) is the Planck function at the dust temperature T small .It is convenient to express this in terms of the Planck mean A105, page 8 of 17 C. P. Dullemond et al.: Razor-thin dust layers: Limits on the vertical shear instability opacity κ P (T small ) defined as with σ SB the Stefan-Boltzmann constant.We can then express q cool,small as q cool,small = 4ρ small κ P (T small ) σ SB T 4 small . (10) In Appendix C we give a convenient approximate expression for κ P (T small ).
The thermal energy in the dust per unit volume of the disk is: with c V,small the specific thermal heat capacity of the dust, c V,small ≲ 10 7 erg g −1 K −1 (Draine & Li 2001).The thermal energy in the gas per unit volume of the disk is: with the specific thermal heat capacity of the gas given by where k B is the Boltzmann constant, µ ≃ 2.3 is the mean molecular weight of the gas in units of the atomic unit mass m u , and γ is the ratio of specific heats.The total thermal energy density is the sum of the two e th = e th,g + e th,small .
For a small-grain dust-to-gas ratio smaller than or equal to 0.01 we can safely approximate this as e th = e th,g .The optically thin radiative cooling time is then assuming T small = T g = T mid .In the optically thin limit this then becomes However, what we need for the analysis of the VSI is the relaxation time, which is shown in Appendix D, to be where the 1.7 is valid for the opacity model of Eq. (C.3) of Appendix C. So far we have not included optical depth effects, and have therefore considered the most VSI-friendly scenario.Optical depth effects can only increase the relaxation time, not shorten  19)) of the disk in units of the Kepler time Ω −1 K for three small-grain dust-to-gas ratios: Normal (Z small = 10 −2 ), depleted by a factor of 10 −1 (Z small = 10 −3 ) and depleted by a factor of 10 −2 (Z small = 10 −4 ).Dotted lines: optically thin approximation (Eq.( 17)).Solid lines: including optical depth effects.In dashed black: upper limit to the relaxation time (Eq.( 7)) for which VSI is operational is shown.The disk and stellar parameters are those of the star Oph 163131 (the fiducial model of Appendix A).
it.We are primarily interested in the regions that are spatially resolvable with ALMA, meaning we are interested in r ≳ 10 au.The optical depth of the disk to its own radiation is moderate to low in these outer regions.Optical depth effects are therefore not expected to play a large role in these regions.But it is not a major effort to include them.In Appendix E we discuss the relaxation time scale in the optically thick regime, and write it as t rad cool,thick given by Eq. (E.1).
Finally, we have to account for the time it takes to transfer heat between the gas and the dust, t dg .This will play a big role for disks around bright stars such as Herbig Ae/Be stars, where it will be the limiting factor of the radiative the cooling.In Appendix F we give an expression for t dg .
We estimate the combined cooling time scale to be the sum of all three time scales 3 : which gives a smooth transition between regions, and ensures that the limiting factor determines the actual relaxation time.In Fig. 10, this relaxation time is shown for the fiducial disk model of Appendix A, for small-grain dust-to-gas ratios of Z small = 10 −2 (no depletion, i.e., X = 0), for Z small = 10 −3 (a factor of 10 depletion of small dust grains, i.e., X = 0.9), and for Z small = 10 −4 (a factor of 100 depletion of small dust grains, i.e., X = 0.99).The dotted lines represent t rad relax,thin + t dg , i.e., without optical depth effects.
In Fig. 11, the same is shown for a 10 times lower disk mass, both in dust and in gas.Because of the lower optical depth, the solid curves are now closer to the optically thin estimates.
One can see that for both the fiducial model and for the 10× lower mass disk, Ω K t relax is well below unity, justifying the locally isothermal appoximation for most applications.However, for the VSI to be operational, the relaxation time has to be below the limit given in Eq. ( 7 For the fiducial model, with small-grain dust-to-gas ratio of 10 −2 , the thermal relaxation time scale is everywhere below this limit, meaning that the disk is prone to the VSI everywhere.However, if dust coagulation converts, say, 90% of the small grains into large dust aggregates (a depletion of 10 −1 , or X = 0.9), leading to a small-grain dust-to-gas ratio of Z small = 10 −2 (1 − X) = 10 −3 , then the Ω K t relax is above the threshold value for r ≳ 50 au.If coagulation converts 99% of the small grains (a depletion of 10 −2 , or X = 0.99), then the curve is everywhere well above the threshold, and the entire disk is VSI-stable.
If we redo our analysis for a brighter star, say a Herbig Ae star, then the disk will be warmer due to the stronger irradiation.This will lower the cooling times and thus make the disk more susceptible to the VSI.We show the resulting cooling time scales for a Herbig Ae star of M = 2.4 M ⊙ and L = 50 L ⊙ , with otherwise the same disk parameters, in Fig. 12.Indeed, the cooling time for a normal dust-to-gas ratio is substantially shorter.A depletion of small grains of a factor of 10 is not sufficient, but a factor of 100 will, again, make the disk stable against the VSI in most of the disk.
A depletion of small dust grain by a factor of 10 or even 100 due to coagulation is not extreme.The fact that most protoplanetary disks look "fat" (geometrically vertically extended) in optical and near-infrared observations is not evidence of a lack of coagulation.In Appendix G we quantify this by computing the optical appearance of our fiducial disk at a wavelength of λ = 0.8 µm for various degrees of small-grain depletion.These images show that the typical appearance of the disk, with its two bright layers separated by a dark lane, is retained even at large degrees of depletion.The required amount of dust coagulation to inhibit the VSI is therefore within the observational constraints.

Earlier work on the stirring of large grains by the VSI
The extreme effectiveness of the VSI to stir up even large dust aggregates to high elevations above the midplane is not a new result, and has been noted by several previous authors.Flock et al. (2017Flock et al. ( , 2020) ) presented detailed 3D radiation hydrodynamical models of protoplanetary disks with dust particle dynamics.They show that 0.1 mm and 1 mm dust grains achieve greater elevations above the midplane than expected from isotropic turbulence.Similar conclusions were also made by Lehmann & Lin (2022), who show the dependency of this effect on Stokes number St 0 and vertically integrated dust-to-gas ratio Z big , although they focus on smaller values of the Stokes number than we explore in this paper.However, we put this into context with recent observational evidence of the extreme vertical geometrical thinness of the big-grain dust layers in (most?) protoplanetary disks.

Effect of dust traps
The fact that our models are for disks without dust traps limits the applicability of the results.If the large dust grains in the outer regions of protoplanetary disks remain at those large distances because they are trapped, then it is rather natural to get Z big ≫ 0.01 in these traps because all the dust elsewhere will radially drift into these traps, enhancing Z there.As argued by Lin (2019) and Lehmann & Lin (2022), this then could naturally push Z big ≳ 0.02 • • • 0.05 which, according to their simulations, strongly suppresses the VSI.
The fact that the upper limit of Z big,VS I for the VSI lies around the same value as the lower limit of Z big,S I for streaming instability (Carrera et al. 2017-04) leads to an interesting speculation.It was shown by Stammler et al. (2019), using an argument related to that of Sect.3, that Z big,S I coincides with an optical depth at millimeter wavelengths of order unity, as appears to be observed in ALMA observations.They argue that dust traps attract more and more dust, until Z big reaches Z big,S I , at which point Z big stabilizes: Any further dust added to the trap will be converted into planetesimals by the streaming instability, keeping Z big = Z big,S I .If Z big,S I > Z big,VS I , then this self-regulating system naturally keeps the disk VSI-stable.This could be another natural explanation for the lack of VSI, but this requires more detailed study of the combined VSI+SI, as in Schäfer et al. (2020).So at this point, this is merely a speculative idea.

Effect of small but non-zero background turbulence
It was noted by Nelson et al. (2013) that the VSI is also damped if the viscosity parameter of the disk α ≳ 4 × 10 −4 , i.e., typically when the disk is turbulent due to the magnetorotational instability (MRI).While MRI turbulence will also stir up large grains away from the midplane, it is far less effective than the VSI.And so, somewhat paradoxically, the existence of weak, but no-zero, turbulence might, by inhibiting the VSI, allow large grains to settle to a thinner layer than is the case for a non-turbulent disk.

3D effects
Since our models are 2D in (r, z), we cannot treat any potential non-axisymmetric modes, such as the formation of long-lived vortices (Lehmann & Lin 2022).The main effect of the VSI acts in the (r, z) directions, however, and is not dependent on the ϕ-direction.Any 3D effects on the large-scale may affect the observational appearance of the disk, but will likely not affect our conclusion that the flat disks seen with ALMA are incompatible with the VSI.

Small-scale modes
The VSI may operate on smaller scales than we can model with our global models, for example, via the parametric instability mechanism described by Cui & Latter (2022).This can have consequences for the dust dynamics and dust growth.If the VSI operates on the large scales as explored in this paper, the observational consequences will be dominated by these largescale motions, even if smaller-scale motions are superposed on them.However, as shown in (Cui & Latter 2022), the small-scale motions excited by the larger ones act as an energy sink to the large-scale motions.The long-term saturation state of the largescale VSI modes may therefore depend on the very small-scale motions that require super-high spatial resolution to resolve.(Cui & Latter 2022) cite a resolution of 300 grid cells per scale height, which is out of reach for global simulations.These considerations show that it remains to be explored to which degree the VSI or any other instabilities are inhibited if a protoplanetary disk is observed to have a very flat midplane dust layer, and what this means for the implied conditions in the disk.

Uncertainties in the thermal relaxation time
Our estimates of the thermal relaxation time suffer from some uncertainties.First, they are very sensitive to the disk temperature profile T (r).This is because radiative cooling goes as ∝T 4 .But in practice this sensitivity is not so severe, because any uncertainty in the irradiation q heat ∝ L * of the disk only enters the temperature as T ∝ q 1/4 heat ∝ L 1/4 * .It does show, however, that for disks around Herbig Ae stars the relaxation times are smaller than for T Tauri stars.
A much more critical uncertainty is the dust opacity at long wavelengths.The popular Bell & Lin (1994) opacity represents a relatively low estimate, lower than what we use in this paper.This opacity model was used by Lin & Youdin (2015) and Pfeil & Klahr (2019) for their relaxation time estimates, which leads to less favorable conditions for the VSI than our estimates.Malygin et al. (2017) use the opacity model of Semenov et al. (2003), which is, for T mid < 100 K, very similar to that of Bell & Lin (1994).The factor of ∼1.5 difference can mostly by explained by the use of a different dust-to-gas ratio, so that the dust opacities are more or less the same.In contrast, Fukuhara et al. (2021) use a simple analytic opacity model (see Ivezic et al. 1997), which is higher than what we use in this paper, and leads to more favorable conditions for the VSI than our estimates (although this effect is limited by the dust-gas coupling time scale that becomes the limiting factor).The comparison of these opacities is shown in Fig. C.1.The opacity at the far-infrared and submillimeter of the dust in protoplanetary disks is notoriously uncertain, and the real opacity is likely somewhere between these two extremes.
Another major uncertainty is the dust-to-gas ratio.In our analysis we kept this constant at a value of Z = 0.01, although we allowed coagulation to convert the small grain population (responsible for the radiative cooling) into a big grain population (responsible for the dust observed with ALMA).However, these big grains can radially drift, leading to a reduction of the dustto-gas ratio in the outer regions.However, since the big grains do not contribute to the cooling, this does not affect our analysis.What matters is only the small grain abundance Z small .Dust coagulation can reduce Z small .What subsequently happens to Z big is irrelevant for the estimation of t relax .

Non-flat rings
It should be noted that there are observed protoplanetary disks for which one or more of the rings do not appear to be geometrically very thin.For instance, Doi & Kataoka (2021) conclude, after detailed radiative transfer modeling, that the inner dust ring of HD 163296 is, in fact, likely to be vertically extended to about a gas pressure scale height.One interpretation of this could be that in this ring the cooling rate is, in fact, not sufficiently reduced by grain growth, so that the VSI is operational.From Fig. 12 it can be seen that for a small-grain depletion somewhere around 3 × 10 −2 , the inner disk regions remain unstable to the VSI while the outer disk regions stabilize against the VSI, which might explain why the inner ring of HD 163296 is vertically more extended than the outer one.A similar point was made by Fukuhara et al. (2021).The disk of HD 163296 is fairly dim at those wavelengths (Garufi et al. 2022), which does seem to point to substantial small-grain depletion.

Conclusions
In this paper we show that the geometrically very thin midplane layers of dust observed in many protoplanetary disks, most strikingly shown in a recent paper by Villenave et al. (2022), are evidence that the VSI is not operational in these outer disk regions.Dust particles with dimensionless stopping times St ≲ 1 would be stirred up by the VSI to a substantial fraction of the gas pressure scale height, even for large particles with St ≃ 0.1 − 1.Only for even larger particles, with St ≳ 1, does the dust layer remain largely unaffected by the VSI.But we show that to have such a layer being marginally optically thick (τ ≳ 0.3) at ALMA wavelengths, as seems to be the case for many such rings, this requires a vertically integrated dustto-gas ratio Z big ≳ 0.08 • • • 0.16, which we consider an unlikely scenario.
Damping or inhibiting the VSI in the outer regions of a protoplanetary disk can be due to an enhanced vertically integrated dust-to-gas ratio of Z big ≳ 0.02 • • • 0.04, as shown by Lin (2019) and Lehmann & Lin (2022), or due to a modest background turbulence (Nelson et al. 2013).
We show that another possible explanation is that dust coagulation has converted more than 90% of the small grains in the disk into big grains (likely the ones that make up the midplane dust layer).In that case, the gas cannot cool fast enough through the thermal emission of the small grains, and the VSI is inhibited, as shown by Lin & Youdin (2015).Small-grain dust depletion of more than 90% by coagulation (X = 0.9) is reasonable during the life time of these disks (Tanaka et al. 2005;Dullemond & Dominik 2005;Birnstiel et al. 2010), and remains consistent with the SEDs of these disks (Dullemond & Dominik 2004a).
Our conjecture is thus that protoplanetary disks that show geometrically thin disks in ALMA observations are stable against the VSI.Bell & Lin (1994) is plotted for comparison, which is used by Lin & Youdin (2015), Pfeil & Klahr (2019) and many other works.Malygin et al. (2017) use the dust opacity model of Semenov et al. (2003), which is, however, very similar to that of Bell & Lin (1994) in this temperature range.Also for comparison, the top two curves in each panel show the mean opacity for the simple κ λ = 3π/2ρ s λ small-grain opacity model used by Fukuhara et al. (2021), for our value of the material density ρ s = 1.48 g cm −3 (Simple) and their value of ρ s = 1.00 g cm −3 (SimpleF).where the effective scattering opacity κ scat,eff ν,small is defined as in Eq. (B.1).
The resulting mean opacities are shown in Fig. C.1.As one can see, for small enough grains and small enough temperatures, the Planck mean opacity of the dust can be well approximated for T ≲ 100 K and a ≲ 10 µm by the following fitting formula: where g is to be interpreted as gram of small-grain dust.The symbol K is the unit of Kelvin.Likewise the Rosseland mean opacity can be approximated as For comparison, the dusty part of the Bell & Lin (1994) opacity is κ R,BellLin (T ) = 2 × 10 −2 (T/K) 2 cm 2 /g ≃ (T/7.07K) 2 cm 2 /g, where we assume a dust-to-gas ratio of 0.01.The equivalent Planck-mean opacity for Bell & Lin would be κ P,BellLin (T ) = (T/4.58K) 2 cm 2 /g.This means that our opacity model is more favorable to the onset of the VSI than the Bell & Lin opacity, implying that in our analysis we need to reduce the small-grain dust more strongly than the analysis done by Lin & Youdin (2015) and Pfeil & Klahr (2019) to suppress the VSI.The primary reason why our opacity model exceeds that of Bell & Lin is the amorphous carbon mixed into the composition, which is responsible for the "antenna effect", which strongly enhances the long-wavelength opacity.If the carbon would be, instead, in the form of organics, this effect would not be seen.
The 25% porosity also increases the opacity a bit.If we use pure pyroxene without porosity, our opacity would exceed that of Bell & Lin by only 18%, which is well within the uncertainty of the dust-to-gas ratio we used to convert the Bell & Lin total opacity to a dust-only opacity.We refer to Woitke et al. (2016) for details on the role of carbon.

Fig. 1 .Fig. 2 .
Fig. 1.Vertical gas velocity v z in the disk model at time t = 300 P orb (r 0 ), where P orb (r 0 ) is the orbital period at r = r 0 .The coordinates are the natural spherical coordinates of the numerical hydrodynamic model: On the horizontal axis the natural logarithm of the spherical radius r in units of r 0 .On the vertical axis the polar angle π/2 − θ.Blue is upward and red is downward.The gray dotted lines show the gas pressure scale height.The blue dots are the initial locations of the particles, where only every 25th of the 2000 particles is shown.The purple box represents the zoom-in view shown in Fig. 4.

Fig. 3 .
Fig.3.For the fiducial model shown in Figs.1 and 2, the vertical gas velocity v z at r = r 0 and z = 0 in units of the local isothermal sound speed as a function of time in units of orbits after the 300th orbit.

Fig. 4 .Fig. 5 .
Fig.4.Snapshots of the location of the particles (blue dots) 2.5 orbits after they were inserted at the midplane into the fully developed VSI hydrodynamic model, for four values of St 0 .The coordinates are the cylindrical radius R in units of, and relative to, r 0 , and the cylindrical vertical height z above the midplane in units of r 0 .The background image shows the vertical gas velocity v z , where blue is upward and red is downward, in the same color scale as in Fig.1.The gray dotted lines show the gas pressure scale height for comparison.The lightblue dashed lines are the conveyor-belt estimate of the maximum vertical height of the dust particles, Eq. (2) with an upper cap at z = h p (which is why the dashed and dotted lines overlap for St 0 = 0.01).

Fig. 6 .
Fig. 6.Spatial distribution of St 0 = 0.1 dust in the hydrodynamic model in which both gas and dust are dynamically modeled as a fluid (Sect.2.3).The white dashed lines mark one gas pressure scale height above/below the midplane.

Fig. 10 .
Fig. 10.Relaxation time (Eq.(19)) of the disk in units of the Kepler time Ω −1 K for three small-grain dust-to-gas ratios: Normal (Z small = 10 −2 ), depleted by a factor of 10 −1 (Z small = 10 −3 ) and depleted by a factor of 10 −2 (Z small = 10 −4 ).Dotted lines: optically thin approximation (Eq.(17)).Solid lines: including optical depth effects.In dashed black: upper limit to the relaxation time (Eq.(7)) for which VSI is operational is shown.The disk and stellar parameters are those of the star Oph 163131 (the fiducial model of Appendix A).
Fig. C.1.Small dust grain mean opacities as a function of temperature for three grain sizes.Left: Rosseland mean.Right: Planck mean.In the left panel the often-used opacity ofBell & Lin (1994) is plotted for comparison, which is used byLin & Youdin (2015),Pfeil & Klahr (2019) and many other works.Malygin et al. (2017) use the dust opacity model ofSemenov et al. (2003), which is, however, very similar to that ofBell & Lin (1994) in this temperature range.Also for comparison, the top two curves in each panel show the mean opacity for the simple κ λ = 3π/2ρ s λ small-grain opacity model used byFukuhara et al. (2021), for our value of the material density ρ s = 1.48 g cm −3 (Simple) and their value of ρ s = 1.00 g cm −3 (SimpleF).