Dust settling against hydrodynamic turbulence in protoplanetary discs

TL;DR

This study uses numerical simulations to explore how hydrodynamic turbulence affects dust settling in protoplanetary discs, revealing conditions under which dust can concentrate sufficiently for planetesimal formation.

Contribution

It demonstrates that dust can settle despite turbulence if certain metallicity and particle size conditions are met, highlighting the role of dust-gas drag forces and metallicity in planetesimal formation.

Findings

Dust particles with Stokes number ~10^{-3} can settle to less than 10% of the gas scale-height.

Super-solar metallicities are necessary for dust clumping via streaming instability.

Dust rings should be thinner than dust gaps if metallicity contrast exceeds gas density contrast.

Abstract

Enhancing the local dust-to-gas ratio in protoplanetary discs is a necessary first step to planetesimal formation. In laminar discs, dust settling is an efficient mechanism to raise the dust-to-gas ratio at the disc midplane. However, turbulence, if present, can stir and lift dust particles, which ultimately hinders planetesimal formation. In this work, we study dust settling in protoplanetary discs with hydrodynamic turbulence sustained by the vertical shear instability. We perform axisymmetric numerical simulations to investigate the effect of turbulence, particle size, and solid abundance or metallicity on dust settling. We highlight the positive role of drag forces exerted onto the gas by the dust for settling to overcome the vertical shear instability. In typical disc models we find particles with a Stokes number $\sim 10^{-3}$ can sediment to $\lesssim 10\%$ of the gas scale-height, provided that $\Sigma_\mathrm{d}/\Sigma_\mathrm{g}\gtrsim 0.02$-$0.05$, where $\Sigma_\mathrm{d,g}$ are the surface densities in dust and gas, respectively. This coincides with the metallicity condition for small particles to undergo clumping via the streaming instability. Super-solar metallicities, at least locally, are thus required for a self-consistent picture of planetesimal formation. Our results also imply that dust rings observed in protoplanetary discs should have smaller scale-heights than dust gaps, provided that the metallicity contrast between rings and gaps exceed the corresponding contrast in gas density.

Equations

$$\rho \equiv {\rho }_{\mathrm{g}}+{\rho }_{\mathrm{d}},$$

$$\mathbf{v}\equiv \frac{{\rho }_{\mathrm{d}}{\mathbf{v}}_{\mathrm{d}}+{\rho }_{\mathrm{g}}{\mathbf{v}}_{\mathrm{g}}}{\rho },$$

$$E\equiv \frac{P}{\gamma -1}+\frac{1}{2}\rho |\mathbf{v}{|}^2+\rho \Phi ,$$

$$P=c_s^2(R){\rho }_{g0}{\left(\frac{{\rho }_{\mathrm{g}}}{{\rho }_{g0}}\right)}^{\gamma }\equiv K(R){\rho }_{\mathrm{g}}^{\gamma },$$

$$c_s^2(R)=c_{s0}^2{\left(\frac{R}{R_0}\right)}^{-q},$$

$$H_{\mathrm{g}}\equiv \frac{c_s}{{\Omega }_{\mathrm{K}}},$$

$$t_{\mathrm{s}}=\frac{{\rho }_{\bullet }s}{\rho c_s}\equiv \frac{{\rho }_0{c}_{s0}}{\rho c_s}\frac{T_{\mathrm{stop}}}{{\Omega }_{K0}}.$$

$$t_{\mathrm{s}}{\Omega }_{\mathrm{K}}\ll 1.$$

$${\mathbf{v}}_{\mathrm{d}}={\mathbf{v}}_{\mathrm{g}}+t_{\mathrm{s}}\frac{\nabla P}{{\rho }_{\mathrm{g}}}$$

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho \mathbf{v}\right)=0,$$

$$\frac{\partial \mathbf{v}}{\partial t}+\mathbf{v}\cdot \nabla \mathbf{v}=-\frac{1}{\rho }\nabla P-\nabla \Phi -2{\Omega }_{K0}\hat{\mathbf{z}}\times \mathbf{v},$$

$$\frac{\partial E}{\partial t}+\nabla \cdot \left[(E+P)\mathbf{v}\right]=\frac{P}{\gamma -1}\mathbf{v}\cdot \nabla \ln c_s^2$$

$$+\frac{\gamma P}{{\rho }_{\mathrm{g}}(\gamma -1)}\nabla \cdot \left(f_{\mathrm{d}}{t}_{\mathrm{s}}\nabla P\right),$$

$$f_{\mathrm{d}}=1-\frac{1}{\rho }{\left(\frac{P}{K}\right)}^{1/\gamma },$$

$$\Phi (R,\varphi ,z)=$$

$${\Phi }_{\ast }(R,z)=-R^2{\Omega }_{\mathrm{K}}^2\left(1-\frac{z^2}{2R^2}\right).$$

$$S\equiv \ln \frac{P^{1/\gamma }}{\rho }=\ln c_s^{2/\gamma }-\ln (1+ϵ)+const.$$

$$N_z^2\equiv -\frac{1}{\rho }\frac{\partial P}{\partial z}\frac{\partial S}{\partial z}=\frac{1}{1+ϵ}\frac{1}{\rho }\frac{\partial P}{\partial z}\frac{\partial ϵ}{\partial z},$$

$$ϵ(R,z)={ϵ}_{\mathrm{mid}}(R)\times \exp \left(-\frac{z^2}{2H_{ϵ}^2}\right).$$

$$\frac{1}{H_{ϵ}^2}\equiv \frac{1}{H_{\mathrm{d}}^2}-\frac{1}{H_{\mathrm{g}}^2}.$$

$$Z\equiv {\frac{{\Sigma }_{\mathrm{d}}}{{\Sigma }_{\mathrm{g}}}|}_0={{ϵ}_{\mathrm{mid}}\frac{H_{\mathrm{d}}}{H_{\mathrm{g}}}|}_0$$

$$R{\Omega }^2$$

$$0$$

$${\rho }_{\mathrm{g}}={\rho }_{\mathrm{g},\mathrm{mid}}\exp \left\{-\frac{z^2}{2H_{\mathrm{g}}^2}-{ϵ}_{\mathrm{mid}}\frac{H_{ϵ}^2}{H_{\mathrm{g}}^2}\left[1-\exp \left(-\frac{z^2}{2H_{ϵ}^2}\right)\right]\right\}.$$

$${\rho }_{\mathrm{g},\mathrm{mid}}(R)=\frac{{\Sigma }_{\mathrm{g}}(R)}{\sqrt{2\pi }{H}_{\mathrm{g}}(R)},$$

$${\Sigma }_{\mathrm{g}}(R)={\Sigma }_{g0}{\left(\frac{R}{R_0}\right)}^{-p},$$

$$\Omega (R,z)={\Omega }_{\mathrm{K}}{\left[1-\frac{3}{2}\frac{z^2}{R^2}+\frac{P}{\rho {\left(R{\Omega }_{\mathrm{K}}\right)}^2}\left(\gamma \frac{\partial \ln {\rho }_{\mathrm{g}}}{\partial \ln R}-q\right)\right]}^{1/2}.$$

$$R\frac{\partial {\Omega }^2}{\partial z}=\frac{c_s^2(R)}{{\left(1+ϵ\right)}^2}\{\frac{\partial ϵ}{\partial R}\frac{\partial \ln {\rho }_{\mathrm{g}}}{\partial z}-\frac{\partial ϵ}{\partial z}\frac{\partial \ln P}{\partial R}$$

$$-\frac{q}{R}\left(1+ϵ\right)\frac{\partial \ln {\rho }_{\mathrm{g}}}{\partial z}\}$$

$$N_z^2\frac{\text{buildrel}>}{\sim }{\left|R\frac{\partial {\Omega }^2}{\partial z}\right|}_{\mathrm{temp}}\equiv \frac{|q|c_s^2}{(1+ϵ)R}\left|\frac{\partial \ln {\rho }_{\mathrm{g}}}{\partial z}\right|$$

$$\frac{1}{\left(1+ϵ\right)}\left|\frac{\partial ϵ}{\partial z}\right|\frac{\text{buildrel}>}{\sim }\frac{|q|}{R}.$$

$${ϵ}_0\frac{\text{buildrel}>}{\sim }{h}_{\mathrm{g}}|q|\frac{H_{\mathrm{d}}}{H_{\mathrm{g}}},$$

$$T_{\mathrm{stop}}\simeq \frac{10^{-3}}{F}\left(\frac{s}{100\mu \mathrm{m}}\right)\left(\frac{{\rho }_{\bullet }}{\mathrm{g}{\mathrm{cm}}^{-3}}\right){\left(\frac{R_0}{20\mathrm{au}}\right)}^{3/2},$$

Full text

1 Introduction

