Dust settling against hydrodynamic turbulence in protoplanetary discs
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 can sediment to of the gas scale-height, provided that —0.05, where 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.
keywords:accretion, accretion discs, hydrodynamics, methods: numerical, protoplanetary discs
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, , where and 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 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 (). 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 (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 : 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, 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 (LY15, Barker & Latter 2015). Such a vertical shear, , 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 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, , where 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 . We neglect disc self-gravity, magnetic fields, and viscosity. Cylindrical and spherical co-ordinates are centred on the star. In the discussions below, we use subscript ‘0’ to denote evaluation at , where is a reference radius. The Keplerian frequency is , where is the gravitational constant. We adopt a rotating frame with angular frequency .
The gas phase has density, pressure, and velocity fields . We model a single species of dust particles as a second, pressureless fluid with mass density and velocity . The fluid approximation for dust is valid for sufficiently small particles (Jacquet et al., 2011).
We also define the total mass density
the centre-of-mass velocity
and an effective energy
where is a constant defined below and is a time-independent gravitational potential (). Our numerical simulations evolve the hydrodynamic set , see §2.3.
2.1 Polytropic gas
We consider a locally polytropic gas with
where is a prescribed sound-speed profile and is the constant polytropic index. We take so the gas is nearly isothermal in its thermodynamic property. We set
where 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
and the gas disc aspect ratio is .
2.2 Tightly coupled dust
The dust fluid interacts with the gas via a drag force characterised by the particle stopping time . We consider particles in the Epstein regime with (Weidenschilling, 1977), where , 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 (Laibe & Price, 2014a). This is the characteristic decay timescale for the relative drift between dust and gas, , due to the mutual drag forces. The Stokes number can be defined as
We parameterise via a dimensionless stopping time such that
Thus 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 as the particle size.
In addition, we consider small particles with
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:
(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 (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
where the dust fraction can be obtained from the equation of state (Eq. 4)
and . We define the dust-to-gas ratio . 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
is the stellar potential in the thin-disc approximation (). 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 (). 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 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 allows us to define the (dimensionless) effective entropy of the system as
(LY17). A non-uniform entropy profile produces buoyancy forces. For this work the relevant quantity is the squared vertical buoyancy frequency
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 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
Here and below subscript ‘mid’ indicates the midplane value. The profile is chosen to minimise radial disc evolution, as described in Appendix B of CL18. The characteristic thickness in the dust-to-gas ratio, , is defined via
Here the dust layer thickness is set to so that initially the dust is vertically well-mixed with the gas (). For we fit the vertical dust density distribution to a Gaussian profile to obtain (see §4.2).
We scale the dust-to-gas ratio by specifying the metallicity at the reference radius
(Johansen et al., 2014). Recall are the surface densities in dust and gas, respectively.
We obtain the disc structure by assuming dynamical equilibrium and axisymmetry:
We set the midplane gas density to
with gas surface density
and we fix . The total density is . Finally, the disc orbital frequency is obtained from Eq. 21,
3.1 Vertical shear
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 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
(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
For the profiles we adopt (Eq. 18), we have . Then, considering the fiducial radius without loss of generality, we expect dust-induced buoyancy becomes important for midplane dust-to-gas ratios
for somewhat settled dust (), where the definition Eq. 19 was used. Thus we can expect dust-induced buoyant stabilisation of the VSI even when . LY17 give a similar estimate and confirm it with linear stability calculations (see their sections 6.1—6.2).
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 and a polar angle range such that , i.e. two gas scale-heights above and below the disc midplane at the reference radius. We fix . The standard resolution is ; with logarithmic and uniform spacing in and , respectively. This corresponds to a resolution of approximately 100 cells per at . We impose reflective boundary conditions in , and set radial boundary ghost zones to their initial values.
4.1 Physical parameters
Our fiducial parameter values are , , and . 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 . Smaller decreases the strength of the VSI turbulence and is expected to help settling.
We vary the stopping time . Larger , corresponding to larger particles, is expected to settle faster and avoid stirring by the VSI, at least initially.
We vary the overall metallicity . Increasing 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 and is slightly more dusty at (Chiang & Youdin, 2010) than our fiducial setup. Epstein drag stopping times in MMSN-like disc models vary as
where is the mass scale relative to the MMSN. At fixed , smaller particles reside at larger distances. In a disc and taking , we find corresponds to particle sizes from at au to at au; or between and for .
We analyse our simulation results using the diagnostics listed below. We denote mean value of a quantity by . Unless otherwise stated, the spatial average is taken over the shell . We also use or to denote the root-mean-squared value of . Time-averaged values are denoted with an additional subscript . Time is reported in units of the orbital period at and is denoted by .
The dust scale-height is obtained by fitting the dust density with
where is the mid-plane dust density. We then normalise by .
The midplane dust-to-gas ratio . We measure both average and maximum values within the sampling region.
The midplane vertical velocity .
Following Stoll et al. (2017), we measure the Reynolds stresses associated with vertical motions via
where denotes the deviation from the initial value. We then define the normalised vertical stress
as a dimensionless measure of the VSI strength, where is the initial density field. Note that in time evolution plots below, the average is also taken over the height of the disc.
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 . Here, we fix the particle size and metallicity such that and , respectively. Larger 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, ; the average and maximum midplane dust-to-gas ratio, and , respectively; the root-mean-squared midplane vertical velocity, ; and the dimensionless vertical stress, ; for , , , and . Fig. 2 compares several snapshots of the dust-to-gas ratio from these runs. It is evident that VSI turbulence hinders dust settling.
For the disc is strictly isothermal and does not develop the VSI. In this laminar disc, dust settling proceeds until the resolution limit () and . Introducing moderate VSI turbulence with prevents the dust from settling below . Consequently only increases to .
In the fully turbulent disc with 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 , with undulations near the surface. The banded structure in with characteristic separation 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 is consistent with Stoll et al. (2017), who found for .
In Fig. 4 we compare the time-averaged dust distribution and turbulence properties from a suite of runs with . The averaging period is taken between when the system enters a quasi-steady state and the end of the simulation. We omit the run here since in that case there is no physical steady state. We find increases rapidly with , but for the dust becomes well-mixed with and ( here). Notice the well-mixed and settled cases are separated by or . The former dependence is consistent with theoretical dust settling models, described below. As and increases with , this suggest that given a fixed particle size there exists a critical temperature gradient, , below which dust-settling is allowed. However, as shown below, vertical stresses and velocities also depends on particle size and metallicity.
(see also Youdin & Lithwick, 2007; Zhu et al., 2015), where the particle stopping time is measured at the disc midplane, and we regard the eddy timescale as a model parameter for simplicity. Setting gives a good match to the measured in the settled cases, while it somewhat underestimates the well-mixed cases. This value of 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 in our simulations. Here we fix and .
Fig. 5 compares the fiducial case with to runs with , , and times that particle size. Increasing allows faster settling, since the settling speed (Takeuchi & Lin, 2002). Larger particles can settle to a thinner layer before VSI-stirring, which begins at independently of . That is, the initial VSI growth rate is insensitive to particle size. However, VSI-stirring weakens with increasing , as ‘bounce’ to smaller values with larger particles after the VSI activates. This is also reflected in the much smaller with and compared to cases with and .
In the run, dust is almost perfectly remixed by the initial VSI. A quasi-steady state is apparently attained, but only until , when we observe spontaneous settling. This is shown in Fig. 6, suggesting that the prior steady state was unstable. The and 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 run attains a brief steady state () before settling; while at the end of the run 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 run. A dust clump develops at the outer boundary, possibly due to boundary conditions, which may be related to the increased dust mixing in . However, the sampling region around 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 . Here the time average is taken between the onset of the VSI and the end of the simulations. Significant settling occurs for , for which drops to and . The measured is consistent with the model described by Eq. 35 if we choose , again similar to eddy timescales found by Stoll & Kley (2016). Although less clear-cut, this figure also indicates that dust settling occurs when or .
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 . Larger corresponds to stronger dust feedback onto the gas. For our fiducial parameter values dust does not settle. We thus consider 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 of a few per cent.
Fig. 9 compares the disc evolution for , , , and . Increasing does not affect the settling timescale. Instead, larger delays the VSI, which allows the dust to settle to a thinner layer. This reflects the stabilising effect of particle back-reaction. For , the VSI remixes the dust at , and settles slowly afterwards. For and , the initial VSI only halts the settling process, and does not increase , unlike the and runs.
For , we find for . 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 and runs, does not reach unity within the simulation timescale, but is still significantly enhanced relative to the run. Comparing these runs to the previous simulations with increased (Fig. 5) show that increasing is more effective in raising the dust-to-gas ratio than increasing particle size. This is not surprising since, for a given a larger implies a larger (Eq. 20); let alone in the present cases where a smaller dust scale-height is reached with increasing compared to the same increase in .
The and cases show much reduced midplane vertical velocities and vertical turbulence compared to , which is consistent with the former cases reaching much smaller values of . The run also displays a slow drop in for , with a corresponding rise in ; while remains relatively constant. This signifies particle settling into midplane unhindered; in fact this settling is stabilising the system. The run attains early on, similar to the case, but it eventually drops to . This also reflects the slow settling observed after the initial VSI in that case.
Fig. 10 shows and at the end of the simulation. Maximum vertical velocities are reduced by at a factor of two from the standard case with (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 near the boundaries . This reflects VSI activity in the dust-poor layers above and below the midplane. The spike at 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 .) Although vertical velocities are reduced towards the midplane, notice that . 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 simulation. We show both the actual vertical shear rate, ; and that associated with the imposed temperature gradient , , 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, . Away from the midplane at , buoyancy frequencies vanish, whereas the temperature-related, destabilising vertical shear rate is larger there. About the midplane (), 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 run here since in that case. For this ratio is , though rising later on, which is consistent with the slow settling. For and this ratio is , 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 . We also plot several cases with .
For we find stabilisation by dust-loading is most dramatic for —; while for dust settles to with little variation with respect to . On the other hand, for dust can settle to with ; suggesting the critical for settling scales as . We find the settled dust scale-heights can be fitted with Eq. 35, but this required , which is smaller than that found by Stoll & Kley (2016). This indicates that increasing the metallicity reduces the turbulence eddy timescales. We find 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 , i.e. the measured vertical turbulence is associated with the gas layers above and below the midplane, hence flattens off. Nevertheless, we find settled cases still obey .
These results indicate that in weakly VSI-turbulent discs, dust settling is sensitive to the local metallicity. This implies that small changes in , 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.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 , 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 and we find is needed for dust to settle down to ; while in less turbulent discs with this only required , 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 becomes comparable to the disc’s aspect-ratio , i.e. a few per cent in PPDs.
Our simulations show that generally depends on , , and . All three parameters affect the VSI turbulence strength (the dimensionless vertical Reynolds stress) and the midplane vertical velocity . We find that the settled and well-mixed cases are separated by or . 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 () and abundance can settle. Recent analyses indicate for VSI turbulence (Barker & Latter, 2015; Latter & Papaloizou, 2018), which is consistent with our simulations. We thus expect particles with to settle, even in the limit . However, if , then 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 via particle growth, raising the local metallicity via radial dust drift, or reducing the radial temperature gradient . 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 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 drops significantly with decreasing , it is still non-zero about the dusty midplane. Such a layered structure in 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.
6.3 Implications for planetesimal formation
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 . If we take and solar metallicity , then by small we mean particles with of . At a distance of au, this roughly corresponds to particle sizes of , assuming an internal density (see Eq. 31).
In the above example, particles smaller than will remain lofted. However, if they can grow to 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 in two regards. First, a larger allows to be attained at larger (see Eq. 20). That is, the required degree of settling is smaller. Second, increasing actively stabilises the disc against the VSI, thereby allow for effective settling.
Coincidentally, non-linear numerical simulations of the SI indicate is required for strong clumping of particles with — (Johansen et al., 2009b; Bai & Stone, 2010a) ; while smaller particles with require (Yang et al., 2017). These particle clumps can then undergo GI. One also requires 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 at au, and assuming . They infer a dimensionless vertical diffusion parameter , assuming passive particles (Fromang & Nelson, 2009). Interestingly, these values are consistent with our simulation222Note that we used while Pinte et al. and Jin et al. assumed and , respectively. However, the destabilising effect of increasing offsets the stabilising effect of decreasing (Nelson et al., 2013). shown in Fig. 8, if we identify as our dimensionless vertical Reynolds stress. Using Eq. 31, corresponds to particle sizes m at au, similar to that adopted in the HL Tau disc model of Jin et al. (2016). For mm-sized particles, requires a disc times more massive than the MMSN; or grains with internal density .
For dusty VSI turbulence we find . This may be used to estimate the local metallicity. For example, let us assume and particle sizes of , and a massive (young) disc with . This gives at au. Then Fig. 14 suggest Pinte et al.’s estimate, , can be achieved with . On the other hand, if 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 to such that . Then the increased (decreased) metallicity at () would stabilise (destabilise) the disc there. Now, if the gas density contrast between and is smaller than that in , i.e. , so that the contrast in stopping times is outweighed by that in metallicity, then dust settling (mixing) should be promoted at (). Ultimately, , assuming similar gas scale heights if . 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 can settle to ; whereas we find such particles remain perfectly mixed with the gas. On the other hand, Stoll & Kley found particles with only settles to , whereas we find particles can already settle to , 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.
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, (essentially the vertical Reynolds stress-to-pressure ratio); particle size or Stokes number, ; and solid abundance (metallicity) on dust settling.
We model VSI turbulence explicitly by imposing a fixed radial temperature gradient . However, we find the resulting stress can also depend on and . Nevertheless, our results indicate that dust settles when the resulting . In our simulations this is equivalent to having a midplane vertical velocity dispersion . 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 and/or , where 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 , corresponding to grain sizes of m at au in MMSN-like discs, we find that enables dust reach 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.
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.
- ALMA Partnership et al. (2015) ALMA Partnership et al., 2015, ApJ, 808, L3
- Alexander et al. (2014) Alexander R., Pascucci I., Andrews S., Armitage P., Cieza L., 2014, Protostars and Planets VI, pp 475–496
- Bai (2016) Bai X.-N., 2016, ApJ, 821, 80
- Bai & Stone (2010a) Bai X.-N., Stone J. M., 2010a, ApJ, 722, 1437
- Bai & Stone (2010b) Bai X.-N., Stone J. M., 2010b, ApJ, 722, L220
- Balbus & Hawley (1991) Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214
- Balsara et al. (2009) Balsara D. S., Tilley D. A., Rettig T., Brittain S. D., 2009, MNRAS, 397, 24
- Barge & Sommeria (1995) Barge P., Sommeria J., 1995, A&A, 295, L1
- Barker & Latter (2015) Barker A. J., Latter H. N., 2015, MNRAS, 450, 21
- Barranco (2009) Barranco J. A., 2009, ApJ, 691, 907
- Blum (2018) Blum J., 2018, Space Sci. Rev., 214, 52
- Chen & Lin (2018) Chen J.-W., Lin M.-K., 2018, MNRAS, 478, 2737
- Chiang (2008) Chiang E., 2008, ApJ, 675, 1549
- Chiang & Youdin (2010) Chiang E., Youdin A. N., 2010, Annual Review of Earth and Planetary Sciences, 38, 493
- Dittrich et al. (2013) Dittrich K., Klahr H., Johansen A., 2013, ApJ, 763, 117
- Dubrulle et al. (1995) Dubrulle B., Morfill G., Sterzik M., 1995, Icarus, 114, 237
- Flock et al. (2017) Flock M., Nelson R. P., Turner N. J., Bertrang G. H.-M., Carrasco-González C., Henning T., Lyra W., Teague R., 2017, ApJ, 850, 131
- Fromang & Lesur (2017) Fromang S., Lesur G., 2017, preprint, (arXiv:1705.03319)
- Fromang & Nelson (2009) Fromang S., Nelson R. P., 2009, A&A, 496, 597
- Fromang & Papaloizou (2006) Fromang S., Papaloizou J., 2006, A&A, 452, 751
- Gammie (1996) Gammie C. F., 1996, ApJ, 457, 355
- Gammie (2001) Gammie C. F., 2001, ApJ, 553, 174
- Gibbons et al. (2012) Gibbons P. G., Rice W. K. M., Mamatsashvili G. R., 2012, MNRAS, 426, 1444
- Goldreich & Ward (1973) Goldreich P., Ward W. R., 1973, ApJ, 183, 1051
- Gonzalez et al. (2017) Gonzalez J.-F., Laibe G., Maddison S. T., 2017, MNRAS, 467, 1984
- Helled et al. (2014) Helled R., et al., 2014, Protostars and Planets VI, pp 643–665
- Jacquet et al. (2011) Jacquet E., Balbus S., Latter H., 2011, MNRAS, 415, 3591
- Jin et al. (2016) Jin S., Li S., Isella A., Li H., Ji J., 2016, ApJ, 818, 76
- Johansen & Youdin (2007) Johansen A., Youdin A., 2007, ApJ, 662, 627
- Johansen et al. (2006) Johansen A., Henning T., Klahr H., 2006, ApJ, 643, 1219
- Johansen et al. (2009a) Johansen A., Youdin A., Klahr H., 2009a, ApJ, 697, 1269
- Johansen et al. (2009b) Johansen A., Youdin A., Mac Low M.-M., 2009b, ApJ, 704, L75
- Johansen et al. (2014) Johansen A., Blum J., Tanaka H., Ormel C., Bizzarro M., Rickman H., 2014, Protostars and Planets VI, pp 547–570
- Kanagawa et al. (2017) Kanagawa K. D., Ueda T., Muto T., Okuzumi S., 2017, ApJ, 844, 142
- Klahr & Hubbard (2014) Klahr H., Hubbard A., 2014, ApJ, 788, 21
- Klahr et al. (2018) Klahr H., Pfeil T., Schreiber A., 2018, Instabilities and Flow Structures in Protoplanetary Disks: Setting the Stage for Planetesimal Formation. p. 138, doi:10.1007/978-3-319-55333-7_138
- Laibe & Price (2014a) Laibe G., Price D. J., 2014a, MNRAS, 440, 2136
- Laibe & Price (2014b) Laibe G., Price D. J., 2014b, MNRAS, 444, 1940
- Latter (2016) Latter H. N., 2016, MNRAS, 455, 2608
- Latter & Papaloizou (2018) Latter H. N., Papaloizou J., 2018, MNRAS, 474, 3110
- Lee et al. (2010) Lee A. T., Chiang E., Asay-Davis X., Barranco J., 2010, ApJ, 718, 1367
- Lesur & Latter (2016) Lesur G. R. J., Latter H., 2016, MNRAS, 462, 4549
- Lesur & Papaloizou (2009) Lesur G., Papaloizou J. C. B., 2009, A&A, 498, 1
- Lesur & Papaloizou (2010) Lesur G., Papaloizou J. C. B., 2010, A&A, 513, A60
- Lin & Youdin (2015) Lin M.-K., Youdin A. N., 2015, ApJ, 811, 17
- Lin & Youdin (2017) Lin M.-K., Youdin A. N., 2017, ApJ, 849, 129
- Lissauer (1993) Lissauer J. J., 1993, ARA&A, 31, 129
- Lyra (2014) Lyra W., 2014, ApJ, 789, 77
- Lyra & Lin (2013) Lyra W., Lin M.-K., 2013, ApJ, 775, 17
- Lyra & Umurhan (2018) Lyra W., Umurhan O., 2018, preprint, (arXiv:1808.08681)
- Manger & Klahr (2018) Manger N., Klahr H., 2018, MNRAS, 480, 2125
- Marcus et al. (2015) Marcus P. S., Pei S., Jiang C.-H., Barranco J. A., Hassanzadeh P., Lecoanet D., 2015, ApJ, 808, 87
- Mignone et al. (2007) Mignone A., Bodo G., Massaglia S., Matsakos T., Tesileanu O., Zanni C., Ferrari A., 2007, ApJS, 170, 228
- Mignone et al. (2012) Mignone A., Flock M., Stute M., Kolb S. M., Muscianisi G., 2012, A&A, 545, A152
- Nakagawa et al. (1986) Nakagawa Y., Sekiya M., Hayashi C., 1986, Icarus, 67, 375
- Nelson et al. (2013) Nelson R. P., Gressel O., Umurhan O. M., 2013, MNRAS, 435, 2610
- Pinilla & Youdin (2017) Pinilla P., Youdin A., 2017, in Pessah M., Gressel O., eds, Astrophysics and Space Science Library Vol. 445, Astrophysics and Space Science Library. p. 91, doi:10.1007/978-3-319-60609-5_4
- Pinte et al. (2016) Pinte C., Dent W. R. F., Ménard F., Hales A., Hill T., Cortes P., de Gregorio-Monsalvo I., 2016, ApJ, 816, 25
- Price & Laibe (2015) Price D. J., Laibe G., 2015, MNRAS, 451, 813
- Raettig et al. (2015) Raettig N., Klahr H., Lyra W., 2015, ApJ, 804, 35
- Raymond et al. (2014) Raymond S. N., Kokubo E., Morbidelli A., Morishima R., Walsh K. J., 2014, Protostars and Planets VI, pp 595–618
- Rice et al. (2004) Rice W. K. M., Lodato G., Pringle J. E., Armitage P. J., Bonnell I. A., 2004, MNRAS, 355, 543
- Richard et al. (2016) Richard S., Nelson R. P., Umurhan O. M., 2016, MNRAS, 456, 3571
- Riols & Lesur (2018) Riols A., Lesur G., 2018, A&A, 617, A117
- Safronov (1969) Safronov V. S., 1969, Evoliutsiia doplanetnogo oblaka.
- Shi & Chiang (2013) Shi J.-M., Chiang E., 2013, ApJ, 764, 20
- Shi et al. (2016) Shi J.-M., Zhu Z., Stone J. M., Chiang E., 2016, MNRAS, 459, 982
- Simon et al. (2016) Simon J. B., Armitage P. J., Li R., Youdin A. N., 2016, ApJ, 822, 55
- Squire & Hopkins (2018) Squire J., Hopkins P. F., 2018, MNRAS, 477, 5011
- Stoll & Kley (2014) Stoll M. H. R., Kley W., 2014, A&A, 572, A77
- Stoll & Kley (2016) Stoll M. H. R., Kley W., 2016, A&A, 594, A57
- Stoll et al. (2017) Stoll M. H. R., Kley W., Picogna G., 2017, A&A, 599, L6
- Takahashi & Inutsuka (2014) Takahashi S. Z., Inutsuka S.-i., 2014, ApJ, 794, 55
- Takeuchi & Lin (2002) Takeuchi T., Lin D. N. C., 2002, ApJ, 581, 1344
- Umurhan et al. (2016) Umurhan O. M., Shariff K., Cuzzi J. N., 2016, ApJ, 830, 95
- Weidenschilling (1977) Weidenschilling S. J., 1977, MNRAS, 180, 57
- Whipple (1972) Whipple F. L., 1972, in Elvius A., ed., From Plasma to Planet. p. 211
- Yang & Johansen (2016) Yang C.-C., Johansen A., 2016, ApJS, 224, 39
- Yang et al. (2017) Yang C.-C., Johansen A., Carrera D., 2017, A&A, 606, A80
- Yang et al. (2018) Yang C.-C., Mac Low M.-M., Johansen A., 2018, ApJ, 868, 27
- Youdin & Goodman (2005) Youdin A. N., Goodman J., 2005, ApJ, 620, 459
- Youdin & Johansen (2007) Youdin A., Johansen A., 2007, ApJ, 662, 613
- Youdin & Lithwick (2007) Youdin A. N., Lithwick Y., 2007, Icarus, 192, 588
- Zhu et al. (2015) Zhu Z., Stone J. M., Bai X.-N., 2015, ApJ, 801, 81
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 (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 and particles with . Fig. 16 compares the simulation with , as presented in the main text, to a disc with . They evolve similarly until , at which point dust in the disc undergoes complete remixing by the VSI; whereas dust continues to settle in the disc. In the latter case at , 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 of order 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, . Numerical simulations find KHI occurs for , where the critical value depends on the dust layer parameters (Lee et al., 2010).
Fig. 17 show the evolution of in selected simulations where dust settles down to a few per cent of . For Lee et al. finds , 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 for , whereas we find —2 in this case. The simulation may thus be marginally unstable to KHI. Full 3D simulations will be required to check this.