A self-consistent model for dust settling and the vertical shear instability in protoplanetary disks
Yuya Fukuhara Satoshi Okuzumi







1Department of Earth and Planetary Sciences, Tokyo Institute of Technology, Meguro, Tokyo 152-8551, Japan
fukuhara.y.ab@m.titech.ac.jp

 protoplanetary disks — hydrodynamics — instabilities — turbulence 

A self-consistent model for dust settling and the vertical shear instability in protoplanetary disks
    Satoshi Okuzumi1
    
====================================================================================================




    The spatial distribution of dust particles in protoplanetary disks affects dust evolution and planetesimal formation processes. 
    The vertical shear instability (VSI) is one of the candidate hydrodynamic mechanisms that can generate turbulence in the outer disk region and affect dust diffusion. 
    Turbulence driven by the VSI has a predominant vertical motion that can prevent dust settling. 
    On the other hand, the dust distribution controls the spatial distribution of the gas cooling rate, thereby affecting the strength of VSI-driven turbulence. 
    Here, we present a semi-analytic model that determines the vertical dust distribution and the strength of VSI-driven turbulence in a self-consistent manner. 
    The model uses an empirical formula for the vertical diffusion coefficient in VSI-driven turbulence obtained from our recent hydrodynamical simulations. 
    The formula returns the vertical diffusion coefficient as a function of the vertical profile of the cooling rate, which is determined by the vertical dust distribution. 
    We use this model to search for an equilibrium vertical dust profile where settling balances with turbulent diffusion for a given maximum grain size.  
    We find that if the grains are sufficiently small, there exists a stable equilibrium dust distribution where VSI-driven turbulence is sustained at a level of α_z ∼ 10^-3, where α_z is the dimensionless vertical diffusion coefficient. 
    However, as the maximum grain size increases, the equilibrium solution vanishes because the VSI can no longer stop the settling of the grains. 
    This runaway settling may explain highly settled dust rings found in the outer part of some protoplanetary disks.








§ INTRODUCTION


The initial stage of planet formation is the formation of kilometer-sized planetesimals from micron-sized dust grains (for reviews, ).
These dust growth and planetesimal formation depend on the vertical distribution of dust particles.
If the dust density at the midplane well exceeds the gas density, the streaming and gravitational instabilities set in (e.g., ), leading to planetesimal formation.
Understanding what processes determine the dust vertical profile in disks is also essential for interpreting recent high-resolution radio observations showing that the dust at different radial locations or in different disks is settled to different extents (, and for a review, ).

In protoplanetary disks, the vertical profile of dust particles critically depends on the intensity of gas disk turbulence causing dust diffusion.
The question then is what mechanisms generate disk turbulence.
Recent theoretical studies have shown that not only the magnetorotational instability (MRI; ) but also some thermo-hydrodynamical instabilities are important turbulence-driving mechanisms (for reviews, ).
Among them, the vertical shear instability (VSI; ) is thought to be the leading mechanism for driving turbulence in outer disk regions where the current radio observations best constrain the degree of dust settling.
The VSI requires rapid cooling of disk gas in addition to a vertically varying gas orbital velocity <cit.>.
Therefore, the VSI tends to operate in the outer disk region with low optical depths <cit.>, where it may dominate over the MRI (e.g., ).
The VSI generates turbulence with a predominant vertical gas motion (e.g., ) that can prevent dust settling <cit.>.
 This turbulence can also stir dust particles significantly <cit.> and thus induce their collisional velocities, leading to suppression of planetesimal formation through coagulation (e.g, ).
On the other hand, VSI-driven turbulence can produce both small short-lived and azimuthally large long-lived vortices <cit.>.
These vortices may lead to dust concentration and subsequent planetesimal formation through gravitational collapse (e.g., ).

Importantly, dust particles are the dominant opacity source and determine the local cooling rates of protoplanetary disks <cit.>. 
Therefore, their spatial distribution controls where in the disks the VSI operates <cit.> and can even affect the dust vertical diffusivity at the midplane  <cit.>.
The consequence of this dust–VSI thermal interaction for dust settling remains unclear.

In this paper, we model the above-mentioned interaction between dust and the VSI to study how dust setting and diffusion balance in VSI-driven turbulence.
To this end, we present a semi-analytic model that determines the vertical dust distribution and the strength of VSI-driven turbulence in a self-consistent manner. 
The model uses an empirical formula for the vertical diffusion coefficient in VSI-driven turbulence obtained from our recent hydrodynamical simulations <cit.>. 
The formula returns the vertical diffusion coefficient as a function of the vertical profile of the cooling rate, which is determined by the vertical dust distribution. 
We use this model to search for an equilibrium vertical dust profile where settling balances with turbulent diffusion for a given grain size.

This paper is organized as follows. 
In section <ref>, we describe our self-consistent model. 
We present the main results in section <ref> and discuss the implication of our study in section <ref>.
Section <ref> presents a summary.




§ METHOD