In the standard ‘core accretion’ scenario, planets are formed by the aggregation of km-sized or larger planetesimals via their mutual gravitational attraction (Safronov, 1969; Lissauer, 1993; Raymond et al., 2014; Helled et al., 2014). The precursor to planet formation is therefore planetesimal formation (Chiang & Youdin, 2010). However, solids in protoplanetary discs (PPDs) begin as micron-sized dust. These small grains can grow to mm-cm sizes by sticking, but collisional growth beyond this size is limited by bouncing or fragmentation (Blum, 2018). Other effects, such as turbulence or radial drift, can also hinder particle growth (Johansen et al., 2014).

One way to circumvent these growth barriers is to invoke the collective gravitational collapse of a population of dust grains (Goldreich & Ward, 1973). A dynamical gravitational instability (GI) in dust requires a high dust-to-gas ratios, $\rho_{\mathrm{d}}/\rho_{\mathrm{g}}\gg 1$, where $\rho_{\mathrm{d}}$ and $\rho_{\mathrm{g}}$ are the dust and gas mass densities, respectively (Chiang & Youdin, 2010; Shi & Chiang, 2013). However, in a newly formed PPD, where dust and gas are well-mixed, we expect $\rho_{\mathrm{d}}/\rho_{\mathrm{g}}\sim 10^{-2}$ everywhere, as typical of the dust-to-gas mass ratio in the interstellar medium. Thus, to trigger GI for planetesimal formation, dust grains must first be concentrated to sufficiently high volume densities relative to the gas.

A promising mechanism for concentrating dust is the streaming instability (SI, Youdin & Goodman, 2005; Youdin & Johansen, 2007). The SI naturally arises from the mutual dust-gas drag in rotating discs (Jacquet et al., 2011), and is driven by the relative drift between dust and gas caused by radial pressure gradients. More recent physical interpretations of the SI have been given by Lin & Youdin (2017, hereafter LY17) and Squire & Hopkins (2018).

In principle, the SI can occur at any dust-to-gas ratios, but for small particles the SI grows slowly under dust-poor conditions (\rho_{\mathrm{d}}\;\raisebox{-3.44444pt}{\buildrel\textstyle<\over{\sim}}\;\rho_{\mathrm{g}}). Furthermore, in this limit the SI does not lead to dust clumping (Johansen & Youdin, 2007). Indeed, numerical simulations show strong clumping, and hence planetesimal formation, only when \rho_{\mathrm{d}}\;\raisebox{-3.44444pt}{\buildrel\textstyle>\over{\sim}}\;\rho_{\mathrm{g}} (e.g. Johansen et al., 2009b; Bai & Stone, 2010a, b; Yang et al., 2017; Simon et al., 2016). How to realise this optimal condition for the SI , at least locally in a PPD, is therefore a crucial issue for planet formation.

Several other dust-concentrating mechanisms may act to seed the SI (Chiang & Youdin, 2010; Johansen et al., 2014). These include dust-trapping by pressure bumps , zonal flows, and vortices (Barge & Sommeria, 1995; Johansen et al., 2009a; Dittrich et al., 2013; Lyra & Lin, 2013; Raettig et al., 2015; Pinilla & Youdin, 2017); gas removal by photoevaporation or disc winds (Alexander et al., 2014; Bai, 2016); other dust-gas drag instabilities (Takahashi & Inutsuka, 2014; Gonzalez et al., 2017); radial pile-up of dust (Kanagawa et al., 2017); and dust settling (Nakagawa et al., 1986; Dubrulle et al., 1995; Takeuchi & Lin, 2002). Some of these processes may directly induce GI in dust without going through SI.

In this work we focus on dust settling. This is perhaps the simplest way to increase $\rho_{\mathrm{d}}/\rho_{\mathrm{g}}$: dust particles sediment towards the disc midplane due to the vertical component of stellar gravity and drag forces from the gas. However, dust settling requires the disc to be sufficiently laminar (Dubrulle et al., 1995; Youdin & Lithwick, 2007). Otherwise, turbulence can stir up the particles. Turbulence may arise from the settling process itself, for example via Kelvin-Helmholtz instabilities associated with the differential rotation between the dust and gas layers (KHI, Johansen et al., 2006; Barranco, 2009; Chiang, 2008; Lee et al., 2010), or weak SI (Johansen & Youdin, 2007; Bai & Stone, 2010a). In these cases, $\rho_{\mathrm{d}}/\rho_{\mathrm{g}}$ may be self-limited and never reach unity.

Furthermore, there are external sources of turbulence in PPDs. Magnetic instabilities (Balbus & Hawley, 1991) and gaseous GI (Gammie, 2001) can have significant impact on dust dynamics (e.g. Rice et al., 2004; Fromang & Papaloizou, 2006; Balsara et al., 2009; Gibbons et al., 2012; Zhu et al., 2015; Shi et al., 2016; Riols & Lesur, 2018; Yang et al., 2018).

In addition, there is a plethora of newly (re)discovered, purely hydrodynamic instabilities that drive turbulence (Fromang & Lesur, 2017; Klahr et al., 2018; Lyra & Umurhan, 2018). These include: the zombie vortex instability (Marcus et al., 2015; Lesur & Latter, 2016; Umurhan et al., 2016), the convective overstability (Klahr & Hubbard, 2014; Lyra, 2014; Latter, 2016), and the vertical shear instability (VSI, Nelson et al. 2013; Lin & Youdin 2015, hereafter LY15; Barker & Latter 2015). These instabilities are relevant to magnetically inactive or non-self-gravitating regions of the disc. Among them, the VSI is particularly relevant to dust settling.

The VSI taps into the free energy associated with a vertical gradient in the gaseous disc’s angular velocity $\Omega$ (LY15, Barker & Latter 2015). Such a vertical shear, $\partial_{z}~{}\Omega\neq~{}0$, arises whenever the gas is baroclinic. This is the case in the outer parts of PPDs where a non-uniform radial temperature profile is set by stellar irradiation (Stoll & Kley, 2014). The VSI is characterised by large-scale vertical motions (Nelson et al., 2013) and the associated stress is much stronger than that in the radial direction (Stoll et al., 2017). It is then natural to question whether or not dust particle can still settle to the point of triggering rapid SI.

Recent numerical simulations have begun to study the influence of the VSI on particle dynamics in PPDs (Stoll & Kley, 2016; Flock et al., 2017). These studies find that small particles with sizes  $\buildrel\textstyle<\over{\sim}$ mm can be efficiently lifted by VSI turbulence. However, these studies model dust as passive particles, i.e. the back-reaction of dust-drag onto the gas is neglected. However, this feedback would induce an effective buoyancy force, as the dust adds weight to the mixture but not pressure (LY17). On the other hand, buoyancy is known to be effective in stabilising the VSI (LY15). It is thus necessary to include particle back-reaction onto the gas to study dust settling against VSI turbulence.

In this work, we study dust settling against the VSI by applying LY17’s thermodynamic model of dusty PPDs, which implicitly includes particle back-reaction. This allows us investigate the effect of the overall solid abundance on dust settling in VSI-active discs. LY17’s framework was previously applied to study dusty disc-planet interaction (Chen & Lin, 2018, hereafter CL18). Here we adapt CL18’s numerical models.

This paper is organised as follows. In §2 we briefly review the governing equations of LY17 and CL18, which forms the basis of this work. We then describe our initial disc models in §3. Our numerical setup is described in §4, along with several diagnostics. We present simulation results in §5. We find that dust settling can overcome the VSI if the local metallicity, $\Sigma_{\mathrm{d}}/\Sigma_{\mathrm{g}}$, where $\Sigma_{\mathrm{d,g}}$ are dust and gas surface densities, respectively, is a few times above the solar value. We interpret and discuss our results in §6, including its implication for planetesimal formation and observations of PPDs. We conclude and summarise in §7.

2 Basic equations

We consider an accretion disc comprised of gas and dust orbiting a central star of mass $M_{*}$. We neglect disc self-gravity, magnetic fields, and viscosity. Cylindrical $(R,\phi,z)$ and spherical $(r,\theta,\phi)$ co-ordinates are centred on the star. In the discussions below, we use subscript ‘0’ to denote evaluation at $(R,z)=(R_{0},0)$, where $R_{0}$ is a reference radius. The Keplerian frequency is $\Omega_{\mathrm{K}}(R)\equiv\sqrt{GM_{*}/R^{3}}$, where $G$ is the gravitational constant. We adopt a rotating frame with angular frequency $\Omega_{\mathrm{K0}}\hat{\bm{z}}$.

The gas phase has density, pressure, and velocity fields $(\rho_{\mathrm{g}},P,\bm{v}_{\mathrm{g}})$. We model a single species of dust particles as a second, pressureless fluid with mass density and velocity $(\rho_{\mathrm{d}},\bm{v}_{\mathrm{d}})$. The fluid approximation for dust is valid for sufficiently small particles (Jacquet et al., 2011).

We also define the total mass density

[TABLE]

the centre-of-mass velocity

[TABLE]

and an effective energy

[TABLE]

where $\gamma$ is a constant defined below and $\Phi$ is a time-independent gravitational potential ($\partial_{t}\Phi=0$). Our numerical simulations evolve the hydrodynamic set $(\rho,\bm{v},E)$, see §2.3.

2.1 Polytropic gas

We consider a locally polytropic gas with

[TABLE]

where $c_{s}(R)$ is a prescribed sound-speed profile and $\gamma$ is the constant polytropic index. We take $\gamma=1.001$ so the gas is nearly isothermal in its thermodynamic property. We set

[TABLE]

where $q$ is the constant temperature profile index. Such a temperature profile is said to be locally or vertically isothermal. We also define the pressure scale height $H_{\mathrm{g}}$

[TABLE]

and the gas disc aspect ratio is $h_{\mathrm{g}}\equiv H_{\mathrm{g}}/R$.

2.2 Tightly coupled dust

The dust fluid interacts with the gas via a drag force characterised by the particle stopping time $\tau_{\mathrm{s}}$. We consider particles in the Epstein regime with $\tau_{\mathrm{s}}=\rho_{\bullet}s/\rho_{\mathrm{g}}c_{s}$ (Weidenschilling, 1977), where $\rho_{\bullet}$, $s$ is the internal density and radius of a grain, respectively.

However, in the one-fluid framework adopted below, it is convenient to define the relative stopping time $t_{\mathrm{s}}=\tau_{\mathrm{s}}\rho_{\mathrm{g}}/\rho$ (Laibe & Price, 2014a). This is the characteristic decay timescale for the relative drift between dust and gas, $\left|\bm{v}_{\mathrm{d}}-\bm{v}_{\mathrm{g}}\right|$, due to the mutual drag forces. The Stokes number can be defined as $t_{\mathrm{s}}\Omega_{\mathrm{K}}.$

We parameterise $t_{\mathrm{s}}$ via a dimensionless stopping time $T_{\mathrm{stop}}$ such that

[TABLE]

Thus $T_{\mathrm{stop}}=\Omega_{\mathrm{K0}}\rho_{\bullet}s/\rho_{0}c_{s0}$ is also the Stokes number for a particle of a given size and internal density at the reference point. For convenience we will also refer to $T_{\mathrm{stop}}$ as the particle size.

In addition, we consider small particles with

[TABLE]

In this limit dust particles are tightly, but not necessarily perfectly, coupled to the gas; and one can relate dust and gas velocities by the terminal velocity approximation:

[TABLE]

(Youdin & Goodman, 2005; Jacquet et al., 2011; Price & Laibe, 2015, LY17). Note that a differential Coriolis force has been neglected in Eq. 9 as it is of order $t_{\mathrm{s}}^{2}$ (see, e.g. Laibe & Price, 2014b). Eq. 9 reflects the fact that particles tend to drift towards pressure maxima (Whipple, 1972; Weidenschilling, 1977).

2.3 Evolutionary equations

With the above definitions and approximations, LY17 and CL18, building upon Laibe & Price (2014a), showed that the dusty gas disc can be described by the following evolutionary equations:

[TABLE]

where the dust fraction $f_{\mathrm{d}}\equiv\rho_{\mathrm{d}}/\rho$ can be obtained from the equation of state (Eq. 4)

[TABLE]

and $\rho_{\mathrm{g}}=(1-f_{\mathrm{d}})\rho$. We define the dust-to-gas ratio $\epsilon\equiv\rho_{\mathrm{d}}/\rho_{\mathrm{g}}=f_{\mathrm{d}}/\left(1-f_{\mathrm{d}}\right)$. A similar one-fluid model was also adopted by Fromang & Papaloizou (2006) to simulate dusty, magnetised discs.

In the rotating frame the total gravitational potential is

[TABLE]

Here,

[TABLE]

is the stellar potential in the thin-disc approximation ($|z|\ll R$). We adopt this approximate form in order specify initial conditions analytically.

Eqs. 10—12 provides a hydrodynamic description of the dusty gas. The dust density is recovered from the equation of state (4), and the dust velocity from the terminal velocity approximation ($\bm{v}_{\mathrm{d}}=\bm{v}+t_{\mathrm{s}}\nabla P/\rho$). Eq. 10—11 physically corresponds to mass and momentum conservation for the mixture.

However, Eq. 12 is not a physical energy equation. It simply stems from dust mass conservation, the equation of state, and the terminal velocity approximation (LY17; CL18). The term $\propto\nabla c_{s}^{2}$ on the right-hand-side models the rapid cooling that enforces a locally (vertically) isothermal temperature profile, and is the term responsible for the VSI. The second term on the right-hand-side models dust-gas decoupling.

2.4 Dusty entropy and buoyancy

Describing the dusty gas in terms of the hydrodynamic variables $(\rho,\bm{v},P)$ allows us to define the (dimensionless) effective entropy of the system as

[TABLE]

(LY17). A non-uniform entropy profile produces buoyancy forces. For this work the relevant quantity is the squared vertical buoyancy frequency

[TABLE]

where the second equality applies to the vertically isothermal temperature profiles we consider.

The notion of a dust-induced buoyancy is most applicable to dust perfectly coupled to the gas. In this limit, dust particles add to the mixture’s inertia but not the pressure. Dust therefore ‘weighs down’ the fluid. Hence there is an associated buoyancy force. Settled dust provides a stabilising force against vertical motions of the VSI since $N_{z}^{2}>0$ in that case (see also §3.2).

However, settling in the first place requires some decoupling between dust and gas. This diminishes the effective buoyancy as dust can slip past the gas. The definition above should thus be interpreted as the maximum buoyancy induced by dust.

3 Disc models

We initialise our discs by first specifying the dust-to-gas ratio distribution

[TABLE]

Here and below subscript ‘mid’ indicates the midplane value. The profile $\epsilon_{\mathrm{mid}}(R)$ is chosen to minimise radial disc evolution, as described in Appendix B of CL18. The characteristic thickness in the dust-to-gas ratio, $H_{\epsilon}$, is defined via

[TABLE]

Here the dust layer thickness $H_{\mathrm{d}}$ is set to $0.99H_{\mathrm{g}}$ so that initially the dust is vertically well-mixed with the gas ($H_{\epsilon}\gg 1$). For $t>0$ we fit the vertical dust density distribution to a Gaussian profile to obtain $H_{\mathrm{d}}$ (see §4.2).

We scale the dust-to-gas ratio by specifying the metallicity $Z$ at the reference radius

[TABLE]

(Johansen et al., 2014). Recall $\Sigma_{\mathrm{d,g}}$ are the surface densities in dust and gas, respectively.

We obtain the disc structure by assuming dynamical equilibrium and axisymmetry:

[TABLE]

Here $\Omega\equiv v_{\phi}/R+\Omega_{\mathrm{K0}}$ is the angular velocity in the inertial frame. For simplicity we set $\gamma=1$ to integrate Eq. 22, making use of Eq. 18, to obtain

[TABLE]

We set the midplane gas density to

[TABLE]

with gas surface density

[TABLE]

and we fix $p=3/2$. The total density is $\rho~{}=~{}\rho_{\mathrm{g}}(1+~{}\epsilon)$. Finally, the disc orbital frequency is obtained from Eq. 21,

[TABLE]

3.1 Vertical shear

The disc rotation possesses vertical shear. This can be seen directly from Eq. 21—22:

[TABLE]

Vertical shear can thus arise from a non-uniform dust-to-gas ratio distributions (the first two terms), and from the imposed radial temperature gradient (the last term). However, LY17 showed that gradients in $\epsilon$ generally do not lead to axisymmetric instabilities in PPDs, unless sharp radial edges are present. For axisymmetric models, as will be adopted below, dust settling in radially smooth discs only contributes to stabilisation via vertical buoyancy forces (cf. the non-axisymmetric KHI, see Chiang, 2008). The VSI in our simulations are associated with the imposed radial temperature gradient.

3.2 Buoyant stabilisation

We expect the VSI to be affected by buoyancy forces when the buoyancy frequency becomes comparable to the vertical shear rate such that

[TABLE]

(LY15), where the right-hand-side is the last term in Eq. 27. This is the only destabilising effect in axisymmetric discs in the limit of vanishing particle size (LY17). Evaluating this inequality using Eq. 17 gives

[TABLE]

For the $\epsilon$ profiles we adopt (Eq. 18), we have $\left|\partial_{z}\epsilon\right|\sim~{}\epsilon_{\mathrm{mid}}/H_{\epsilon}$. Then, considering the fiducial radius without loss of generality, we expect dust-induced buoyancy becomes important for midplane dust-to-gas ratios

[TABLE]