In this section, we describe our self-consistent model that determines the dust vertical distribution and VSI-driven turbulence intensity simultaneously (see figure <ref> for an overview).
We assume a protoplanetary disk consisting of gas (section <ref>) and dust (section <ref>).
The dust grains dominate the disk's opacity (section <ref>) and their distribution determines the disk's local cooling rate distribution (section <ref>).
Because the VSI requires rapid gas cooling, the cooling rate distribution in turn determines the location where the linear VSI operates (section <ref>), which we call the VSI-unstable layer[<cit.> called these the VSI zones.].
The VSI produces turbulent gas motion propagating across the boundary between the VSI-unstable and stable regions. 
We calculate the VSI-driven vertical diffusion coefficient from the thicknesses of the VSI-unstable and -stable layers using an empirical formula based on our previous hydrodynamical simulations (section <ref>).
The vertical diffusion coefficients can be used to update the vertical dust distribution. 
Iterating these calculating steps, we determine the vertical dust distribution and the strength of VSI-driven turbulence in a self-consistent manner (section <ref>).
In the following subsections, we describe the model in more detail.

    
    

 §.§ Disk model

    
    We consider an axisymmetric disk around a solar-mass star.
    We adopt the cylindrical coordinate system (R, z), where R and z are the distance from the central star and the height from the midplane, respectively.
    The gas surface density is given by
    
    Σ_ g(R) =  2-β_ΣMdisk/2π R_c^2R/R_c^-β_Σexp[-R/R_c^2-β_Σ] ,

    where Mdisk is the total mass of the gas disk, R_c is the characteristic radius, and β_Σ is a dimensionless number characterizing the radial slope of the gas surface density. 
    We fix the gas disk parameters to Mdisk = 0.01, R_c = 100, and β_Σ = 1.
    
    Assuming that the disk is optically thick to stellar radiation and that the stellar luminosity is equal to the solar luminosity, the temperature of the disk is given by
    
    T(R) = T_0R/1 au^q,

    with T_0 = 130  K and q = -3/7  <cit.>. 
    We assume that the disk is vertically isothermal, ignoring warmer surface layers that are optically thin to the starlight <cit.>. 
    
    From vertical hydrostatic equilibrium, the gas density is given by
    
    ρ_ g(R, z) = Σ_ g/√(2π)H_ gexp(-z^2/2H_ g^2).

    where H_ g = c_ s/Ω_ K is the gas scale height with c_ s and Ω_ K being the isothermal sound speed and Keplerian frequency, respectively.
    The isothermal sound speed is given by c_ s = √(k_ BT/m_ g), where k_ B is the Boltzmann constant and m_ g is the mean molecular mass of the gas.
    The Keplerian frequency is given by Ω_ K = √(GM_*/R^3), where G is the gravitational constant and M_* is the mass of the central star.
    In this study, m_ g and M_* are taken to be 2.3m_ p and 1, respectively, where m_ p is the proton mass.

    
    

 §.§ Dust model

    
    We here describe the dust model that we use to calculate the cooling timescale (see section <ref>). 
    The ratio between the dust surface density Σ_ d and Σ_ g is fixed to the interstellar dust abundance of 1%, whereas the local dust-to-gas ratio is allowed to vary with z considering dust settling.
 
    We consider two grain size-distribution models. 
    The first model is the single-sized model where all dust particles are assumed to have equal size a_1. 
    Formally, the size distribution for the single-sized model can be written as 
    
    dN_ d(a)/da = 3Σ_ d/4πρ_ inta_1^3δ(a-a_1)

    where dN_ d(a)/da is the number surface density per unit particle size a, Σ_ d is the total dust mass surface density, ρ_ int is the grains' internal density, and δ is the delta function.
    Equation (<ref>) satisfies the normalization
    
    Σ_ d = ∫ m_ ddN_ d(a)/da da,
    
    where m_ d = (4π/3)ρ_ inta^3 is the particle mass.
    The second model is the power-law model where the grain size distribution is given by
    
    dN_ d(a)/da =
            {[ (12+3p)Σ_ d/4πρ_ int(a_ max^4+p-a_ min^4+p)a^p,                            a_ min < a < a_ max,;                                              0,                                      otherwise, ].

    where p (≠ -4) is the slope of size distribution, and amin and amax are the minimum and maximum particle sizes, respectively.
    We fix p = -3.5 and amin = 1.
    The maximum particle size serves as a free parameter in this study (see also section <ref>).
    Equation (<ref>) also fulfills equation (<ref>).

    Dust particles settle toward the midplane owing to stellar gravity and diffuse away from the midplane owing to turbulence.
    We assume that the turbulent diffusion coefficient is constant in the vertical direction.
    This approach is valid if the VSI is the dominant source of disk turbulence and determines dust vertical diffusion.
    Even if there is the linear-VSI stable layer near the midplane, VSI-driven turbulence can penetrate the stable midplane layer and thereby produce nearly constant turbulent vertical intensity in the vertical direction <cit.>.
    This situation can be realized when the stable midplane layer thickness is thinner than two gas scale heights, or the thickness of the unstable layer above the stable midplane layer is thicker than a few gas scale heights <cit.>.
    In this study, the results for all parameter ranges satisfy these conditions.
    We note that, however, the diffusion coefficient may not be uniform in the vertical direction when these conditions are broken.
    
    Assuming the balance between settling and diffusion, the vertical distribution of the particles can be written as 
    <cit.>
    
    dn_ d(a,z)/da = C_d(a) exp[-z^2/2H_ g^2- St_ mid(a)/α_zexpz^2/2H_ g^2-1],

    where dn_ d(a,z)/da is the particle number density per unit radius at height z, St_ mid(a) is the Stokes number of the particles at the midplane, α_z is a dimensionless parameter that characterizes the level of the dust vertical diffusion caused by turbulence, and C_d(a) is the normalized constant determined by the condition dN_ d(a)/da = ∫ (dn_ d(a,z)/da)dz.
    Because most dust particles lie at the region of z ≪ H_ g, the exponential factor in equation (<ref>) can be approximated as exp[-z^2/(2H_ d^2)], yielding <cit.>
    
    C_d(a) = 1/√(2π)H_ d(a,α_z)dN_ d(a)/da.

    Here, H_ d(a,α_z) is the scale height of particles with size a given by <cit.>
    
    H_ d(a,α_z) = [1+St_ mid(a)/α_z]^-1/2H_ g.

    The Stokes number is the product of the stopping time and Keplerian frequency. Assuming that the particle radius is smaller than the mean free path of the disk gas molecules, gas drag onto the particles follows Epstein's law, which gives <cit.>
    
    St_ mid(a) = πρint a/2Σ_ g.


    The vertical–size distribution dn_ d(a,z)/da gives the collisional heat transfer and opacity, which control the cooling timescale, as a function of z. The mean travel length of gas molecules colliding with dust particles ℓgd is given by 
    
    ℓgd = (∫π a^2dn_ d/da  da)^-1.

    This determines the timescale of collisional heat transfer (see section <ref>).

    To calculate the cooling time, we use the Plank mean opacity and Rosseland mean opacity defined by
    
    κ_ P = ∫_0^∞κ_ g,λ B_λ (T) dλ/∫_0^∞ B_λ (T) dλ,

    and
    
    1/κ_ R = ∫_0^∞1/κ_ g,λ B_λ (T)/ T dλ/∫_0^∞ B_λ (T)/ T dλ,

    respectively, where κ_ g,λ is the wavelength-dependent opacity per gas mass, λ is the wavelength, and B_λ (T) is the Planck function.
    
    The opacity per gas mass is related to the size distribution of the dust particles as (e.g., )
    
    κ_ g,λ = 1/ρ_ g∫κ_ d,λ(a)m_ ddn_ d/da da,

    where κ_ d,λ(a) is the opacity per dust mass.
    Assuming that a dust particle is an uniform sphere, κ_ d,λ(a) can be written as <cit.>
    
    κ_ d,λ(a) = π a^2/m_ d×{[                   24nkx/(n^2+2 )^2,                              x ≤ 1,; min{8kx/3n(n^3-(n^2-1 )^3/2),0.9 },                              x > 1, ].

    where x=2π a/λ is the size parameter and n and k are the real and imaginary parts of the complex refractive index, respectively.
    We calculate n and k at each wavelength using the Bruggeman mixing rule, assuming that the dust is a mixture of silicate and ice with a mass mixing ratio of 1:1 with no porosity.
    The values of n and k for silicate and ice are taken from <cit.> and <cit.>, respectively.
    The internal density of the dust grains is then ρ_ int=1.46   g cm^-3.

    
    

 §.§ Cooling model

    
    The spatial profile of the gas disk cooling rate depends on the size and spatial distribution of the dust particles <cit.>.
    Using the dust model described in section <ref>, we approximate the local cooling (thermal relaxation) timescale τ_ relax as <cit.>
    
    τ_ relax(z) = max{τ_ diff(z), τ_ coll(z), τ_ emit(z)},

    where τ_ diff, τ_ coll, and τ_ emit are the timescales of radiative diffusion, collisional heat transfer, and radiative cooling, respectively.

    The cooling time in the optically thick regime is dominated by the radiative diffusion timescale τ_ diff as <cit.>
    
    τ_ diff = 1/D̅k^2,

    where D̅ and k are the effective energy diffusion coefficient and the wavenumber of the perturbation, respectively.
    The effective energy diffusion coefficient is given by
    
    D̅ = λ_ rc/κ_ R(T) ρ_ g4η/1+3η

    where λ_ r is the flux limiter, c is the light speed, and η is the ratio between radiation energy density E_ r to combined radiation and internal energy density E_ int.
    Here, η is given by η = E_ r/(E_ r + E_ int) with E_ r = 4σ_ SB T^4/c and E_ int = ρ_ g C_V T, where σSB is the Stefan–Boltzmann constant and C_V = 5k_ B/2m_ g is the specific heat at constant volume.
    The flux limiter is fixed to λ_ r = 1/3, which is the value for the optically thick limit <cit.>.
    For k, we assume that the thermal perturbation is determined by the wavenumber of the VSI-driven turbulence structure.
    The VSI unstable modes emerge when a radial perturbation wavenumber is larger than a vertical one (e.g., ), and typically takes ∼ 20/H_ g by previous simulations of VSI-driven turbulence[ Note that the radial wavenumber of the VSI-driven turbulence structure can depend on the density and opacity. The previous hydrodynamical simulations including radiative transport by <cit.> have shown that an increase in the density leads to a smaller radial wavelength of turbulence structure.] (see appendix of ).
    Therefore, we set k=20/H_ g.

    The timescale of collisional heat transfer is given by <cit.>
    
    τcoll = ℓgd/vth,

    where vth is the mean relative velocity between the gas molecules and dust particles. 
    The relative velocity vth can be approximated as the mean thermal speed of the molecules, 
    
    vth = √(8 T/π m_ g).

    We note that this collisional timescale is an approximation to the actual thermal accommodation timescale of the gas molecules and dust particles, which differs by a factor of γ/(γ-1) <cit.>, where γ is the heat capacity ratio.
    
    The radiative cooling timescale τemit in the optically thin limit is given by <cit.>
    
    τemit = C_V/16κ_ P(T)σSB T^3.

    Both ℓgd and κ_ P depend on the local size distribution of the dust particles (see section <ref>). 

    
    

 §.§ Defining the VSI-unstable layer

    
    The cooling timescale determines the VSI-unstable layer because rapid cooling reduces buoyancy preventing instability growth (e.g., ). 
    Following <cit.>, the criterion for the VSI can be expressed as
    
    τ_ relax(z) ≲τ_ crit.

    Here, τ_ crit is the vertically global critical cooling timescale defined by
    
    τ_ crit = H_ g/R|q|/γ-1Ω_ K^-1,

    with γ = 1.4.
    
    By applying the VSI criterion of equation (<ref>) to each point in the disk, we search for the VSI-unstable and -stable layers.
    We apply this vertically global criterion instead of the local criterion [equation (4) in ].
    This is because the intensity of VSI-driven turbulence is tightly correlated with the thicknesses of the VSI-unstable and -stable layers predicted by the global criterion (, see also section <ref>).
    The unstable and stable layers exist on the upper (z>0) and lower (z<0) halves of the disk, symmetrically concerning the midplane (z=0, see the middle panel of figure <ref>).
    For the upper halves, we determine the height of the unstable layer's upper boundary z_ u and the height of the midplane stable layer's upper boundary z_ s.
    When the midplane stable layer is absent, we set z_ s = 0.
    Following <cit.>, we also define the thicknesses of the linearly stable and unstable layers as Δ L_ s = 2z_ s and Δ L_ u = 2z_ u-Δ L_ s, respectively (see figure <ref>).

    
    

 §.§ Estimating dust vertical diffusion coefficient

    
    The thicknesses of the VSI-unstable and -stable layers can be used to determine the vertical diffusion of dust particles <cit.>.
    VSI-driven turbulence generally exhibits the vertical gas motion that is uniform in the vertical direction (e.g., ).
    Therefore, we estimate the dimensionless dust vertical diffusion coefficient of VSI-driven turbulence α_z, VSI as
    
    α_z, VSI = ⟨ v_z^2 ⟩/c_ s^2·τ_ corrΩ_ K,

    where ⟨ v_z^2 ⟩ is the time-averaged squared gas vertical velocity at the midplane and τ_ corr is the correlation time of VSI-driven turbulence.
    The correlation time of VSI-driven turbulence is still indeterminate, ranging from τ_ corrΩ_ K∼ 0.2 to ∼ 20 (, see also section 4.1 of ).
    In this study, we fix τ_ corrΩ_ K = 1.0.

    

    Global axisymmetric simulations by <cit.> show that the mean squared gas vertical velocity ⟨ v_z^2 ⟩ in VSI-driven turbulence is tightly correlated with the thicknesses of the VSI-unstable and -stable layers.
    We estimate ⟨ v_z^2 ⟩ from Δ L_ u and Δ L_ s, using an empirical formula provided by <cit.>,
    
    ⟨ v_z^2 ⟩/c_ s^2= f_ T(Δ L_ u, Δ L_ s) + f_ pT(Δ L_ u, Δ L_ s),

    where f_ T and f_ pT represent the dependence of Δ L_ u and Δ L_ s on ⟨ v_z^2 ⟩.
    For f_ T and f_ pT, we use equations (17) and (19) of <cit.>, which represent the sharp decrease of ⟨ v_z^2 ⟩ in Δ L_ u≲ 2H_ g and Δ L_ s≳ 2H_ g.
    In figure <ref>, We plot α_z, VSI [equation (<ref>)] as a function of Δ L_ u and Δ L_ s.
    When Δ L_ u≳ 2H_ g and Δ L_ s≲ 2H_ g, equation (<ref>) predicts α_z, VSI≈ 2× 10^-3; otherwise, α_z, VSI approaches zero.

    
    

 §.§ Calculation procedure

    
    The vertical dust distribution depends on the turbulent diffusion coefficient α_z [equation (<ref>)]. 
    On the other hand, if we assume that the VSI is the main driver of disk turbulence, its diffusion strength α_z, VSI depends on the disk's cooling structure [equation (<ref>)] and in turn on vertical dust distribution. 
    In general, α_z, VSI produced by VSI-driven turbulence under a given vertical dust distribution does not necessarily match α_z required to maintain the dust distribution. 
    We search for an equilibrium state where the two turbulence strengths match in the following steps.
    
    
        
  * For a given grain size distribution and radial position R, use equation (<ref>) to generate vertical dust profiles for various values of trial turbulence strengths α_z = α_z, assume. The trial values are generated by dividing the range 10^-8≤α_z, assume≤ 10^-2 into grids of an equal logarithmic interval of 0.05.
        
  * For each value of α_z, assume, use the cooling model [equation (<ref>)], VSI criterion [equation (<ref>)], and the empirical formula for VSI-driven turbulence strength [equation (<ref>)] to evaluate α_z, VSI.
        
  * Interpolate α_z, VSI as a smooth function of α_z, assume, and search for self-consistent equilibrium solutions where the predicted VSI-driven turbulence strength α_z, VSI equals the assume turbulence strength α_z, assume. 
        Below, we denote the vertical diffusion coefficient for an equilibrium solution by α_z, equi. If α_z, VSI<α_z, assume for all values of α_z, assume, VSI-driven turbulence cannot be strong enough to sustain the assumed dust vertical distribution, meaning that it cannot stop grain settling. We call this runaway dust settling.
    
    For every set of grain size and R, we repeat this procedure and search for equilibrium solutions.

    
    We ignore that VSI-driven turbulence revives before the dust settles toward the midplane in the runaway fashion.
    This assumption is valid if the dust-settling timescale is longer than the growth timescale of turbulence.
    The settling timescale can be estimated as ∼ z/|v_z| ≈ St_ mid^-1Ω_ K^-1 (e.g., ), where v_z ≈ - St_ midΩ_ Kz is the vertical velocity of dust particles for the terminal velocity approximation.
    The typical growth timescale of VSI-driven turbulence is ∼ 10^2Ω_ K^-1 (e.g., ), which is longer than the settling timescale for large dust grains with St_ mid≳ 10^-2.
    However, for small grains with St_ mid≲ 10^-2, the growth timescale is shorter than the settling timescale, suggesting re-development of VSI-driven turbulence.
    In this case, weak turbulence may operate before the dust settles.

    
    

 §.§ Parameter choices

    
    Our model involves two main parameters: the size of dust grains a_1 for the single-sized model or the maximum size of dust grains a_ max for the power-law size distribution model, and the distance from the central star R.
    We divide the parameter ranges  1  < a < 1 cm or 1  < a_ max < 1 cm and 1 au < R < 100 au into grids of an equal logarithmic interval of 0.1.
    



§ RESULTS



In this section, we use the model described in section <ref> to search for equilibrium solutions where dust settling balances with VSI-driven turbulent diffusion.
In section <ref>, we describe the properties of the equilibrium solutions and their dependence on the dust grain size and radial distance for the single-sized model.
In section <ref>, we show the equilibrium solutions for the power-law size distribution model.

    
    

 §.§ Equilibrium solutions for the single-sized model

    

    

    First, we use the single-sized model at R=10 au to illustrate the properties of the equilibrium solutions.
    Following the procedure described in section <ref>, we search for the equilibrium vertical diffusion coefficient α_z, equi for different values of single grain size a_1.
    The upper panel of figure <ref> shows α_z, equi for each a_1.
    We find that when a_1≲ 10, there exist two equilibrium solutions.
    The equilibrium solutions vanish when the grain size exceeds ∼ 10.
    We also plot in the lower panel of figure <ref> the dust scale height with equilibrium solutions defined by H_ d,equi(a) ≡ H_ d(a,α_z, equi) as a function of a_1.

    Figure <ref> illustrates two cases at R = 10  au with  a_1 = 10 and 100, representing the cases with two and no equilibrium solutions, respectively.
    The upper panel plots α_z, VSI as a function of α_z, assume.
    For a_1 = 10, α_z, VSI drops sharply from ≈ 2× 10^-3 to ≪ 10^-3 as α_z, assume falls below ∼ 10^-4.
    A decrease in α_z, assume leads to dust depletion at high altitudes, yielding a decrease in the thickness of the VSI-unstable layer.
    This decrease in the unstable layer's thickness results in a strong suppression of VSI-driven turbulence around Δ L_ u≈ 2H_ g (see figure <ref>).
    For a_1 = 100, α_z, VSI takes lower values of ≲ 10^-10.
    This is because the VSI-unstable layer is thin (Δ L_ u≲ 2H_ g) for any values of α_z, assume due to dust vertical settling, resulting in suppression of VSI-driven turbulence.
    
    Equilibrium solutions correspond to the points where α_z, VSI = α_z, assume.
    For a_1=10, two equilibrium solutions exist at α_z, equi≈ 2× 10^-3 and 2× 10^-5.
    For a_1=100, no equilibrium solution exists because α_z, VSI<α_z, assume for all values of α_z, assume.
    This implies that the vertical diffusion by VSI-driven turbulence alone cannot sustain dust vertical distribution of any scale height in the latter case.

    To understand the results presented in the upper panel of figure <ref> in terms of the competition between vertical settling and turbulent diffusion, we consider a time-dependent toy model where the vertical distribution of dust grains evolves through settling and diffusion. 
    The equation that describes the evolution of the squared mean of the grains' vertical positions, ⟨ z^2 ⟩, is given by (for a derivation, see appendix <ref>)
    
    1/2d⟨ z^2⟩/dt = - St_ midΩ_ K⟨ z^2⟩ + D_z(α_z, VSI)(1-⟨ z^2⟩/H_ g^2),

    where D_z(α_z) = α_zc_ sH_ g is the dust vertical diffusion coefficient.
    Below we simply call ⟨ z^2⟩^1/2 the dust scale height. 
    In the right-hand side of equation (<ref>), the first term represents settling toward the midplane at the terminal vertical velocity - St_ midΩ_ Kz. 
    The second term corresponds to the vertical diffusion by VSI-driven turbulence in a disk of gas scale height H_ g, with the factor (1-⟨ z^2⟩/H_ g^2) guaranteeing that the dust scale height never exceeds H_ g in the limit of strong diffusion <cit.>.
    It can be shown from equation (<ref>) that our single-sized model assumes
    
    ⟨ z^2 ⟩ = [H_ d(a,α_z, assume)]^2 = [1+ St_ mid(a)/α_z, assume]^-1 H_ g^2

    for an arbitrary value of α_z, assume.
    Substituting this into equation (<ref>), we find that d⟨ z^2⟩/dt=0 if α_z, assume=α_z, VSI, meaning that the dust settling and diffusion indeed balance in the equilibrium solutions defined in section <ref>. 
    This toy model also allows us to predict what would occur when α_z, assume and α_z, VSI are unequal. 
    For example, if α_z, assume > α_z, VSI (i.e., if VSI-driven turbulence is not strong enough to sustain the vertical dust distribution), equation (<ref>) shows that settling dominates over diffusion and the dust scale height decreases.

    

    We now use equation (<ref>) to explain why either two or no equilibrium solutions emerge depending on the value of a_1.
    The lower panel of figure <ref> shows d⟨ z^2⟩/dt as a function of ⟨ z^2⟩ for a_1 = 10 and 100.
    In general, there exists one equilibrium solution for a constant value of D_z because d⟨ z^2⟩/dt is then a monotonically decreasing function of ⟨ z^2⟩, i.e., diffusion  (settling) tends to dominate at small (large) dust scale heights.
    This is illustrated by the dashed curves, which show d⟨ z^2⟩/dt in the case where turbulence with α_z, VSI≈ 2× 10^-3 is present independently of ⟨ z^2⟩. 
    Therefore, multiple or no equilibrium solution occurs only when D_z is a function of ⟨ z^2⟩.  
    For a_1 = 10, the two equilibrium solutions have dust scale heights of ⟨ z^2⟩^1/2∼ H_ g and 0.3H_ g. 
    Comparison with the dashed curve shows that the solution with higher ⟨ z^2⟩^1/2 corresponds to the equilibrium solution for constant diffusion coefficient with α_z, VSI≈ 2× 10^-3.     
    The solution with lower ⟨ z^2⟩^1/2 emerges because VSI-driven turbulent diffusion is suppressed for small dust scale heights, i.e., for small α_z,  assume, as mentioned earlier. 
    For a_1 = 100, d⟨ z^2⟩/dt is negative for all dust scale heights, indicating runaway dust settling.
    
    Equation (<ref>) can also be used to predict whether the equilibrium solutions are stable or not against a small perturbation in the dust scale height.
    In general, equilibrium solutions can be classified into stable and unstable ones.
    In this context, stable (unstable) solutions are the ones for which d⟨ z^2⟩/dt becomes negative (positive) as we slightly increase ⟨ z^2⟩ from the equilibrium value, so that the perturbation decays (diverges). 
    In other words, stable (unstable) solutions have a negative (positive) slope of d⟨ z^2⟩/dt as a function of ⟨ z^2⟩. 
    In the case of a_1 = 10 μ m, the equilibrium solution with the larger ⟨ z^2⟩ is stable, whereas that with the smaller ⟨ z^2⟩ is unstable (see figure <ref>).

     

    

    Next, we describe how the stable equilibrium solution varies with radial distance R. 
    Figure <ref> shows α_z, equi for the stable solution as a function of R and a_1 for the single-sized model. 
    The black area in the R–a_1 plane indicates the parameter space where no equilibrium solution exists and runaway dust settling is expected.
    The figure shows that equilibrium solutions only exist beyond R ∼ 3  au. 
    At R ≲ 3  au, the midplane region is significantly optically thick to its own thermal emission, yielding a thick VSI-stable layer (Δ L_ s≳ 2H_ g) around the midplane. 
    We map in figure <ref> the unstable and stable layers for a_1 = 10 with α_z, assume = 2× 10^-3 from the single-sized model at the R–z plane. 
    This figure indicates that the thickness of the stable midplane layer increases sharply as R decreases from ∼ 3 au.
    This thick stable layer prevents VSI-driven turbulence from developing (α_z, VSI≪ 10^-3; see figure <ref>), resulting in runaway dust settling (α_z, VSI<α_z, assume for all values of α_z, assume).  
    For the cases of large dust, optical depth decreases and the VSI-stable layer becomes thin.
    However, because the dust settles vertically, the VSI-unstable layer is also thin (Δ L_ u≲ 2H_ g), leading to suppression of VSI-driven turbulence.
    The stable equilibrium solutions are generally accompanied by fully developed VSI-driven turbulence with α_z, VSI≈ 2× 10^-3. 
    This is reasonable because any vulnerable turbulence whose strength decreases steeply with decreasing dust scale height would lead to an unstable solution.
    
    At R ≳ 3  au, equilibrium solutions exist for sufficiently small a_1 as we already showed in the particular case of R = 10  au. 
    The maximum grain size for equilibrium decreases with increasing R.
    This is because the dust number density decreases as R increases: decreasing dust density leads to a longer cooling time, which makes the VSI-unstable layers vertically thinner ().

    

    Figure <ref> shows the radial profiles of the dust scale height from the stable equilibrium solutions, H_ d,equi(a_1), for different values of a_1.
    We set H_ d,equi(a_1) to zero if runaway settling occurs.
    As mentioned above, the stable solutions are generally accompanied by fully developed VSI-driven turbulence with α_z, equi≈ 2× 10^-3. 
    For reference, the dotted lines in the figure show the dust scale height for α_z, equi = 2× 10^-3. 
    Because the stable solutions require small a_1, they always lead to a thick dust disk with H_ d,equi(a_1) ≈ H_ g.

    Because the dust scale height is a function of St_ mid(a_1), we can expect that the value of St_ mid(a_1) rather than a_1 more directly determines whether there exist equilibrium solutions or not. 
    To confirm this expectation, we overplot in figure <ref> contours of constant St_ mid(a_1). 
    We find that the equilibrium solution at R ≳ 3  au vanishes for St_ mid(a_1) ≳ 3× 10^-4. 

    
    

 §.§ Equilibrium solutions for the power-law size distribution model

    

    The grain size distribution affects the cooling timescale profile.
    Of the three timescales determining the cooling timescale (see section <ref>), the timescales of collisional heat transfer τ_ coll and radiative diffusion τ_ diff determine the thicknesses of the unstable and midplane stable layers, respectively.
    For the power-law size distribution with a slope of p=-3.5, the largest grains dominate the dust mass budget, whereas the smallest grains dominate the total geometric cross-section and thereby control collisional heat transfer.
    Therefore, even if large grains exist, small grains can still sustain the unstable layer.

        

    Figure <ref> illustrates how the small grains in the power-law size distribution model extend the parameter space where strong VSI-driven turbulence is sustained. 
    This figure shows α_z, equi for the stable equilibrium solutions as a function of R and a_ max for the power-law size distribution model.
    At R = 10 au, the equilibrium solutions exist up to a_ max = 100 μ m. 
    This is in contrast to the single-sized model, which only allows equilibrium solutions up to 10 at the same radial position (see figure <ref>).
    In terms of the Stokes number, the power-law size distribution model allows the equilibrium solutions up to St_ mid(a_ max) ≈ 10^-3–10^-2 depending on R. 
    As in the single-sized model, the stable solution is accompanied by fully developed VSI-driven turbulence with α_z, equi≈ 2× 10^-3.
    At R≲ 3  au or at sufficiently large a_ max, the equilibrium solution vanishes and runaway settling (α_z, assume > α_z, VSI) occurs for the same reason as in the single-sized model.

    

    Figure <ref> plots the dust scale height for the stable equilibrium solution as a function of R from the power-law size distribution model with different values of a_ max. 
    As in figure <ref>, we set H_ d, equi = 0 for runaway settling.
    For a_ max=10, VSI-driven turbulence maintains a thick dust layer of H_ d,equi≳ 0.8H_ g at 2 au≲ R ≲ 40 au.
    As a_ max increases, the region with the stable equilibrium solution shrinks, and the equilibrium dust scale height within that region decreases because the settling velocity increases while keeping α_z, equi≈ 2× 10^-3.
    For a_ max > 1  mm, the equilibrium solution vanishes at all R.
    



§ DISCUSSION


    
    

 §.§ Implications for disk observations

    

    In section <ref>, we have shown that the dust grain size and radial distance determine whether the VSI can sustain high dust diffusion.
    VSI-driven turbulence can maintain a thick dust layer of H_ d≈ H_ g within an annular region with a moderate optical depth and small grains.
    However, when both the dust size and radial distance are large, VSI-driven turbulence is suppressed, resulting in runaway dust settling leading to H_ d≈ 0.
    
    These differences in dust settling levels may provide an explanation for varying degrees of dust settling as inferred from recent radio interferometric observations of some protoplanetary disks.
    The disks around HL Tau and Oph 163131 exhibit well-defined dust gaps, limiting the dust scale height to ≲ 0.1H_ g at R ≳ 10 au <cit.> and at R≈ 100 au <cit.>, respectively.
    These estimated dust scale heights are in line with runaway dust settling in the outer disk regions predicted in this study.
    In contrast, the disk around HD 163296 has two major dust rings at 70 and 100 au that show high and low degrees of dust settling with H_ d/H_ g≳ 0.8 and ≲ 0.1, respectively <cit.>.
    This difference in dust diffusion by radial distance is consistent with the trend in the radial profile of dust scale heights predicted by our model (see figures <ref> and <ref>).
    Therefore, we hypothesize that VSI-driven turbulence dominates vertical dust diffusion in these disks, and that this turbulence is suppressed in the outer disk regions.
    Testing this hypothesis requires detailed modeling of gas, dust, and cooling rate profiles in these disks.
    
    
    

 §.§ Implications for disk and dust evolution

    

    In section <ref>, we have shown that there exists a parameter space where VSI-driven turbulence is too weak to stop dust settling. 
    This occurs either when the VSI-stable layer at the midplane is too thick or when the VSI-unstable layer is too thin for the VSI to operate. 
    The former condition is met at small orbital radii (R ≲ 3 au), whereas the latter condition is met when the grains grow beyond a certain size.
    
    This runaway settling is beneficial for dust growth and planetesimal formation.
    A high degree of dust settling generally promotes planetesimal formation through the streaming and gravitational instabilities (e.g., ).
    Weak turbulence leading to high dust concentration is also conducive to dust growth through coagulation without collisional fragmentation and erosion (e.g., ).
    
    Figures <ref> and <ref> imply that strong dust diffusion by the VSI tends to occur in an annular region where the optical depth is moderate (see also figure <ref>). 
    The vertically extended dust in this annulus may block the radiation of the central star and thereby cast a shadow beyond the annulus <cit.>. 
    The shadowing generally results in a significant drop in temperature and therefore can affect the chemical evolution in that region <cit.>.

    The inner and outer edges of the VSI-turbulence zone also serve as potential sites for planetesimal formation. 
    A steep radial change in turbulence viscosity at the edges may trigger the Rossby wave instability <cit.> and create a long-lived vortex, promoting dust concentration and subsequent planetesimal formation through gravitational collapse (e.g., ).
    The sharp drop in turbulent diffusivity at the inner edge of the VSI-driven turbulence zone can also trigger a runaway pile-up of dust grains <cit.> because the increase of dust-to-gas ratio at the midplane reduces the radial drift velocity of dust.
    Moreover, if the turbulent viscosity caused by the VSI dominates gas disk accretion, the outer edge of VSI-driven turbulence may create a local maximum in the radial profile of the gas pressure because the turbulent viscosity decreases sharply toward the outside.
    This is similar to the mechanism by which a pressure maximum is generated near the dead-zone inner edge of MRI (e.g., ).
    The pressure maximum can trap dust particles <cit.>.
    All these dust concentration processes could lead to planetesimal formation via the streaming and the gravitational instabilities (e.g., ).
    We plan to quantify the impact of these dust concentration processes on planetesimal formation by conducting hydrodynamical simulations around the edge of the VSI-driven turbulence zone.

    
    

 §.§ Effects of minimum grain size, grain size distribution, and disk mass

    
    
    The results presented in section <ref> depend on the minimum grain size and slope of the size distribution because the small grains dominate the cooling rate.
    In this study, we have fixed the minimum size of grains and slope of grain size distribution to a_ min = 1 and p=-3.5, respectively.
    Smaller dust grows quickly through Brownian motion <cit.>, but its size limit is uncertain, approximately ranging from 0.1 to 1.
    As a_ min decreases, the cooling timescale decreases, leading to a vertically more extended VSI-unstable layer.
    The grain size distribution limited by the radial drift can also become significantly steeper (e.g., ), with p reaching approximately -2.5.
    As p increases, cooling would be less efficient because the number of the smallest grains is smaller, leading to more suppressed the VSI.
    
    Furthermore, variations in the disk mass can alter the cooling rate profile.
    In this study, we have fixed the disk mass to M_ disk=0.01M_⊙.
    The disk mass can depend on disk age <cit.>; in particular, the disk mass of the young disk around HL Tau can be estimated as M_ disk∼ 0.1M_⊙ <cit.>.
    The massive disk would extend the parameter space where the equilibrium solutions exist because the larger disk mass can make cooling more efficient, leading to the radial and vertical expansion of the VSI-unstable layers (see section 4.3 and figure 10 of ).
    However, as the disk mass increases, the region with no equilibrium solutions due to the thick stable layer around the midplane, corresponding to the regions with R≲ 3   au for figures <ref> and <ref>, would be more extended to larger radial distance.
    This is because as the dust density increases, optical depth increases, leading to inefficient cooling.
    
    
    

 §.§ Limitations of the model

    
    So far we have assumed the local model in the radial direction, meaning that gas and dust remain stationary radially.
    However, VSI-driven turbulence can diffuse dust radially, altering its spatial distribution <cit.>.
    This effect implies that the edges of the VSI-driven turbulence zone evolve with dust.
    Quantifying this effect should be studied in hydrodynamical simulations that include the dynamic and thermal coupling between gas and dust.

    In the region where dust settles in a runaway fashion, other mechanisms that drive disk turbulence may also contribute to the vertical diffusion of dust grains.
    The MRI can drive gas turbulence in both the inner and outer regions of protoplanetary disks because of the high ionization of gas (for reviews, ).
    The streaming instability caused by strong dust settling can also trigger turbulent gas motion <cit.>.
    Turbulence driven by these mechanisms can maintain a vertical dust distribution that is balanced between turbulent diffusion and settling.
    Moreover, the dust vertical diffusion due to these other mechanisms may revive VSI-driven turbulence through changes in the vertical profiles of dust particles and the cooling rate.
    To understand this, we should investigate which mechanism dominates turbulence generation in each region of protoplanetary disks.
    
    Furthermore, this study ignores the effects of dust, magnetic field, and vertical thermal structure.
    Dust would increase the effective buoyancy frequency of the gas <cit.> that prevents the growth of the linear VSI.
    The gas–dust drag force can also suppress the VSI and dust vertical diffusion <cit.>.
    Moreover, magnetic fields threading the global disk may suppress the VSI either directly through magnetic tension or indirectly through MRI turbulence <cit.>.
    The roles of magnetic fields in the VSI suppression can be positive or negative depending on non-ideal magnetohydrodynamical effects (ambipolar diffusion, Ohmic resistivity, and Hall effect; ).
    Additionally, <cit.> recently found that vertical thermal stratification with a colder interior and a hotter surface can suppress VSI-driven turbulence around the midplane.
    They may change the levels of the equilibrium vertical dust profile.




§ SUMMARY


We have searched for the equilibrium vertical dust profile where settling balances with diffusion caused by VSI-driven turbulence.
We construct the semi-analytic model that determines the vertical profile of dust grains and the intensity of VSI-driven turbulence in a self-consistent manner (figure <ref>).
Our key findings are summarized as follows.

    
  * We find that there exist equilibrium solutions where dust settling balances with VSI-driven turbulent diffusion for small grains (figure <ref>). If we assume that all grains have equal size, there exist two equilibrium solutions when the single grain size is smaller than 10 at 10 au. If the grain size exceeds ∼ 10, the equilibrium solutions vanish.
    
  * For the cases of small grains, two equilibrium solutions are classified into stable and unstable ones [equation (<ref>) and figure <ref>]. The stable ones correspond to the dust scale height of α_ z,equi≈ 2× 10^-3, where α_ z,equi is the equilibrium dimensionless vertical diffusion coefficient. For the cases of large grains, no equilibrium solutions indicate runaway dust settling.
    
  * The existence of the equilibrium solutions depends on the radial distance R as well as dust size (figure <ref>). The equilibrium solutions only exist beyond R∼ 3 au because the midplane region at R≲ 3 au is optically thick, yielding a thick VSI-stable layer (figure <ref>). The maximum grain size that allows for the equilibrium solutions also decreases with increasing R. Because the equilibrium solutions require a small grain size, they lead to a thick dust disk with H_ d,equi≈ 0.8H_ g (figure <ref>), where H_ d,equi is the dust scale height with the equilibrium solutions.
    
  * If the particle size distribution is assumed to follow a power law, the small grains extend the parameter space where strong VSI-driven turbulence is sustained (figure <ref>). At 10  au, the equilibrium solutions exist up to a_ max = 100, where a_ max is the maximum size of dust grains. For a_ max=10 in the entire disk, VSI-driven turbulence maintains a thick dust layer of H_ d,equi≳ 0.8H_ g at 2 au≲ R ≲ 50 au (figure <ref>).


Our results suggest that dust diffusion by VSI-driven turbulence has different levels depending on the radial distance.
This variation may explain the different degrees of dust settling inferred from observations of some protoplanetary disks.
This implies that VSI-driven turbulence plays a dominant role in vertical dust diffusion within these disks.
Testing this hypothesis requires a more quantitative investigation of these disks' cooling rate structure.



    We thank Tomohiro Ono for discussions that motivated this project.
    We also thank Akimasa Kataoka, Kiyoaki Doi, Hidekazu Tanaka, Mario Flock, and Takahiro Ueda for the useful discussions of applications to observed protoplanetary disks.
    We appreciate the anonymous referee for comments that greatly helped improve the manuscript.
    This work was supported by JSPS KAKENHI Grant Numbers 
    JP20H01948, JP20H00182, JP22KJ1337, JP23H01227, and JP23K25923.




apj







§ VERTICAL DIFFUSION OF DUST PARTICLES


The evolution equation of the squared mean of the grains’ vertical positions [equation (<ref>)] is derived by considering the flux of dust in the vertical direction.
The vertical diffusion of dust density ρ_ d in a turbulent protoplanetary disk is given by <cit.>

    ρ_ d/ t = - / z(ρ_ dv_z) + / z[ρ_ gD_z/ z(ρ_ d/ρ_ g)],

where v_z is the vertical velocity of dust and D_z is the vertical turbulent diffusion coefficient of dust.
Expanding the second term of the right side in equation (<ref>) and assuming vertical hydrostatic equilibrium for gas [equation (<ref>)] yield <cit.>

    ρ_ d/ t = - / z(ρ_ dv_z) + / z(D_zρ_ d/ z) + / z(D_zz/H_ g^2ρ_ d).

Next, we multiply both sides of this equation by z^2 and integrate it from z = -∞ to z =∞.
The partial integrals are performed repeatedly and surface terms are dropped, resulting in

    1/2/ t∫_-∞^∞ρ_ dz^2 dz = ∫_-∞^∞ρ_ dv_zz dz + D_z( ∫_-∞^∞ρ_ d dz - 1/H_ g^2∫_-∞^∞ρ_ dz^2 dz ).

Assuming the terminal vertical velocity v_z = - St_ midΩ_ Kz and a constant dust surface density defined by Σ_ d≡∫ρ_ d dz, equation (<ref>) can be written as 

    1/2d⟨ z^2⟩/dt = -  St_ midΩ_ K⟨ z^2⟩ + D_z(1-⟨ z^2⟩/H_ g^2),

where ⟨ z^2 ⟩ is the ensemble average of z^2 for dust particles defined by

    ⟨ z^2 ⟩≡∫_-∞^∞ρ_ dz^2 dz/∫_-∞^∞ρ_ ddz.

If we assume the single-sized model and ρ_ d∝exp(-z^2/2H_ d^2), equation (<ref>) results in ⟨ z^2 ⟩^1/2 = H_ d, where H_ d is the dust scale height.


We use the Stokes number at the midplane St_ mid to calculate the vertical velocity.
In reality, the Stokes number increases as z increases, causing the small grains to settle toward the midplane faster.
However, in this study, the Stokes numbers of the smallest grains exceed unity only at z ≳ 4H_ g.
This height is higher than the height of the VSI-unstable layer (z ≲ 3H_ g, see also figure <ref> in section <ref>).
Therefore, the effect of ignoring the regime for larger Stokes numbers than unity has little impact on our results.