for somewhat settled dust (H_{\mathrm{d}}\;\raisebox{-3.44444pt}{\buildrel\textstyle<\over{\sim}}\;H_{\mathrm{g}}), where the definition Eq. 19 was used. Thus we can expect dust-induced buoyant stabilisation of the VSI even when $\epsilon_{0}\sim h_{\mathrm{g}}|q|\ll 1$. LY17 give a similar estimate and confirm it with linear stability calculations (see their sections 6.1—6.2).

The VSI is sensitive to dust-loading because the temperature-induced, destabilising vertical shear is weak, so even a small amount of dust -induced buoyancy can stabilise it (Nelson et al. 2013; LY15; LY17). We give a numerical example of this effect in Appendix A.

4 Numerical simulations

We carry out global, axisymmetric numerical simulations using the pluto hydrodynamics code (Mignone et al., 2007, 2012), modified to evolve the dust-free model equations 10—12, as described in CL18. We configure pluto with piecewise linear reconstruction , the HLLC Riemann solver, and second order Runge-Kutta time integration. We found this default setup to be robust111Test simulations with piecewise parabolic interpolation, a Roe solver, and third order Runge-Kutta time integration produced numerical artefacts near the vertical boundaries., but its diffusive nature likely suppresses the streaming instability.

We adopt spherical co-ordinates with $r\in[0.6,1.4]R_{0}$ and a polar angle range such that $\tan{\left(\pi/2-\theta\right)}\in[-2,2]h_{\mathrm{g0}}$, i.e. two gas scale-heights above and below the disc midplane at the reference radius. We fix $h_{\mathrm{g0}}=0.05$. The standard resolution is $N_{r}\times N_{\theta}=1920\times 480$; with logarithmic and uniform spacing in $r$ and $\theta$, respectively. This corresponds to a resolution of approximately 100 cells per $H_{\mathrm{g}}$ at $R=R_{0}$. We impose reflective boundary conditions in $\theta$, and set radial boundary ghost zones to their initial values.

4.1 Physical parameters

Our fiducial parameter values are $q=1$, $T_{\mathrm{stop}}=10^{-3}$, and $Z=10^{-2}$. As shown below, VSI turbulence prevents dust settling for this setup. We are interested in how these parameters can be varied to enable settling:

We vary the imposed radial temperature gradient $q~{}\in~{}[0,1]$. Smaller $|q|$ decreases the strength of the VSI turbulence and is expected to help settling. 2. 2.

We vary the stopping time $T_{\mathrm{stop}}\in[10^{-3},10^{-2}]$. Larger $T_{\mathrm{stop}}$, corresponding to larger particles, is expected to settle faster and avoid stirring by the VSI, at least initially. 3. 3.

We vary the overall metallicity $Z\in[10^{-2},10^{-1}]$. Increasing $Z$ is expected to stabilise the VSI and favour settling.

For comparison, the oft-used Minimum Mass Solar Nebula (MMSN) has a shallower temperature profile with $q~{}=~{}0.43$ and is slightly more dusty at $Z=0.015$ (Chiang & Youdin, 2010) than our fiducial setup. Epstein drag stopping times in MMSN-like disc models vary as

[TABLE]

where $F$ is the mass scale relative to the MMSN. At fixed $T_{\mathrm{stop}}$, smaller particles reside at larger distances. In a $F=1$ disc and taking $\rho_{\bullet}=1\,\mathrm{gcm}^{-3}$, we find $T_{\mathrm{stop}}=10^{-3}$ corresponds to particle sizes from $\sim 280\mu\mathrm{m}$ at $10$au to $\sim 9\mu\mathrm{m}$ at $100$au; or between $\sim 3\mathrm{mm}$ and $90\mu\mathrm{m}$ for $T_{\mathrm{stop}}=10^{-2}$.

4.2 Diagnostics

We analyse our simulation results using the diagnostics listed below. We denote mean value of a quantity $Q$ by $\left\langle Q\right\rangle$. Unless otherwise stated, the spatial average is taken over the shell $\left|r-R_{0}\right|\leq H_{\mathrm{g0}}$. We also use $\mathrm{rms}(Q)$ or $\sqrt{\left\langle Q^{2}\right\rangle}$ to denote the root-mean-squared value of $Q$. Time-averaged values are denoted with an additional subscript $t$. Time is reported in units of the orbital period at $R_{0}$ and is denoted by $P_{0}$.

• •

The dust scale-height $H_{\mathrm{d}}$ is obtained by fitting the dust density with

[TABLE]

where $\rho_{\mathrm{d}}(R,0)$ is the mid-plane dust density. We then normalise $H_{\mathrm{d}}$ by $H_{\mathrm{g}}$.

• •

The midplane dust-to-gas ratio $\epsilon_{0}\equiv\left(\rho_{\mathrm{d}}/\rho_{\mathrm{g}}\right)_{0}$. We measure both average and maximum values within the sampling region.

• •

The midplane vertical velocity $v_{\mathrm{z0}}$.

• •

Following Stoll et al. (2017), we measure the Reynolds stresses associated with vertical motions via

[TABLE]

where $\Delta$ denotes the deviation from the initial value. We then define the normalised vertical stress

[TABLE]

as a dimensionless measure of the VSI strength, where $\rho_{\mathrm{i}}$ is the initial density field. Note that in time evolution plots below, the average is also taken over the height of the disc.

5 Results

5.1 Lifting dust particles by the VSI

We first illustrate the effect of VSI turbulence on dust settling by varying the radial temperature gradient via the power-law index $q$. Here, we fix the particle size and metallicity such that $T_{\mathrm{stop}}=10^{-3}$ and $Z=0.01$, respectively. Larger $|q|$ leads to stronger vertical shear and hence more vigorous VSI turbulence.

Fig. 1 shows the time evolution of the average dust scale-height normalised by the gas scale-height, $\left\langle H_{\mathrm{d}}/H_{\mathrm{g}}\right\rangle$; the average and maximum midplane dust-to-gas ratio, $\left\langle\epsilon_{0}\right\rangle$ and $\mathrm{max}(\epsilon_{0})$, respectively; the root-mean-squared midplane vertical velocity, $\sqrt{\left\langle v_{\mathrm{z0}}^{2}\right\rangle}$; and the dimensionless vertical stress, $\mathcal{R}_{z\phi}$; for $q=0$, $0.3$, $0.5$, and $1$. Fig. 2 compares several snapshots of the dust-to-gas ratio from these runs. It is evident that VSI turbulence hinders dust settling.

For $q=0$ the disc is strictly isothermal and does not develop the VSI. In this laminar disc, dust settling proceeds until the resolution limit ($\simeq 10^{-2}H_{\mathrm{g}}$) and $\epsilon_{0}\sim 1$. Introducing moderate VSI turbulence with $q=0.5$ prevents the dust from settling below $\sim 0.1H_{\mathrm{g}}$. Consequently $\epsilon_{0}$ only increases to $O(0.1)$.

In the fully turbulent disc with $q=1$ we find no signs of dust settling over the simulation timescale. Fig. 2 show in this case that dust remains well-mixed with the gas throughout the vertical disc column with $\rho_{\mathrm{d}}/\rho_{\mathrm{g}}\sim 10^{-2}$, with undulations near the surface. The banded structure in $\rho_{\mathrm{d}}/\rho_{\mathrm{g}}$ with characteristic separation $\sim H_{\mathrm{g}}/2$ is caused by the advection of dust particles by the large-scale up/down VSI motions , which have no vertical structure. This is shown in Fig. 3. The average vertical stress $\mathcal{R}_{z\phi}\simeq 5\times 10^{-3}$ is consistent with Stoll et al. (2017), who found \mathcal{R}_{z\phi}\;\raisebox{-3.44444pt}{\buildrel\textstyle<\over{\sim}}\;0.01 for $|z|\leq 2H_{\mathrm{g}}$.

In Fig. 4 we compare the time-averaged dust distribution and turbulence properties from a suite of runs with $q\in[0.1,1]$. The averaging period is taken between when the system enters a quasi-steady state and the end of the simulation. We omit the $q=0$ run here since in that case there is no physical steady state. We find $\left\langle H_{\mathrm{d}}/H_{\mathrm{g}}\right\rangle_{t}$ increases rapidly with $q$, but for q\;\raisebox{-3.44444pt}{\buildrel\textstyle>\over{\sim}}\;0.5 the dust becomes well-mixed with $\left\langle H_{\mathrm{d}}/H_{\mathrm{g}}\right\rangle_{t}\simeq 1$ and $\left\langle\epsilon_{0}\right\rangle_{t}\simeq 0.01$ ($=Z$ here). Notice the well-mixed and settled cases are separated by $\sqrt{\left\langle v_{\mathrm{z0}}^{2}\right\rangle}_{t}\sim\sqrt{T_{\mathrm{stop}}}c_{s0}$ or $\mathcal{R}_{z\phi\mathrm{,}t}\sim T_{\mathrm{stop}}$. The former dependence is consistent with theoretical dust settling models, described below. As $\left\langle v_{\mathrm{z0}}^{2}\right\rangle_{t}$ and $\mathcal{R}_{z\phi\mathrm{,}t}$ increases with $|q|$, this suggest that given a fixed particle size there exists a critical temperature gradient, $|q|_{*}$, below which dust-settling is allowed. However, as shown below, vertical stresses and velocities also depends on particle size and metallicity.

For comparison, we also plot in Fig. 4 the expected dust scale-height based on Dubrulle et al. (1995)’s model,

[TABLE]

(see also Youdin & Lithwick, 2007; Zhu et al., 2015), where the particle stopping time $\tau_{\mathrm{s}}$ is measured at the disc midplane, and we regard the eddy timescale $t_{\mathrm{eddy}}$ as a model parameter for simplicity. Setting $\Omega_{\mathrm{K}}t_{\mathrm{eddy}}=0.15$ gives a good match to the measured $\left\langle H_{\mathrm{d}}/H_{\mathrm{g}}\right\rangle_{t}$ in the settled cases, while it somewhat underestimates the well-mixed cases. This value of $t_{\mathrm{eddy}}$ is similar to that found by Stoll & Kley (2016).

5.2 Effect of particle size

We now examine the effect of particle size, which is parameterised by $T_{\mathrm{stop}}$ in our simulations. Here we fix $q=1$ and $Z=0.01$.

Fig. 5 compares the fiducial case with $T_{\mathrm{stop}}=10^{-3}$ to runs with $3$, $5$, and $10$ times that particle size. Increasing $T_{\mathrm{stop}}$ allows faster settling, since the settling speed $\left|v_{z,d}\right|\sim c_{s}T_{\mathrm{stop}}$ (Takeuchi & Lin, 2002). Larger particles can settle to a thinner layer before VSI-stirring, which begins at $t\sim 100P_{0}$ independently of $T_{\mathrm{stop}}$. That is, the initial VSI growth rate is insensitive to particle size. However, VSI-stirring weakens with increasing $T_{\mathrm{stop}}$, as $\left\langle H_{\mathrm{d}}/H_{\mathrm{g}}\right\rangle$ ‘bounce’ to smaller values with larger particles after the VSI activates. This is also reflected in the much smaller $\mathcal{R}_{z\phi}$ with $T_{\mathrm{stop}}=0.005$ and $0.01$ compared to cases with $T_{\mathrm{stop}}=0.001$ and $0.003$.

In the $T_{\mathrm{stop}}=0.003$ run, dust is almost perfectly remixed by the initial VSI. A quasi-steady state is apparently attained, but only until $t\sim 800P_{0}$, when we observe spontaneous settling. This is shown in Fig. 6, suggesting that the prior steady state was unstable. The $T_{\mathrm{stop}}=0.005$ and $T_{\mathrm{stop}}=0.01$ runs behave similarly. In both cases dust begins to settle soon after being remixed by the VSI. However, they also exhibit hints of spontaneous settling. The $T_{\mathrm{stop}}=0.005$ run attains a brief steady state (150P_{0}\;\raisebox{-3.44444pt}{\buildrel\textstyle<\over{\sim}}\;t\;\raisebox{-3.44444pt}{\buildrel\textstyle<\over{\sim}}\;250P_{0}) before settling; while at the end of the $T_{\mathrm{stop}}=0.01$ run $\left\langle\epsilon_{0}\right\rangle$ suddenly increases. We discuss in §6.2 how spontaneous settling may be possible due to particle back-reaction.

Fig. 7 shows the final dust-to-gas ratio and vertical velocity distribution in the $T_{\mathrm{stop}}=0.01$ run. A dust clump develops at the outer boundary, possibly due to boundary conditions, which may be related to the increased dust mixing in R\;\raisebox{-3.44444pt}{\buildrel\textstyle>\over{\sim}}\;1.1R_{0}. However, the sampling region around $R\simeq R_{0}$ is unaffected, and particles there have settled to a thin layer with little vertical motion, due to self-stabilisation by dust-induced buoyancy (LY17).

Fig. 8 compares an extended set of simulations with varying $T_{\mathrm{stop}}$. Here the time average is taken between the onset of the VSI and the end of the simulations. Significant settling occurs for T_{\mathrm{stop}}\;\raisebox{-3.44444pt}{\buildrel\textstyle>\over{\sim}}\;0.005, for which $\left\langle H_{\mathrm{d}}/H_{\mathrm{g}}\right\rangle_{t}$ drops to $\sim 0.1$ and $\left\langle\epsilon_{0}\right\rangle_{t}\sim O(0.1)$.

The measured $\left\langle H_{\mathrm{d}}/H_{\mathrm{g}}\right\rangle_{t}$ is consistent with the model described by Eq. 35 if we choose $\Omega_{\mathrm{K}}t_{\mathrm{eddy}}=0.2$, again similar to eddy timescales found by Stoll & Kley (2016).

Although less clear-cut, this figure also indicates that dust settling occurs when \sqrt{\left\langle v_{\mathrm{z0}}^{2}\right\rangle}_{t}\;\raisebox{-3.44444pt}{\buildrel\textstyle<\over{\sim}}\;\sqrt{T_{\mathrm{stop}}}c_{s0} or \mathcal{R}_{z\phi\mathrm{,}t}\;\raisebox{-3.44444pt}{\buildrel\textstyle<\over{\sim}}\;T_{\mathrm{stop}}.

5.3 Effect of dust abundance

Dust-induced buoyancy is expected to help settling by stabilising the VSI (LY17). Here we investigate this effect more precisely by varying the solid abundance or metallicity $Z=\Sigma_{\mathrm{d}}/\Sigma_{\mathrm{g}}$. Larger $Z$ corresponds to stronger dust feedback onto the gas. For our fiducial parameter values $(q,T_{\mathrm{stop}},Z)=(1,10^{-3},10^{-2})$ dust does not settle. We thus consider $Z\in[0.01,0.1]$ in this section. In Appendix A we demonstrate the VSI is particularly sensitive to dust-loading in weakly turbulent discs, in which case the stabilising effect of feedback manifests even with $\epsilon_{0}$ of a few per cent.

Fig. 9 compares the disc evolution for $Z=0.01$, $0.03$, $0.05$, and $0.1$. Increasing $Z$ does not affect the settling timescale. Instead, larger $Z$ delays the VSI, which allows the dust to settle to a thinner layer. This reflects the stabilising effect of particle back-reaction. For $Z=0.03$, the VSI remixes the dust at $t\sim 200P_{0}$, and settles slowly afterwards. For $Z=0.05$ and $Z=0.1$, the initial VSI only halts the settling process, and does not increase $\left\langle H_{\mathrm{d}}/H_{\mathrm{g}}\right\rangle$, unlike the $Z=0.01$ and $Z=0.03$ runs.

For $Z=0.1$, we find $\left\langle\epsilon_{0}\right\rangle\sim 2$ for t\;\raisebox{-3.44444pt}{\buildrel\textstyle>\over{\sim}}\;300P_{0}. This is expected to develop the SI, but our numerical setup and resolution is insufficient to resolve the SI for the small particles considered here (Yang & Johansen, 2016). In the $Z=0.05$ and $Z=0.03$ runs, $\left\langle\epsilon_{0}\right\rangle$ does not reach unity within the simulation timescale, but is still significantly enhanced relative to the $Z=0.01$ run. Comparing these runs to the previous simulations with increased $T_{\mathrm{stop}}$ (Fig. 5) show that increasing $Z$ is more effective in raising the dust-to-gas ratio than increasing particle size. This is not surprising since, for a given $H_{\mathrm{d}}$ a larger $Z$ implies a larger $\epsilon_{0}$ (Eq. 20); let alone in the present cases where a smaller dust scale-height is reached with increasing $Z$ compared to the same increase in $T_{\mathrm{stop}}$.

The $Z=0.05$ and $Z=0.1$ cases show much reduced midplane vertical velocities and vertical turbulence compared to $Z=0.01$, which is consistent with the former cases reaching much smaller values of $\left\langle H_{\mathrm{d}}/H_{\mathrm{g}}\right\rangle$. The $Z=0.1$ run also displays a slow drop in $\mathcal{R}_{z\phi}$ for t\;\raisebox{-3.44444pt}{\buildrel\textstyle>\over{\sim}}\;400P_{0}, with a corresponding rise in $\left\langle\epsilon_{0}\right\rangle$; while $\left\langle H_{\mathrm{d}}/H_{\mathrm{g}}\right\rangle$ remains relatively constant. This signifies particle settling into midplane unhindered; in fact this settling is stabilising the system. The $Z=0.03$ run attains $\mathcal{R}_{z\phi}\simeq 5\times 10^{-3}$ early on, similar to the $Z=0.01$ case, but it eventually drops to $\mathcal{R}_{z\phi}\simeq 2\times 10^{-3}$. This also reflects the slow settling observed after the initial VSI in that case.

Fig. 10 shows $\rho_{\mathrm{d}}/\rho_{\mathrm{g}}$ and $v_{z}$ at the end of the $Z=0.1$ simulation. Maximum vertical velocities are reduced by at a factor of two from the standard case with $Z=0.01$ (Fig. 3), while midplane values are reduced by an order of magnitude (Fig. 9). Similar to the large-particle simulation (Fig. 7), the dust midplane is quiescent, although not completely inactive. This dense dust layer also decouples the upper, dust-free region from that in the lower disc.

Fig. 11 shows the corresponding vertical stress and vertical velocity profiles. Vertical stresses generally increase away from the midplane, reaching $\mathcal{R}_{z\phi}\sim 10^{-3}$ near the boundaries . This reflects VSI activity in the dust-poor layers above and below the midplane. The spike at $z=0$ is associated with the fact that the dust layer rotates closer to the Keplerian speed compared to the initially gas-dominated state. (Recall from Eq. 33 that $R_{z\phi}\propto\Delta v_{\phi}$.) Although vertical velocities are reduced towards the midplane, notice that $v_{\mathrm{z0}}\neq 0$. Such residual midplane vertical motions are

due to disturbance by bulk vertical motions from the VSI-active layers (the midplane dust layer in itself is expected to be stable, see below). This explains why the dust layer maintains a finite thickness and does not settle further.

5.3.1 Dust-induced buoyancy

Fig. 12 compares the vertical buoyancy frequency and vertical shear profiles at the end of the $Z=0.1$ simulation. We show both the actual vertical shear rate, $\left|R\partial_{z}\Omega^{2}\right|$; and that associated with the imposed temperature gradient , $\left|R\partial_{z}\Omega^{2}\right|_{\mathrm{temp}}$, as estimated from Eq. 28. This quantity is the driver of the underlying VSI (LY17) and is essentially constant in time. On the other hand, the total vertical shear displayed in Fig. 12 is either dominated by the settled dust in the midplane, or by VSI modes away from it.

The total vertical shear exceeds buoyancy throughout the disc column. However, recall in axisymmetric discs that dust-induced vertical shear does not cause instabilities (LY17). Thus in our simulations instability is associated with vertical shear induced by the imposed temperature gradient, $\left|R\partial_{z}\Omega^{2}\right|_{\mathrm{temp}}$. Away from the midplane at |z|\;\raisebox{-3.44444pt}{\buildrel\textstyle>\over{\sim}}\;0.1H_{\mathrm{g}}, buoyancy frequencies vanish, whereas the temperature-related, destabilising vertical shear rate is larger there. About the midplane ($|z|\sim 0.1H_{\mathrm{g}}$), however, this destabilising effect is overwhelmed by the stabilising buoyancy. This explains the ‘layered’ structure seen in Fig. 10: a quiet midplane with VSI-active layers above and below it.

Fig. 13 compares the stabilising buoyancy to the destabilising temperature-related vertical shear, averaged over the dust layer, for several cases. We omit the $Z=0.01$ run here since $N_{z}^{2}\sim 0$ in that case. For $Z=0.03$ this ratio is \;\raisebox{-3.44444pt}{\buildrel\textstyle<\over{\sim}}\;1, though rising later on, which is consistent with the slow settling. For $Z=0.05$ and $Z=0.1$ this ratio is $\gg 1$, and the VSI is strongly stabilised by buoyancy within the dust layer.

5.3.2 Dust settling by dust-loading

We now examine the critical metallicity for dust settling to overcome the VSI. Fig. 14 compare quasi-steady states as a function of $Z$. We also plot several cases with $q=0.5$.

For $q=1$ we find stabilisation by dust-loading is most dramatic for $Z\simeq 0.03$—$0.05$; while for $Z>0.05$ dust settles to \left\langle H_{\mathrm{d}}/H_{\mathrm{g}}\right\rangle_{t}\;\raisebox{-3.44444pt}{\buildrel\textstyle<\over{\sim}}\;0.1 with little variation with respect to $Z$. On the other hand, for $q=0.5$ dust can settle to \left\langle H_{\mathrm{d}}/H_{\mathrm{g}}\right\rangle_{t}\;\raisebox{-3.44444pt}{\buildrel\textstyle<\over{\sim}}\;0.1 with Z\;\raisebox{-3.44444pt}{\buildrel\textstyle>\over{\sim}}\;0.02; suggesting the critical $Z$ for settling scales as $\sim 1/|q|$. We find the settled dust scale-heights can be fitted with Eq. 35, but this required $\Omega_{\mathrm{K}}t_{\mathrm{eddy}}\simeq 0.04$, which is smaller than that found by Stoll & Kley (2016). This indicates that increasing the metallicity reduces the turbulence eddy timescales. We find $\sqrt{\left\langle v_{\mathrm{z0}}^{2}\right\rangle}_{t}\sim\sqrt{T_{\mathrm{stop}}}c_{s0}$ continues to distinguish between settled and well-mixed cases as we vary the metallicity.

In the settled cases, we find the dust layer is too thin to affect the average vertical stress $\mathcal{R}_{z\phi\mathrm{,}t}$, i.e. the measured vertical turbulence is associated with the gas layers above and below the midplane, hence $\mathcal{R}_{z\phi\mathrm{,}t}$ flattens off. Nevertheless, we find settled cases still obey \mathcal{R}_{z\phi\mathrm{,}t}\;\raisebox{-3.44444pt}{\buildrel\textstyle<\over{\sim}}\;T_{\mathrm{stop}}.

These results indicate that in weakly VSI-turbulent discs, dust settling is sensitive to the local metallicity. This implies that small changes in $Z$, for example due to radial drift and pile-up of dust, could have significant impact on dust settling in disc regions with shallow temperature profiles (see §6 for a discussion).

6 Discussion

6.1 Overcoming VSI turbulence through particle back-reaction

Our simulations indicate that particle back-reaction onto the gas is important for dust settling in discs prone to the VSI. This implies that previous studies using passive particles may have overestimated the minimum particle sizes that can settle in VSI-active discs (Stoll & Kley, 2016; Flock et al., 2017).

We find that raising the metallicity $Z=\Sigma_{\mathrm{d}}/\Sigma_{\mathrm{g}}$, which increases the particle back-reaction, helps dust settling. In this case the settling speed is unaffected. Instead, VSI stirring is increasingly delayed due to stabilisation by dust-induced buoyancy forces. This results in thinner dust layers before the VSI halts the settling process, which are even more difficult to stir up. (The same effect occurs for increasing particle size.)

For $T_{\mathrm{stop}}=10^{-3}$ and $q=1$ we find Z\;\raisebox{-3.44444pt}{\buildrel\textstyle>\over{\sim}}\;0.05 is needed for dust to settle down to H_{\mathrm{d}}\;\raisebox{-3.44444pt}{\buildrel\textstyle<\over{\sim}}\;0.1H_{\mathrm{g}}; while in less turbulent discs with $q=0.5$ this only required Z\;\raisebox{-3.44444pt}{\buildrel\textstyle>\over{\sim}}\;0.02, slightly above the solar value. This suggests that in weakly turbulent discs the VSI is sensitive to the disc metallicity. This is consistent with the argument presented in §3.2 and the simulation in Appendix A, which show that the VSI can be stabilised when the dust-to-gas ratio $\epsilon$ becomes comparable to the disc’s aspect-ratio $h_{\mathrm{g}}$, i.e. a few per cent in PPDs.

Our simulations show that $H_{\mathrm{d}}$ generally depends on $q$, $T_{\mathrm{stop}}$, and $Z$. All three parameters affect the VSI turbulence strength $\mathcal{R}_{z\phi}$ (the dimensionless vertical Reynolds stress) and the midplane vertical velocity $v_{\mathrm{z0}}$. We find that the settled and well-mixed cases are separated by $\mathcal{R}_{z\phi}\sim T_{\mathrm{stop}}$ or $\sqrt{\left\langle v_{\mathrm{z0}}^{2}\right\rangle}\sim\sqrt{T_{\mathrm{stop}}}c_{s}$. The latter condition is consistent with the settling model of Dubrulle et al. (1995), provided an appropriate eddy timescale is chosen (see Eq. 35).

We can apply the above results to assess whether particles of a given $T_{\mathrm{stop}}$ ($>0$) and abundance $Z$ can settle. Recent analyses indicate $\left|v_{z}\right|\sim h_{\mathrm{g}}|q|c_{s}$ for VSI turbulence (Barker & Latter, 2015; Latter & Papaloizou, 2018), which is consistent with our simulations. We thus expect particles with T_{\mathrm{stop}}\;\raisebox{-3.44444pt}{\buildrel\textstyle>\over{\sim}}\;h_{\mathrm{g}}^{2}q^{2} to settle, even in the limit $Z\to 0$. However, if T_{\mathrm{stop}}\;\raisebox{-3.44444pt}{\buildrel\textstyle<\over{\sim}}\;h_{\mathrm{g}}^{2}q^{2}, then Z\;\raisebox{-3.44444pt}{\buildrel\textstyle>\over{\sim}}\;h_{\mathrm{g}}|q| is needed to weaken the VSI and enable dust settling (see §3.2).

6.2 Self-sustained settling

We propose particle back-reaction can lead to self-sustained dust settling, as illustrated in Fig. 15. In this picture, dust settles if it can induce an effective buoyancy that stabilises the VSI, which favours further settling, and the cycle repeats. This feedback cycle may be initiated by increasing $T_{\mathrm{stop}}$ via particle growth, raising the local metallicity $Z$ via radial dust drift, or reducing the radial temperature gradient $|q|$. The last trigger is less likely in PPDs at large radii because the local temperature profile is fixed by stellar irradiation.

We suggest this positive feedback loop is responsible for spontaneous settling observed in some simulations (e.g. Fig. 6). The quasi-steady balance between VSI-stirring and dust settling is delicate because the VSI is sensitive to buoyancy. If the system is perturbed such that $Z$ increases locally, then that region may begin to settle and enter the feedback loop.

However, self-sustained settling does not continue indefinitely. As dust settles to the midplane, it leaves behind dust-free, VSI-turbulent layers above and below it. Vertical motion from the active layers continue to disturb/distort the midplane dust layer, eventually halting further sedimentation.

Indeed, we find that although $\left|v_{z}\right|$ drops significantly with decreasing $|z|$, it is still non-zero about the dusty midplane. Such a layered structure in $\left|v_{z}\right|$ is characteristic of settled dust in VSI-prone discs. In fact, a similarly layered structure is found in the linear VSI with a dusty midplane (LY17, see their section 6.3). This picture is analogous to layered accretion discussed by Gammie (1996); although in that case it is the radial mass accretion rate that is minimised in the magnetic dead zones near the midplane. Dust settling in the dead zone is halted by residual activity induced by the turbulent active layers (Fromang & Papaloizou, 2006). We suggest a similar effect is present in our simulations: the VSI active layers induce residual vertical motions at the disc midplane to maintain a finite dust layer thickness.

In reality, self-sustained settling may also be terminated by other processes, such as the KHI (Chiang, 2008; Lee et al., 2010), GI (Goldreich & Ward, 1973), or SI (Youdin & Goodman, 2005). However, these effects are not captured by our numerical models.

6.3 Implications for planetesimal formation

Efficient planetesimal formation via rapid SIs require a midplane dust-to-gas ratio of order unity or larger (Youdin & Goodman, 2005; Johansen et al., 2009b).

However, our results indicate small particles would remain lofted by VSI turbulence. Here, what constitutes as ‘small’ depends on the turbulence strength and local metallicity. In the outer parts of PPDs, VSI turbulence is primarily related to the local temperature gradient $|q|$. If we take $|q|\simeq 1$ and solar metallicity $Z=0.01$, then by small we mean particles with $T_{\mathrm{stop}}$ of $O(10^{-3})$. At a distance of $20$au, this roughly corresponds to particle sizes of $O(10^{2}\mu\rm{m})$, assuming an internal density $\rho_{\bullet}=1\,\rm{g}\rm{cm}^{-3}$ (see Eq. 31).

In the above example, particles smaller than $\sim 100\mu\rm{m}$ will remain lofted. However, if they can grow to $\sim$mm sizes amidst the turbulence, e.g. by sticking (Blum, 2018), or if the local metallicity increases (e.g. due to a radial pile-up of dust, Kanagawa et al., 2017); then these particles can settle against the VSI and eventually undergo SI.

Raising the local metallicity is an effective way to achieve a local dust-to-gas ratio $\rho_{\mathrm{d}}/\rho_{\mathrm{g}}\sim 1$ in two regards. First, a larger $Z$ allows $\epsilon_{0}\sim 1$ to be attained at larger $H_{\mathrm{d}}$ (see Eq. 20). That is, the required degree of settling is smaller. Second, increasing $Z$ actively stabilises the disc against the VSI, thereby allow for effective settling.

Coincidentally, non-linear numerical simulations of the SI indicate Z\;\raisebox{-3.44444pt}{\buildrel\textstyle>\over{\sim}}\;0.02 is required for strong clumping of particles with $T_{\mathrm{stop}}\sim 0.01$—$0.1$ (Johansen et al., 2009b; Bai & Stone, 2010a) ; while smaller particles with $T_{\mathrm{stop}}\sim 10^{-3}$ require Z\;\raisebox{-3.44444pt}{\buildrel\textstyle>\over{\sim}}\;0.04 (Yang et al., 2017). These particle clumps can then undergo GI. One also requires Z\;\raisebox{-3.44444pt}{\buildrel\textstyle>\over{\sim}}\;0.03 in order to overcome the KHI associated with dust-induced vertical shear, which would overturn the dust layer (Chiang, 2008; Lee et al., 2010).

We conclude that in the outer parts of protoplanetary discs, small particles must have super-solar metallicities in order to form planetesimals self-consistently : dust can then settle against VSI turbulence, trigger efficient streaming instability, clump, and finally undergo gravitational collapse.

6.4 Application to observations

Our simulations may aid the interpretation of PPD observations. As a concrete example, we consider the PPD around HL Tau (ALMA Partnership et al., 2015). This inclined disc displays well-defined dust rings at all azimuth, which implies that the dust has settled.

For the HL Tau disc, Pinte et al. (2016) estimates $H_{\mathrm{d}}\sim 0.07H_{\mathrm{g}}$ at $100$au, and $\epsilon_{0}\sim 0.2$ assuming $Z=0.01$. They infer a dimensionless vertical diffusion parameter $\alpha\simeq 3\times 10^{-4}$, assuming passive particles (Fromang & Nelson, 2009). Interestingly, these values are consistent with our $T_{\mathrm{stop}}=0.01$ simulation222Note that we used $(q,h_{\mathrm{g}})=(1,0.05)$ while Pinte et al. and Jin et al. assumed $(q,h_{\mathrm{g}})=(0.7,0.1)$ and $(0.43,0.096)$, respectively. However, the destabilising effect of increasing $|q|$ offsets the stabilising effect of decreasing $h_{\mathrm{g}}$ (Nelson et al., 2013). shown in Fig. 8, if we identify $\alpha$ as our dimensionless vertical Reynolds stress. Using Eq. 31, $T_{\mathrm{stop}}=0.01$ corresponds to particle sizes $\sim 90\mu$m at $100$au, similar to that adopted in the HL Tau disc model of Jin et al. (2016). For mm-sized particles, $T_{\mathrm{stop}}=0.01$ requires a disc $\sim 10$ times more massive than the MMSN; or grains with internal density $\sim 0.1\mathrm{gcm}^{-3}$.

For dusty VSI turbulence we find $H_{\mathrm{d}}=H_{\mathrm{d}}(q,T_{\mathrm{stop}},Z)$. This may be used to estimate the local metallicity. For example, let us assume $q\simeq 0.5$ and particle sizes of $\sim 100\mu\rm{m}$, and a massive (young) disc with $F=10$. This gives $T_{\mathrm{stop}}\simeq 10^{-3}$ at $100$au. Then Fig. 14 suggest Pinte et al.’s estimate, $H_{\mathrm{d}}\simeq 0.07H_{\mathrm{g}}$, can be achieved with $Z\simeq 0.03$. On the other hand, if $Z\sim 10^{-3}$ so that particle back-reaction is negligible, then VSI turbulence would be too strong to allow dust to settle (see, e.g. Flock et al., 2017).

The interplay between VSI turbulence, dust settling, and the local metallicity may lead to a radially-dependent dust scale-height, as follows. Suppose dust drifts radially, say from $R_{2}$ to $R_{1}$ such that $Z(R_{1})>Z(R_{2})$. Then the increased (decreased) metallicity at $R_{1}$ ($R_{2}$) would stabilise (destabilise) the disc there. Now, if the gas density contrast between $R_{1}$ and $R_{2}$ is smaller than that in $Z$, i.e. $Z_{1}/Z_{2}>\Sigma_{\mathrm{g1}}/\Sigma_{\mathrm{g2}}$, so that the contrast in stopping times is outweighed by that in metallicity, then dust settling (mixing) should be promoted at $R_{1}$ ($R_{2}$). Ultimately, $H_{\mathrm{d}}(R_{1})<H_{\mathrm{d}}(R_{2})$ , assuming similar gas scale heights if $R_{1}>R_{2}$. That is, dust rings should have smaller scale-heights compared to dust gaps333 However, if the gas density in the dust ring is also increased significantly, then the decreased stopping times may negate the positive effect of increased metallicity on settling.. This can be tested with future observations of the radial and vertical dust structure in PPDs.

6.5 Caveats and outlooks

Below we discuss several caveats of our modelling and directions for future work.

We assumed axisymmetry throughout. This simplification overestimates the strength of VSI turbulence (Manger & Klahr, 2018). For example, Stoll & Kley (2016) found that passive particles with $T_{\mathrm{stop}}=0.002$ can settle to $H_{\mathrm{d}}\simeq 0.66H_{\mathrm{g}}$; whereas we find such particles remain perfectly mixed with the gas. On the other hand, Stoll & Kley found particles with $T_{\mathrm{stop}}=0.02$ only settles to $H_{\mathrm{d}}\simeq 0.4H_{\mathrm{g}}$, whereas we find $T_{\mathrm{stop}}=0.01$ particles can already settle to H_{\mathrm{d}}\;\raisebox{-3.44444pt}{\buildrel\textstyle<\over{\sim}}\;0.1H_{\mathrm{g}}, which could be due to the stabilising effect of particle back-reaction in our case.

Axisymmetric models also suppress vortex formation due to the VSI (Richard et al., 2016; Manger & Klahr, 2018), as well as the KHI due to the vertical shear between the dust and gas layers (Chiang, 2008). Vortices are known to act as effective dust traps (Lyra & Lin, 2013), although they can also source hydrodynamic turbulence (Lesur & Papaloizou, 2009, 2010). The KHI may become relevant if the dust layer becomes sufficiently thin and/or if the overall metallicity is large (see Appendix B). It will be necessary to perform full 3D simulations to investigate how these additional processes affect dust settling.

We adopted initial disc models that minimised radial disc evolution in order to focus on vertical dust settling (§2, see also CL18). However, in PPDs dust typically drifts inward. This may be important for dust settling because radial pile-up of dust, should it occur, would stabilise the VSI and induce dust settling. Future work should explore different initial dust and gas density profiles.

The dust-gas model of LY17 only allows for a single species of dust particles. In reality, there will be a distribution of particle sizes, each with a different metallicity. Particles that are larger or more abundant will settle first, which is expected to stabilise the system. The settling of smaller particles may then ensue. To investigate this scenario, one would need a dust-gas framework that allows for multiples species of dust grains.

7 Summary

In this paper we study the sedimentation of small dust particles in the outer parts of protoplanetary discs subject to hydrodynamic turbulence driven by the vertical shear instability (VSI). We perform 2D, axisymmetric numerical simulations to investigate the effect of VSI turbulence strength, $\mathcal{R}_{z\phi}$ (essentially the vertical Reynolds stress-to-pressure ratio); particle size or Stokes number, $T_{\mathrm{stop}}$; and solid abundance (metallicity) $\Sigma_{\mathrm{d}}/\Sigma_{\mathrm{g}}$ on dust settling.

We model VSI turbulence explicitly by imposing a fixed radial temperature gradient $q$. However, we find the resulting stress $\mathcal{R}_{z\phi}$ can also depend on $T_{\mathrm{stop}}$ and $\Sigma_{\mathrm{d}}/\Sigma_{\mathrm{g}}$. Nevertheless, our results indicate that dust settles when the resulting \mathcal{R}_{z\phi}\;\raisebox{-3.44444pt}{\buildrel\textstyle<\over{\sim}}\;T_{\mathrm{stop}} . In our simulations this is equivalent to having a midplane vertical velocity dispersion \sqrt{\left\langle v_{\mathrm{z0}}^{2}\right\rangle}\;\raisebox{-3.44444pt}{\buildrel\textstyle<\over{\sim}}\;\sqrt{T_{\mathrm{stop}}}c_{s}. The latter condition is consistent with the classic dust settling model of Dubrulle et al. (1995).

We find that dust settling is favoured by

1. weaker VSI turbulence, corresponding to shallower radial temperature profiles;

2. larger particles; and

3. increased local metallicity.

We expect dust can settle if T_{\mathrm{stop}}\;\raisebox{-3.44444pt}{\buildrel\textstyle>\over{\sim}}\;h_{\mathrm{g}}^{2}q^{2} and/or \Sigma_{\mathrm{d}}\;\raisebox{-3.44444pt}{\buildrel\textstyle>\over{\sim}}\;h_{\mathrm{g}}|q|\Sigma_{\mathrm{g}}, where $h_{\mathrm{g}}$ is the disc aspect-ratio.

Our simulations show that particle back-reaction plays an important role for dust settling against the VSI. This is because dust-loading, specifically vertical gradients in the dust-to-gas ratio, introduces a buoyancy force (LY17), which is known to be effective in stabilising the VSI (LY15). Thus settled dust is strongly self-stabilised against the VSI. This results in a quiescent dusty midplane with residual vertical velocities induced by VSI turbulence from the gaseous layers above and below it, which maintains the dust layer thickness.

We find dust settling can overcome the VSI by increasing the local metallicity. For particles with $T_{\mathrm{stop}}=10^{-3}$, corresponding to grain sizes of $\sim 100\mu$m at $20$au in MMSN-like discs, we find that \Sigma_{\mathrm{d}}\;\raisebox{-3.44444pt}{\buildrel\textstyle>\over{\sim}}\;0.04\Sigma_{\mathrm{g}} enables dust reach $\rho_{\mathrm{d}}\sim\rho_{\mathrm{g}}$ at the disc midplane. Such metallicities are also expected to trigger planetesimal formation via the streaming instability (Johansen et al., 2009b; Yang et al., 2017). Our results are thus consistent with the notion that planetesimals can form from small dust particles in the outer parts of protoplanetary discs, provided that the corresponding solid abundance is a few times larger than the solar value.

Acknowledgements

I thank the anonymous referee for key suggestions that lead to several improvements to this paper. I am indebted to Eugene Chiang for valuable comments on an early draft of this work, especially for motivating §3.2. I also thank Guillaume Laibe for useful discussions.

This work is supported by the Theoretical Institute for Advanced Research in Astrophysics (TIARA) based in the Academia Sinica Institute for Astronomy and Astrophysics (ASIAA) and the Ministry of Science and Education (MOST) grant 107-2112-M-001-043-MY3. Simulations were carried out on the TIARA High Performance Computing cluster and the TAIWANIA cluster hosted by the National Center for High-performance Computing (NCHC). I am grateful to the NCHC for computer time and facilities.

Appendix A Dust settling with negligible back-reaction

It is common to model dust in PPDs as passive particles when the dust-to-gas ratio $\epsilon\ll 1$ (e.g. Stoll & Kley, 2016; Flock et al., 2017). Here we demonstrate that this condition is problem-dependent, and it alone may be insufficient to justify neglecting particle back-reaction. Specifically, we show that feedback has significant impact in weakly VSI-turbulent discs.

We consider a disc with $q=0.3$ and particles with $T_{\mathrm{stop}}=10^{-3}$. Fig. 16 compares the simulation with $Z=10^{-2}$, as presented in the main text, to a disc with $Z=10^{-3}$. They evolve similarly until $t\sim 200P_{0}$, at which point dust in the $Z=10^{-3}$ disc undergoes complete remixing by the VSI; whereas dust continues to settle in the $Z=10^{-2}$ disc. In the latter case $\left\langle\epsilon_{0}\right\rangle\simeq 0.04\ll 1$ at $t=200P_{0}$, but this is apparently sufficient to overcome the VSI. As explained in §3.2, this is because the VSI is sensitive to buoyancy forces induced by vertical gradients in the dust-to-gas ratio. Even an $\epsilon$ of order $h_{\mathrm{g}}\ll 1$ can provide effective stabilisation against the VSI.

Appendix B Kelvin-Helmholtz instabilities

Axisymmetric simulations preclude Kelvin-Helmholtz instabilities associated with the vertical shear between a dusty midplane and the dust-free gas above and below it (Chiang & Youdin, 2010; Lee et al., 2010). Nevertheless, we can use the axisymmetric profiles reached in our simulations to examine stability against the KHI. In the literature this is usually assessed through the Richardson number, $\mathrm{Ri}=N_{z}^{2}/(\partial_{z}v_{\phi})^{2}$. Numerical simulations find KHI occurs for \mathrm{Ri}\;\raisebox{-3.44444pt}{\buildrel\textstyle<\over{\sim}}\;\mathrm{Ri}_{*}, where the critical value $\mathrm{Ri}_{*}$ depends on the dust layer parameters (Lee et al., 2010).

Fig. 17 show the evolution of $\mathrm{Ri}$ in selected simulations where dust settles down to a few per cent of $H_{\mathrm{g}}$. For $Z=0.01$ Lee et al. finds $\mathrm{Ri}_{*}\sim 0.1$, which is smaller than that the corresponding simulations. These cases are therefore stable against the KHI. However, extrapolating the numerical results of Lee et al. gives $\mathrm{Ri}_{*}\simeq 2$ for $Z=0.1$, whereas we find $\mathrm{Ri}=1$—2 in this case. The $Z=0.1$ simulation may thus be marginally unstable to KHI. Full 3D simulations will be required to check this.