© ESO, 2016

1. Introduction

The study of gas drag on mm to dm-sized particles (pebbles) is essential to understanding the formation of planets. Vertical sedimentation owing to drag on small particles in the protoplanetary disc is necessary for the creation of a dense midplane of solids from which larger objects can grow (Youdin & Lithwick 2007). Conversely, the same drag force is also responsible for the rapid radial migration of pebbles in the midplane (on 100 yr timescales in the terrestrial region, Weidenschilling 1977). This is the main barrier for continued growth to larger than cm-to-m sizes by collisions (Brauer et al. 2008; Birnstiel et al. 2012), unless particles can have extremely low internal densities (Kataoka et al. 2013; Krijt et al. 2015).

The radial drift hurdle can be avoided through two mechanisms that also critically rely on gas drag. Firstly, pebbles can be concentrated hydrodynamically, so that the resulting clouds collapse gravitationally to planetesimals of ~100 km in size (for recent reviews on different planetesimal formation models, see Johansen et al. 2014; Chiang & Youdin 2010). Secondly, large planetesimals can accrete the remaining drifting pebbles and grow to planetary sizes (Lambrechts & Johansen 2012, 2014; Guillot et al. 2014).

Not only the drag on the particles is important, but also the backreaction of the particles on the gas. Initially, it was proposed that a secular instability on a settled dust layer could lead to local particle pileups (Goodman & Pindor 2000). The pileup would originate from a process resembling plate drag, where the drag force is assumed to collectively act on a monolithic particle midplane. This assumption is nevertheless questionable (Youdin & Chiang 2004) and numerical studies (Weidenschilling 2006) have not recovered the instability proposed by Goodman & Pindor (2000). Nevertheless, this work paved the way for a further investigation into the role of the backreaction force from gas drag. Youdin & Goodman (2005) identified a linear instability in the disc midplane. Their breakthrough result demonstrated that infinitesimal perturbations grow on an orbital time-scale when the dust-to-gas ratio is around unity or higher. This instability leads to spontaneous particle clumping, triggering the gravitational collapse that results in the formation of planetesimals. In a series of papers (Youdin & Johansen 2007; Johansen & Youdin 2007; Johansen et al. 2009, 2012) the linear and nonlinear evolution of this instability were numerically investigated in detail. These results were independently confirmed and further explored by several other groups (Bai & Stone 2010a,b,c; Miniati 2010; Kowalik et al. 2013).

Several criteria for the streaming instability to achieve particle clumping have been identified:

• a disc with slightly supersolar dust-to-gas ratio (Johansenet al. 2009; Bai &Stone 2010c);

• particles of Stokes number τ_f ~ 0.05–0.5, approximately between mm and dm in size (Johansen & Youdin 2007; Bai & Stone 2010a,c; Carrera et al. 2015); and

• low radial pressure support in the disc (Bai & Stone 2010c).

Further investigations are moving towards a more global understanding of the effects of the streaming instability, by expanding the simulation domain in the azimuthal (Kowalik et al. 2013) or vertical direction (Yang & Johansen 2014). Additionally the streaming instability is placed in a larger context by incorporating magnetized turbulence (Johansen et al. 2007), dust coagulation models (Dra¸żkowska & Dullemond 2014), or vortex formation (Raettig et al. 2015).

In this paper we take a step back and study the general process of particle sedimentation in flows with a dust loading comparable to the gas density. The aim is twofold. Firstly, we hope to gain theoretical insight into particle sedimentation and more complex drag instabilities, such as the streaming instability and the photoelectric instability (Lyra & Kuchner 2013). Secondly, the sedimentation of particles is accessible to laboratory experiments, thus allowing for a potential experimental confirmation of a particle drag instability.

Of specific interest is the question whether any particle clumping even occurs in a mass-loaded particle rain. From previous analytic work on the streaming instability, we do not expect a linear instability to be present, because of the lack of rotation in a pure sedimentation problem. This removes the Coriolis force which is deemed necessary for the streaming instability to operate (Jacquet et al. 2011; Youdin & Goodman 2005). Nevertheless, in our physical set up (described in Sect. 2) a nonlinear drafting instability is clearly present. The results are described in Sect. 3. The implications of this instability are discussed for particle sedimentation in protoplanetary discs, chondrule formation, and the envelopes of giant planets (Sect. 4). We also place our results in the context of planned laboratory experiments (Sect. 5.1). We summarize our findings in Sect. 6.

2. Mass-loaded particle rain

2.1. Model equations

We study the differential motion between particles initially moving with terminal velocity and stationary gas in hydrostatic balance. The dynamics of the gas component are described by $\begin{array}{ccc}& & \\ & & \end{array}$where ρ is the gas density, u the gas velocity, g the gravitational acceleration, P the pressure, ϵ = ρ_p/ρ is the local dust-to-gas ratio, v the particle fluid velocity, and ν the viscosity. The drag term from the particles onto the gas depends on the friction time of the particle (in the Epstein drag regime, Epstein 1924), $\begin{array}{ccc}{\mathit{t}}_{\mathrm{f}}\mathrm{=}\frac{{\mathit{\rho }}_{\mathrm{•}}\mathit{R}}{\mathit{\rho }{\mathit{v}}_{\mathrm{th}}}\mathit{,}& & \end{array}$(3)where R,ρ_• are the radius and solid density of the particle. The thermal velocity v_th is approximately equal to the local, isothermal gas sound speed ${\mathit{v}}_{\mathrm{th}}\mathrm{=}\sqrt{\mathrm{8}\mathit{/}\mathit{\pi }}{\mathit{c}}_{\mathrm{s}}$.

We investigate a regime where we assume an efficient coupling between particles and gas. In this case, the viscous diffusion time for momentum transport between an average particle pair located a distance l_pair apart is shorter than the friction time of a single particle, $\begin{array}{ccc}{\mathit{t}}_{\mathit{\nu ,}\mathrm{pair}}\mathrm{=}\frac{{\mathit{l}}_{\mathrm{pair}}^{\mathrm{2}}}{\mathit{\nu }}\mathrm{=}\frac{{\mathit{n}}_{\mathrm{p}}^{\mathrm{-}\mathrm{2}\mathit{/}\mathrm{3}}}{\mathit{\nu }}\mathit{<}{\mathit{t}}_{\mathrm{f}}\mathit{,}& & \end{array}$(4)with n_p the particle number density. The particles can then be described by a pressureless fluid (a formal derivation can be found in Jacquet et al. 2011) with continuity and momentum equation $\begin{array}{ccc}& & \\ & & \end{array}$In the remainder of the paper we employ friction units: the friction time t_f as time unit and the friction length ${\mathit{l}}_{\mathrm{f}}\mathrm{=}\mathit{g}{\mathit{t}}_{\mathrm{f}}^{\mathrm{2}}$ as length unit. Velocities can then be expressed in units of terminal velocity v_f = gt_f. The criterion expressed in Eq. (4), for example, reduces to ${\mathit{n}}_{\mathrm{p}}^{\mathrm{\prime }\mathrm{-}\mathrm{2}\mathit{/}\mathrm{3}}\mathit{/}{\mathit{\nu }}^{\mathrm{\prime }}\mathit{<}\mathrm{1}$. We preserve the prime notation in the following sections to explicitly denote quantities expressed in friction units.

The use of the friction time as a unit of time is possible for the sedimentation problem, because there are no rotation terms that would necessarily introduce the additional timescale of the orbital Keplerian frequency, as is for example the case for the streaming instability (Youdin & Goodman 2005).

Expressed in friction units, the model equations leave us with only three free dimensionless parameters:

• the viscosity, ${\mathit{\nu }}^{\mathrm{\prime }}\mathrm{=}\mathit{\nu }\mathit{/}\mathrm{\left(}{\mathit{g}}^{\mathrm{2}}{\mathit{t}}_{\mathrm{f}}^{\mathrm{3}}\mathrm{\right)}$, which is the inverse of the Reynolds number (Re) in terms of the terminal velocity and the friction length;

• the sound speed, ${\mathit{c}}_{\mathrm{s}}^{\mathrm{\prime }}\mathrm{=}{\mathit{c}}_{\mathrm{s}}\mathit{/}\mathrm{\left(}\mathit{g}{\mathit{t}}_{\mathrm{f}}\mathrm{\right)}$, which is the inverse of the Mach number (Ma) in terms of the terminal velocity; and

• the dust-to-gas ratio, ϵ.

The latter is arguably the most important, because we aim to work in the incompressible limit (Ma ≪ 1) and we face a lower bound on the viscosity imposed by the requirement for numerical stability.

A major benefit from this choice of units is that our calculations do not require us to specify a particle size (Eq. (3)). Thus our results can be freely scaled to the desired particle size in the context of protoplanetary discs (Sect. 4) or a laboratory setting (Sect. 5.1).

Fig. 1 Development of particle swarms by a drafting instability. Displayed is the evolution of the particle density in the two-dimensional simulation run2 (see Table 1). The left-most panel illustrates the initial conditions: a stratified gas column in the vertical directions with particles sedimenting at terminal velocity, placed randomly throughout the simulation domain (1 lf wide and 20 lf high, note that the figure aspect ratio is enlarged in the x-direction). In the other panels (time t = 10,15 tf), regions marginally overdense with particles locally break the stratification equilibrium and accelerate downwards while in particle-poor regions, a deceleration from terminal velocity occurs. This leads to continued particle pileups through drafting, originating from gas being dragged by the particles. The particle and fluid components unmix and remain in this state. At the end of the simulation, when most particles have sedimented out of the simulation domain, the maximal particle density has increased by a factor 10.

2.2. Numerical implementation

In our numerical simulations, performed with the Pencil Code¹, we do not employ the particle fluid description used for the analytical calculations for our main results, but instead use a Lagrangian super-particle approach. Particles are implemented as super-particles, which represent swarms of physical particles. This is important, especially in the nonlinear regime where it is desirable to allow particle trajectory crossing and steep density gradients (Youdin & Johansen 2007). We have nevertheless used the particle fluid approach numerically as well to verify some of our work in the linear regime. The assignment of drag forces on the particles and on the gas is further described in Youdin & Johansen (2007) and in Johansen & Youdin (2007). Drag is calculated on particles assuming a constant friction time. We have also made use of the block domain decomposition for particle-load balancing amongst processors (Johansen et al. 2011).

For the initial condition, we set up a gas column in hydrostatic equilibrium, taking the drag from the particles on the gas into account. Particles are typically distributed randomly, with the possibility of adding perturbations (discussed in more detail in Sect. 3.1) and initiated with their vertical velocity equal to the terminal velocity. The gas stratification then takes the form $\begin{array}{ccc}\mathit{\rho }\mathrm{=}\mathrm{(}{\mathit{\rho }}_{\mathrm{p}}\mathrm{+}{\mathit{\rho }}_{\mathrm{b}}\mathrm{)}\exp \left(\mathrm{-}\frac{\mathit{g}}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}\mathit{z}\right)\mathrm{-}{\mathit{\rho }}_{\mathrm{p}}\mathit{,}& & \end{array}$(7)see Appendix A for more details. This initial condition works well, but we nevertheless find that in heavily elongated simulation domains, the stratification of the gas combined with a uniform particle distribution triggers a vertical gas density wave, that dissipates over time. There is also a small upwards advection of gas as the top of the domain gets cleared of particles. These effects restrict our simulation domain in practice to approximately 20  l_f in the vertical direction (at c_s/v_f = 10).

A full list of the performed simulations can be found in Tables 1 and 2. Below, we describe the nominal numerical set up in detail.

The boundary condition is set to be periodic in the horizontal direction, for both the particle and gas component. In the vertical direction, particles are removed from the simulation when crossing the edge of the simulated domain. The vertical boundary condition on the gas is symmetric (vanishing first derivative) in all quantities, except for the vertical gas velocity which is antisymmetric (vanishing value). This boundary condition effectively puts a solid surface at the bottom of the simulation domain on which the gas is supported.

We use an ideal gas equation of state with adiabatic index γ = 5/3. The sound speed is set to be c_s = 10v_f (unless mentioned otherwise) to approach the incompressible limit. The Pencil Code is a code optimized for both subsonic and mildly transonic flows, but we found a Mach number of 0.1 sufficient to probe the incompressible regime of interest.

We found a choice of 16 superparticles per grid cell is sufficiently high to model a coherent fluid and reduce particle noise (see also Appendix C). We have consistently used a physical viscosity treatment but, for runs with extended domain and high mass loading (run1.01, run1.02), we added sixth order “artificial viscosity” (Haugen & Brandenburg 2004). We employed the minimal amount necessary to prevent numerical artefacts from developing. The grid Reynolds number is minimally 32 times smaller than the Reynolds number in friction units. Most simulations were performed in two dimensions, but we have verified our results in three dimensions as well (run3d.4, Fig. 2), showing little difference. Nevertheless, these numerically expensive 3D runs are of interest for further study.

3. Numerical results showing spontaneous particle concentrations

3.1. Nonlinear behaviour

Fig. 2 Top panel: evolution of the particle density (normalized to the background dust-to-gas ratio ϵ0), for three different resolutions in 2D simulations and a 3D simulation (run3d.4, gray). Bottom panel: evolution of the instability based on the horizontal gas dispersion, $\sqrt{\mathrm{⟨}{\mathit{u}}_{\mathit{x}}^{\mathrm{2}}\mathrm{⟩}}$, (as an alternative tracer to the maximal particle density, which is an intrinsically noisy variable). We used a similar color coding as for the top panel. Note for the runs with ϵ0 = 4, we displaced the dashed curves for clarity with a factor 10 downwards.

From previous analytic studies of the streaming instability (Youdin & Goodman 2005; Jacquet et al. 2011), it is well known that rotation is an essential ingredient for the linear phase of particle clumping. The interpretation proposed by Jacquet et al. (2011) is that the Coriolis force is necessary to create a pressure maximum supported by geostrophic balance. We have repeated the analysis in Appendix A for completeness. This demonstrates that the sedimentation model without rotation, as discussed in this paper, is not expected to show a linear instability, under the assumption of incompressible gas.

Our numerical results nevertheless show that even a minimal disturbance of the sedimenting particle component with ϵ = 1 leads to spontaneous clumping of material, resulting in the particles unmixing and sedimenting out of the fluid. Figure 1 illustrates the process. Particles located in initially weakly overdense regions sediment faster, dragging the gas along, resulting in a drafting effect which pulls in more particles. For clarity, the drafting mentioned here is the result of the collective motion of a swarm of particles and the resulting gas drag-back reaction, and not from individual particle slipstreaming. At the same time, the opposite occurs in regions less dense than the mean particle density. These effects amplify each other, which results in dense swarms of particles forming, which can undergo secondary instabilities. For example, one can see in the bottom of the last panel of Fig. 1 a particle cloud resembling a characteristic Rayleigh-Taylor mushroom, which we will discuss in more detail below.

Fig. 3 Correlation between the local dust-to-gas ratio ϵ and the deviations in the velocity of the gas uz (full line) and particles vz (dotted line). To aid interpretation, we subtracted the mean velocity of the sedimenting particles ⟨vz⟩ from the particle velocity and a small artificial upwards mean gas velocity ⟨uz⟩ from the gas velocity. Displayed are different times t = 0,5,10 tf, in respectively black, blue, and red, based on the region between z = 0–10 lf in run2. The standard deviations on the binned averages are relatively large, σ ≈ 0.03–0.04 vf for the gas and particle velocities between t = 5–10 tf. The black dashed line corresponds to α = 0.02 in the toy model (Sect. 3.5).

3.2. Particle noise perturbation

We begin by considering the evolution of the sedimenting particles, when they are placed randomly in the simulation domain. This corresponds to a noisy initial condition for the particle distribution, with maximal relative changes in the particle density on the order of 50% for our nominal 2D resolution. For more detail, see Appendix C.

Figure 2 shows the time evolution of the maximal particle density in the simulated domain (for different background dust-to-gas ratios). Particle overdensities reach 10 times the average value, although they are still growing slowly towards the end of the simulations when the particles fall out of the box. Over time, particles sediment out of the simulated domain so, at late times, fewer and fewer particles are traced. For clarity, we also show the evolution of the horizontal gas velocity dispersion, which is a less noisy measurement than the particle density. This illustrates the exponential nature of the instability as well.

Our numerical results stand in contrast to the stable state predicted by linear stability analysis. It appears that the main driver for the particle clumping is an imbalanced stratification. We recall that the particle loading of the sedimenting particles is taken into account when setting up the equilibrium state (Eq. (7)). This balance can be broken along the x-direction by regions with a higher, or less high, particle density compared to the mean value. Apparently, the fluid and particles have no means of finding a global equilibrium state in the x-direction in the response to the particle fluctuation. Instead, the particle component breaks into dense swarms.

This interpretation is supported by the correlation between overdense regions and their increased settling speeds shown in Fig. 3. We have binned the surrounding gas and particle velocity in the grid cell for every particle in the simulation of run2, revealing that, on average, gas and particles sediment about 1% faster or slower for order unity fluctuations in the dust-to-gas ratio.

From an inspection of our numerical results at different resolutions and spatial scales, we find that the instability tends to originate on the smallest available scale, near the grid scale in simulations with minimal viscosity (the dependency on the viscosity is discussed in Sect. 3.6 in more detail). This makes it computationally challenging to characterize the instability, as increased resolutions do not necessarily better resolve the characteristic scale of interest. Instead, the instability takes place a little faster and remains quantitatively similar (Fig. 2).

The view that the drafting instability arises from an imbalanced pressure stratification suggests that the instability is nonlinear, since the particle fluctuations that drive the horizontal imbalance must be seeded from the initial condition. This implies that the growth rate of the instability should decrease with an increasing particle number, as randomly placed particles have a decreased effective density fluctuation with increased particle number. In Appendix C, we show that indeed the fastest growth is found when reducing the particle number to just four per cell. However, we are not able to completely shut off the instability at larger particle number (64 particles per grid cell), indicating that the instability may operate even in the limit of very high particle number.

In the next sections, we investigate how the sedimenting particles react to different perturbations of the system to gain further insight in the nonlinear phase. We then propose that a toy model that can capture the dependency of the instability on the metallicity and Reynolds number (Sects. 3.5–3.6).

Fig. 4 Left panel: sedimenting particles with a kz = π/ 2-mode in the particle density. Represented is a snapshot of runRT at t = 5tf. No instabilities develop, suggesting that the drafting instability is not directly related to the Rayleigh-Taylor instabilities that should occur in this set up. Right panel: evolution of the sedimenting particle wave at an arbitrary placed slice at x = −0.23. Different curves represent different times (t = 0,1,2,...,6tf), that for clarity are offset by 0.2. The particle density ρp(z) is given in black, while the red dashed line gives the gas velocity uz(z). Particle ridges do not approach each other, as indicated by the blue line that tracks the position of a point advected with velocity vf. The pressure that is supporting the mass-loaded stratification adapts to compensate for small changes in the particle density. This allows uz(z) to remain zero while the particles sediment.

3.3. Wave perturbation

We now verify numerically that the drafting instability is not related to either the Rayleigh-Taylor instability or the Kelvin-Helmholtz instability by perturbing the system with either a purely vertical or purely horizontal mode in the particle distribution.

We first present results of vertical wave perturbation of the particle density, which can be seen in Fig. 4. Within the time that particles sediment out of the domain, no instabilities can be detected. The right panel of Fig. 4 shows that the particles simply sediment at terminal velocity. At the same time, the gas component does not react to the perturbation. We have verified that this result even stands when feeding additional particle noise to the simulation. We do not see any Rayleigh-Taylor-like instabilities (that occur in a hydrostatic fluid with a dense layer on a lighter one, Drazin & Reid 2004). This experiment also demonstrates that the mechanism that concentrates particles is at least two-dimensional. In 1D simulations, with only a vertical perturbation, the ridges of increased particle density do not approach each other.

Fig. 5 Sedimenting particles with a kx = π/ 2-mode in the particle density. We show different details (left) from a larger simulated domain (right), at different times (t = 1,6,9,12tf), based on runKH. Initially, we added noise in the particle velocity to perturb the boundary between particle-rich and particle-poor regions. We believe that we see a “parasitic” Kelvin-Helmholtz instability appear over time.

We also experiment with a horizontal wave in the particle distribution. To perturb the interface between the particle-rich and particle-poor region, we displace the particles initially by giving them a random velocity kick of v = 0.1v_f. The resulting evolution is shown in Fig. 5. The particle dense columns supersediment at an accelerated rate. This differential velocity between particle-poor and particle-rich columns drive an instability reminiscent of the Kelvin-Helmholtz instability (Drazin & Reid 2004). However, in this case, the denser fluid that is used in the classical description of the instability is replaced with a fluid that contains an overdensity in particles. These particle-loaded Kelvin-Helmholtz instabilities have been studied before in the context of molecular clouds (Hendrix & Keppens 2014). A characteristic feature of the Kelvin-Helmholtz instability is the emergence of v-shaped wings. These can be identified in the bottom panel of Fig. 5. Similar features are also seen in simulations of the early linear evolution of the streaming instability in unstratified discs (see, for example, Fig. 2 in Johansen & Youdin 2007). Therefore these wings might be a general feature of particle-gas instabilities.

We expect that this parasitic Kelvin-Helmholtz instability operates in the fully mixed state of our noise simulations when dense regions start sedimenting out. It may thus play an important role in the late nonlinear evolution. However, since we only see this type of behaviour resembling Kelvin-Helmholtz instabilities in the case of a large perturbation, we do not believe it to be the origin of the instability in the initial noise runs.

Fig. 6 The evolution of particle droplets, at different times (t = 0,5,10tf), resulting in the emergence of characteristic Rayleigh-Taylor mushroom clouds. The color bar gives the color scale for the shown particle density. Results from runEGG.

3.4. Eggbox perturbation

Finally, we also perturb the system with an eggbox-like perturbation, to investigate the formation and evolution of particle swarms, which we will term particle droplets. The initial particle density perturbation is of the form $\begin{array}{ccc}{\mathit{\rho }}_{\mathrm{p}}\mathrm{\left(}\mathit{x,z}\mathrm{\right)}\mathrm{=}\mathit{A}\sin \mathrm{\left(}{\mathit{k}}_{\mathrm{x}}\mathit{x}\mathrm{\right)}\sin \mathrm{\left(}{\mathit{k}}_{\mathrm{z}}\mathit{z}\mathrm{\right)}\mathit{,}& & \end{array}$(8)with A the amplitude of the two-dimensional perturbation. This initial condition can be inspected in Fig. 6. The subsequent panels show the evolution of the inserted particle droplets. Initially they go through a phase of contraction without altering the amplitude. This can be seen in further detail in the vertical slices in Fig. 7. The droplets remain in terminal velocity. Intriguingly, the wave steepening is scale-independent, as can be seen in the right panel of Fig. 7, where we have decreased the size of the droplets by a factor of 2. Subsequently, the particle swarm transforms into a characteristic mushroom cloud, which is reminiscent of those seen in the standard Rayleigh-Taylor instability.

Fig. 7 Slices through the evolution of the particle droplets, along the z-direction, similar to Fig. 4. The droplets sediment at terminal velocity, and the wave steepening is independent of the initial droplet scale. The left panel shows run runEGG with λx = 1,λz = 4 and the right panel runEGG2 with λx = 0.5,λz = 2.

Figure 8 illustrates how nonlinear drafting results in particle concentrations. Initially, the droplets concentrate the material by collapsing onto themselves. However, in the subsequent evolution, it is clear that the overdense regions drag gas downwards with their sedimentation flow. Simultaneously, gas moves upwards in between the denser areas. This is an aspect of the imbalanced stratification that we discussed in Sect. 3.2.

3.5. Toy model of the drafting instability

The origin of the instability can be grasped from a simplified stability analysis, which is based on the observation that regions overdense in particles sediment faster than regions that are underdense in particles. Equations (5) and (6) describe the behavior of the particles, which only depends on the gas through the gas drag term. We now assume the gas velocity can be written as $\begin{array}{ccc}\mathit{u}\mathrm{=}\mathit{\alpha }\mathrm{(}\mathit{ϵ}\mathrm{-}{\mathit{ϵ}}_{\mathrm{0}}\mathrm{)}\mathit{v}\mathit{.}& & \end{array}$(9)Here, α is a proportionality constant that can be determined numerically and which contains the dependence of the instability on the viscosity and the Mach number. The quantities u and v are the vertical gas and particle velocity. Effectively, we assume that the linearized gas velocity δu can be expressed proportional to the particle density perturbation δu = α(v_f/ρ₀)δρ_p (see Appendix B for more details). Physically, it expresses the observation that the gas follows overdense particle regions but, locally, momentum is conserved by the gas becoming buoyant and moving in the opposite direction in underdense regions. This assumption allows us to reduce the equations to one dimension, even though the above approximation implicitly assumes two or three dimensions to be present to allow gas to move freely and not be trapped as in our stable one-dimensional experiment (Fig. 8). We will also approximate the gas density to be constant.

The equilibrium state corresponds to pure sedimentation with ϵ = ϵ₀ and the particles having the terminal velocity v = −gt_f. The dispersion relation for Fourier modes of the form ∝exp⁽ωt−ikx⁾ becomes $\begin{array}{ccc}{\mathit{\omega }}^{\mathrm{\prime }}\mathrm{=}\mathrm{+}\mathit{i}{\mathit{k}}^{\mathrm{\prime }}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\left(\mathrm{-}\mathrm{1}\mathrm{\pm }\sqrt{\mathrm{1}\mathrm{-}\mathrm{4}\mathit{\alpha }{\mathit{ϵ}}_{\mathrm{0}}{\mathit{k}}^{\mathrm{\prime }}\mathit{i}}\right)\mathit{,}& & \end{array}$(10)where ω′ is the growth amplitude, k′ = 2π/λ′ the wave number of wavelength λ′ in friction units. The last term of Eq. (10) is always positive² for ϵ₀> 0, resulting in the exponential growth of the instability. The fastest growing modes are those with the shortest wavelength. This result, although derived somewhat differently (see Appendix B), is identical to that of plate drag (Goodman & Pindor 2000; Chiang & Youdin 2010; Jacquet et al. 2011). The real part of the dispersion relation is illustrated in Fig. B.1.

Fig. 8 The particle and gas evolution around particle droplets, at times t = 5,10tf. Streamlines show the gas velocity, while the particle density is color coded as in Fig. 6. Results from runEGG.

For large k, corresponding to short wavelengths ($\mathit{\lambda }\mathrm{\ll }\mathrm{8}\mathit{\pi \alpha }{\mathit{ϵ}}_{\mathrm{0}}{\mathit{t}}_{\mathit{f}}^{\mathrm{2}}\mathit{g}$), the growth rate can be approximated by $\begin{array}{ccc}{\mathit{\omega }}_{\mathrm{grow}}^{\mathrm{\prime }}\mathrm{\approx }\frac{\mathrm{1}}{\sqrt{\mathrm{2}}}\sqrt{\mathit{\alpha }{\mathit{ϵ}}_{\mathrm{0}}{\mathit{k}}^{\mathrm{\prime }}}\mathit{,}& & \end{array}$(11)by series expansion to leading order. Alternatively, we can stop expressing this result in friction units but as $\begin{array}{ccc}{\mathit{\omega }}_{\mathrm{grow}}\mathrm{\approx }\sqrt{\mathit{\alpha \pi }\frac{{\mathit{ϵ}}_{\mathrm{0}}\mathit{g}}{\mathit{\lambda }}}\mathrm{·}& & \end{array}$(12)Interestingly, in this asymptotic limit case, the growth rate no longer depends on the particle size (or more accurately the friction time), but only on the spatial scale and the dust-to-gas ratio. On larger scales, the limit expression of the growth rate takes the form $\begin{array}{ccc}{\mathit{\omega }}_{\mathrm{grow}}^{\mathrm{\prime }}\mathrm{=}{\left(\mathit{\alpha }{\mathit{ϵ}}_{\mathrm{0}}{\mathit{k}}^{\mathrm{\prime }}\right)}^{\mathrm{2}}\mathit{,}& & \end{array}$(13)to leading order. Therefore, the growth rate of the drafting instability rapidly decreases with increased spatial scale. In this branch, the growth rate does depend on the particle size, $\begin{array}{ccc}{\mathit{\omega }}_{\mathrm{grow}}\mathrm{=}\mathrm{(}\mathit{\alpha ϵkg}{\mathrm{)}}^{\mathrm{2}}{\mathit{t}}_{\mathrm{f}}^{\mathrm{3}}\mathit{.}& & \end{array}$(14)

Fig. 9 Viscosity dependency on the growth rate of the drafting instability. Dashed lines give the growth time scaling where ω ∝ ν− 1/3, as in the toy model. At high viscosities over ν′ ≈ 10-3, we no longer identify growth in the fluid velocity dispersion. Numerical results from runv2, runv3, run2, runv5, and runv6.

3.6. Dependence on the Reynolds number

If the toy model is right, then the growth rate scales as ${\mathit{\omega }}_{\mathrm{grow}}^{\mathrm{\prime }}\mathrm{\propto }{\mathit{k}}^{{}^{\mathrm{\prime }}\mathrm{2}}$ on length scales above the characteristic scale: $\begin{array}{ccc}{\mathit{\lambda }}_{\mathrm{knee}}\mathrm{=}{\mathrm{2}}^{\mathrm{4}\mathit{/}\mathrm{3}}\mathit{\pi \alpha }{\mathit{ϵ}}_{\mathrm{0}}{\mathit{l}}_{\mathrm{f}}\mathit{,}& & \end{array}$(15)obtained from balancing the fast and slow growth branches (Eqs. (11) and (13)).

Viscous damping has a similar quadratic dependence on k′. Therefore a viscosity cut-off exists: at ν larger than ν_crit growth of the instability is terminated. We find the critical viscosity by equating the large-scale growth timescale (Eq. (13)) with the viscous timescale (λ^′2/ν′), $\begin{array}{ccc}{\mathit{\nu }}_{\mathrm{crit}}^{\mathrm{\prime }}\mathrm{=}\mathrm{2}\mathit{\pi }\mathrm{(}\mathit{\alpha ϵ}{\mathrm{)}}^{\mathrm{2}}\mathit{.}& & \end{array}$(16)The determination of α allows us to scale the growth rates of the toy model. From Fig. 9, we find, numerically, that above viscosities of around ν′ ~ 10^-3, the instability does not show up. A critical viscosity of ${\mathit{\nu }}_{\mathrm{crit}}^{\mathrm{\prime }}\mathrm{~}{\mathrm{10}}^{-3}$ would correspond to α ~ 10^-2. This type of estimate is an approximate agreement with the seen correlation between particle density variations ϵ and the gas fluid velocity in Fig. 3.

For viscosities below the viscosity cut-off on the small-scale branch, the largest growing wavenumber scales as $\begin{array}{ccc}{\mathit{k}}_{\mathrm{crit}}^{\mathrm{\prime }}\mathrm{=}{\left(\mathrm{2}{\mathit{\pi }}^{\mathrm{2}}\frac{\mathit{\alpha }{\mathit{ϵ}}_{\mathrm{0}}}{{\mathit{\nu }}^{\mathrm{\prime }\mathrm{2}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}& & \end{array}$(17)by setting the viscous timescale equal to small-scale growth rate (Eq. (11)). Therefore the growth rate scales with the viscosity as ω′ ∝ ν^′−1/3. This indeed agrees with the results shown in Fig. 9.

3.7. Dependence on initial dust-to-gas ratio

The main free parameter in our model is the initial dust-to-gas ratio, also called the metallicity, when setting up the equilibrium stratification ϵ₀. Evidently in the limit of negligible dust loading, we do not expect any dust clumping. We therefore study the dependency of the growth rate of the instability on a lower than unity initial dust-to-gas ratio (Fig. 10). To measure slower growth rates (and possibly the saturation of the instability), we need to extend the vertical domain, which we achieve by numerically scaling the system (runs run1-4.0.1, see Table 2).

We find that the instability does not vanish, even at a ten times reduced metallicity. The growth rate is slower, and there seems to be a longer dormant phase before particle concentrations settle in. The reason for this delay for the instability to kick in is unclear. Between ϵ₀ = 0.1 to ϵ₀ = 1 the growth rates approximately scale proportionally to $\sqrt{{\mathit{ϵ}}_{\mathrm{0}}}$, as expected from the toy model (Eq. (12)).

At even lower metallicities, ϵ₀ = 0.05, we do not recover the instability. Growth rates become too slow to identify any particle clumping in the simulation (run4.01).

Fig. 10 Long term evolution of the dust-to-gas ratio max(ϵ) (black), combined with the horizontal $\sqrt{\mathrm{⟨}{\mathit{u}}_{\mathit{x}}^{\mathrm{2}}\mathrm{⟩}}$ (red) and vertical $\sqrt{\mathrm{⟨}{\mathit{u}}_{\mathit{z}}^{\mathrm{2}}\mathrm{⟩}}$ (blue). The grey dashed lines represent the growth rate scaling with $\sqrt{{\mathit{ϵ}}_{\mathrm{0}}}$, as found from the toy model (Eq. (11)). The growth rate decreases with dust-to-gas ratio, and we are unable to measure growth rates below ϵ0 = 0.05. The small oscillations seen in vertical velocity dispersion for the low metallicity runs are the result of vertical waves, which are due to a slight imbalance in the initial condition in elongated domains. The figure is based on the simulations run1.01, run2.01, run3.01, and run4.01 (in order of decreasing metallicity).

4. Enhanced particle concentrations in protoplanetary discs

In this section, we rescale our simulation results to the context of the protoplanetary disc. Because the drafting effect seems to prefer higher dust-to-gas ratios than the percentage level initially expected in a protoplanetary disc, we consider particle settling in a disc that has already undergone some grain growth, resulting in an already partly settled particle midplane (Sects. 4.1 and 4.2). In mass loaded regions, the swarms created by the drafting instability may aid the formation of chondrules, which we explore in Sect. 4.3. Finally we comment on the relevance of the drafting instability in the possibly highly dust-enriched envelopes around protoplanets (Sect. 4.4).

4.1. Rescaling friction units

The friction units employed in our simulations can be readily rescaled to a protoplanetary disc setting, for a given particle size expressed in Stokes number $\begin{array}{ccc}{\mathit{\tau }}_{\mathrm{f}}\mathrm{=}{\mathit{t}}_{\mathrm{f}}{\mathrm{\Omega }}_{\mathrm{K}}\mathit{,}& & \end{array}$(18)where Ω_K is the Keplerian frequency (for the definition of the friction time t_f, see Eq. (3)). From this definition, a time t_f corresponds to a fraction of a Keplerian timescale, ${\mathit{t}}_{\mathrm{f}}\mathrm{=}{\mathrm{\Omega }}_{\mathrm{K}}^{-1}{\mathit{\tau }}_{\mathrm{f}}$. Similarly, the friction length l_f can be expressed as $\begin{array}{ccc}{\mathit{l}}_{\mathrm{f}}\mathrm{=}\mathit{g}{\mathit{t}}_{\mathrm{f}}^{\mathrm{2}}\mathrm{=}\mathit{z}{\mathit{\tau }}_{\mathrm{f}}^{\mathrm{2}}\mathit{,}& & \end{array}$(19)when the gravity is expressed as Ω²z. Here, z is the height above the midplane. The friction length in the minimum mass solar nebula (MMSN, Hayashi 1981) at the top of a particle layer of thickness H_p can be written as $\begin{array}{ccc}{\mathit{l}}_{\mathrm{f}}\mathrm{\approx }\mathrm{37}\mathrm{\times }{\mathrm{10}}^{\mathrm{3}}\left(\frac{{\mathit{H}}_{\mathrm{p}}\mathit{/}\mathit{H}}{\mathrm{0.1}}\right){\left(\frac{{\mathit{\tau }}_{\mathrm{f}}}{\mathrm{0.1}}\right)}^{\mathrm{2}}{\left(\frac{\mathit{r}}{\mathrm{5}\hspace{0.17em}\mathrm{AU}}\right)}^{\mathrm{5}\mathit{/}\mathrm{4}}\hspace{0.17em}\mathrm{km}\mathit{.}& & \end{array}$(20)This scale strongly depends on the particle size (τ_f = 0.1 corresponds to a 2 cm particle at an orbital distance of r ≈ 5 AU). Here, we have assumed that the particle scale height is a constant fraction of the gas scale height H.

The ratio of the terminal velocity to the sound speed, the Mach number $\begin{array}{ccc}\mathit{Ma}\mathrm{=}\frac{{\mathit{v}}_{\mathrm{f}}}{{\mathit{c}}_{\mathrm{s}}}\mathrm{=}\mathrm{0.01}\left(\frac{{\mathit{H}}_{\mathrm{p}}\mathit{/}\mathit{H}}{\mathrm{0.1}}\right)\left(\frac{{\mathit{\tau }}_{\mathrm{f}}}{\mathrm{0.1}}\right)\mathit{,}& & \end{array}$(21)reveals the incompressible nature of particle sedimentation.

We can ignore the overall rotation of the protoplanetary disc for the small scales that we consider here. The Rossby number Ro = v_f/ (Ω_Kl_f) takes the form: Ro ~ 1/(Ω_Kt_f), when using friction scales. Therefore, for particles with small Stokes number τ_f = Ω_Kt_f ≪ 1, rotation is not important, and the rotation-free assumption is valid.

The kinematic molecular viscosity depends on the gas mean free path λ in the midplane of the protoplanetary disc as $\begin{array}{ccc}\mathit{\nu }\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{c}}_{\mathrm{s}}\mathit{\lambda }\mathit{.}& & \end{array}$(22)The viscosity can then be expressed in friction units as $\begin{array}{ccc}\frac{\mathit{\nu }}{{\mathit{g}}^{\mathrm{2}}{\mathit{t}}_{\mathrm{f}}^{\mathrm{3}}}\mathrm{=}\mathrm{1.6}\mathrm{\times }{\mathrm{10}}^{-6}{\left(\frac{{\mathit{H}}_{\mathrm{p}}\mathit{/}\mathit{H}}{\mathrm{0.1}}\right)}^{-2}{\left(\frac{{\mathit{\tau }}_{\mathrm{f}}}{\mathrm{0.1}}\right)}^{-3}{\left(\frac{\mathit{r}}{\mathrm{5}\hspace{0.17em}\mathrm{AU}}\right)}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{·}& & \end{array}$(23)This value does not differ greatly from the nominal value probed in our numerical work ($\mathit{\nu }\mathit{/}\mathrm{\left(}{\mathit{g}}^{\mathrm{2}}{\mathit{t}}_{\mathrm{f}}^{\mathrm{3}}\mathrm{\right)}\mathrm{=}{\mathrm{10}}^{-4}$, see also the list of simulations in Table 1). The strong scaling with orbital radius becomes much weaker if one considers particles of constant radius, as opposed to constant Stokes number.

Finally, we also verify that the viscous particle coupling criterion, given by Eq. (4), holds in the MMSN, $\begin{array}{ccc}\frac{{\mathit{t}}_{\mathit{\nu }\mathit{,}\mathrm{pair}}}{{\mathit{t}}_{\mathrm{f}}}\mathrm{\approx }\mathrm{3.0}\mathrm{\times }{\mathrm{10}}^{-4}{\left(\frac{{\mathit{ϵ}}_{\mathrm{0}}}{\mathrm{0.1}}\right)}^{\mathrm{-}\mathrm{2}\mathit{/}\mathrm{3}}\left(\frac{{\mathit{\tau }}_{\mathrm{f}}}{\mathrm{0.1}}\right){\left(\frac{\mathit{r}}{\mathrm{5}\hspace{0.17em}\mathrm{AU}}\right)}^{\mathrm{-}\mathrm{31}\mathit{/}\mathrm{6}}\mathit{,}& & \end{array}$(24)where ϵ₀ is the approximate mean dust-to-gas ratio in the particle midplane.

4.2. Applying the toy model

With the help of the toy model, we can attempt to further constrain where in the protoplanetary disc the drafting instability can occur. Because growth rates decrease rapidly at large scales, we only expect the instability to take place on the small scale branch, below the characteristic scale λ_knee. From Eq. (15), we get $\begin{array}{ccc}{\mathit{\lambda }}_{\mathrm{knee}}& & \mathrm{\approx }\mathrm{290}\mathrm{(}{\frac{\mathit{\alpha }}{\mathrm{0.01}}}^{\mathrm{)}}\mathrm{(}{\frac{{\mathit{ϵ}}_{\mathrm{0}}}{\mathrm{0.1}}}^{\mathrm{)}}\left(\frac{{\mathit{H}}_{\mathrm{p}}\mathit{/}\mathit{H}}{\mathrm{0.1}}\right){\mathrm{(}{\frac{{\mathit{\tau }}_{\mathrm{f}}}{\mathrm{0.1}}}^{\mathrm{)}}}^{\mathrm{2}}{\mathrm{(}{\frac{\mathit{r}}{\mathrm{5}\hspace{0.17em}\mathrm{AU}}}^{\mathrm{)}}}^{\mathrm{5}\mathit{/}\mathrm{4}}\hspace{0.17em}\mathrm{km}\mathit{.}\end{array}$(25)The dust-to-gas ratio of ϵ₀ = 0.1 will be relevant in a midplane layer of solids with H_p/H = 0.1, when the overall metallicity of the protoplanetary disc is the canonical Z = 0.01. However, we have chosen to keep dust-to-gas ratio ϵ₀ and the height of the particle layer H_p as independent quantities, because we do not necessarily want to study the conditions of a particle layer that has settled to equilibrium.

Fig. 11 Relevant length scales for the drafting instability in the minimum mass solar nebula: the friction length (lf), the upper scale for fast growth (λknee), and the scale at which viscosity dominates (λvisc). We consider particles located at a particle scale height above the midplane, with Hp/H = 0.1 and the midplane dust-to-gas ratio is 0.1. We take the toy model parameter to be α = 0.01. Particles are assumed to be 1 cm in radius, or 1 mm in a gas depleted disc with 10 times lower gas surface density (in which case the curves remain the same, but viscous scale λvisc,d is now the red dotted line). The instability likely does not operate in the inner (<5 AU) of the protoplanetary disc, unless particles are large or significant pre-concentration occurs.

A lower limit on the scale of the instability is set by viscous damping of the instability at scales below the knee, $\begin{array}{ccc}{\mathit{\lambda }}_{\mathrm{visc}}& \mathrm{=}& {\left(\frac{{\mathrm{2}}^{\mathrm{2}}\mathit{\pi }}{\mathit{\alpha }{\mathit{ϵ}}_{\mathrm{0}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}{\left(\frac{\mathit{\nu }}{{\mathit{g}}^{\mathrm{2}}{\mathit{t}}_{\mathrm{f}}^{\mathrm{3}}}\right)}^{\mathrm{2}\mathit{/}\mathrm{3}}{\mathit{l}}_{\mathrm{f}}\\ & \mathrm{\approx }& \end{array}$(26)We note that this scale, as opposed to the friction length l_f (Eq. (20)), does not depend on the particle size.

In Fig. 11, we illustrate the different relevant scales presented in Eqs. (20), (25) and (26), as function of orbital radius. The instability would operate on a scale of the order of 10⁴ km in the outer parts of the protoplanetary disc for particles of cm in size, assuming an MMSN model. We note that Fig. 11 shows the scaling for an assumed constant particle size, as opposed to Eqs. (20)–(26), which assume a constant Stokes number.

4.3. Chondrules

Chondrules are mm-sized inclusions found in primitive meteorites, which originate from the asteroid belt. It is generally accepted that a chondrule is the product of a flash heating event. The exact nature of chondrule precursors is unknown. However the heating events likely occurred in particle swarms at least 100 to 1000 km wide, with a local number density of about ~10 m^-3. In this way, the loss of light isotopes (isotopic fractionation) is prevented by exchanging vapour from chondrule to chondrule (Cuzzi & Alexander 2006). This scenario requires local chondrule densities more than 100 times above a dust-to-gas ratio of unity. Even higher concentrations might be necessary to explain the retention of sodium (Alexander et al. 2008).

Such high chondrule densities are surprising, since small particles are hard to concentrate to the midplane. Even in the absence of other forms of turbulence, particles sediment to a midplane with a dust-to-gas ratio not higher than approximately unity, because of the stirring caused by the streaming instability (Bai & Stone 2010a). However, the isotopic constraints on the need to concentrate chondrules weaken if the gas at the chondrule formation sites had a non-solar composition. The atmospheres around planetary embryos have been proposed as such locations (Morris et al. 2012). Nevertheless, in this scenario, pre-clumping of solids by a factor of at least 10 over midplane densities remains necessary and the shock waves invoked to melt chondrules lead, in fact, to destructive collisions (Jacquet & Thompson 2014).

Small particles are difficult to concentrate in the inner protoplanetary disc because of the strong sensitivity of the preferential scale of the instability on particle size (Eqs. (20) and (25)), as can be seen in Fig. 11. Nevertheless the connection to chondrule formation is tantalizing, especially because if clumping conditions are met, the drafting effect only weakly depends on particle size and efficiently clusters particles down to very small sizes (Eq. (12)). This is different from, for example, the streaming instability that has a preferred particle size, somewhat above that of chondrules for nominal metallicities (Carrera et al. 2015).

The drafting instability could operate on such small scales if some form of pre-concentration of solids were to occur. These types of enhanced particle densities could possibly occur near the Kolmogorov scale of the disc turbulence (Cuzzi et al. 2001). Alternatively, particle concentrations near sublimation lines can dramatically peak (Ros & Johansen 2013). An increase in the dust-to-gas ratio can also occur by the accretion of gas onto the star, which depletes the disc relative to the MMSN (Bitsch et al. 2015). Alternatively, growth rates could be increased if the unknown chondrule precursors are much larger than the chondrules they are turned into after the heating event. Even so, it remains to be seen if drafting instabilities can push particle concentrations to the desired high levels, even in such favourable instances.

4.4. Planetary atmospheres

The drafting instability might be important in the atmospheres of giant planets. The opacity in the outer envelope, which regulates the transport of heat, comes from the dust component. Under standard assumed opacities, it is difficult to cool the envelope and trigger runaway gas accretion (Ikoma et al. 2000; Piso & Youdin 2014). However, clumping of solids and the growth of the accreted dust could significantly reduce the opacity in the upper atmosphere. The friction length for particles sedimenting in a planetary atmosphere is given by $\begin{array}{ccc}{\mathit{l}}_{\mathrm{f}\mathit{,}\mathrm{plan}}& \mathrm{=}& \frac{\mathit{GM}}{{\mathit{r}}_{\mathrm{B}}^{\mathrm{2}}}{\mathit{t}}_{\mathrm{f}}^{\mathrm{2}}\\ & \mathrm{\approx }& \end{array}$(27)where ${\mathit{r}}_{\mathrm{B}}\mathrm{=}\mathit{GM}\mathit{/}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}$ is the thermal Bondi radius of a planet with mass M, which corresponds to the outer edge of the atmosphere. Here, we have considered particles on top of the envelope, but deeper in the planet the friction time shrinks owing to the increase in density. Applying the toy model, we estimate the knee scale in the upper envelope at $\begin{array}{ccc}{\mathit{\lambda }}_{\mathrm{knee}\mathit{,}\mathrm{plan}}\mathrm{\approx }\mathrm{570}\left(\frac{\mathit{\alpha }}{\mathrm{0.01}}\right)\left(\frac{{\mathit{ϵ}}_{\mathrm{0}}}{\mathrm{1}}\right){\left(\frac{\mathit{R}}{\mathrm{1}\hspace{0.17em}\mathrm{mm}}\right)}^{\mathrm{2}}{\left(\frac{\mathit{M}}{\mathrm{5}\hspace{0.17em}{\mathrm{M}}_{\mathrm{E}}}\right)}^{-1}{\left(\frac{\mathit{r}}{\mathrm{5}\hspace{0.17em}\mathrm{AU}}\right)}^{\mathrm{5}}\hspace{0.17em}\mathrm{km}\mathit{,}& & \end{array}$(28)which is above the damping viscosity scale at $\begin{array}{ccc}{\mathit{\lambda }}_{\mathrm{visc}\mathit{,}\mathrm{plan}}\mathrm{\approx }\mathrm{13}{\left(\frac{\mathit{\alpha }}{\mathrm{0.01}}\right)}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{3}}{\left(\frac{{\mathit{ϵ}}_{\mathrm{0}}}{\mathrm{1}}\right)}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{3}}{\left(\frac{\mathit{M}}{\mathrm{5}\hspace{0.17em}{\mathrm{M}}_{\mathrm{E}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}{\left(\frac{\mathit{r}}{\mathrm{5}\hspace{0.17em}\mathrm{AU}}\right)}^{\mathrm{2}}\hspace{0.17em}\mathrm{km}\mathit{.}& & \end{array}$(29)The formation of ice giants and super-Earths might be paired with significant amounts of dust in their low-mass gaseous envelopes (Lee et al. 2014). These scaling relations argue that order-of-unity mass loading of atmospheres will lead to clumping and the breakup of the dust component, providing an upper limit on the dust opacity.

5. Future outlook

5.1. Laboratory experiments

This study supports ongoing work to investigate drag instabilities in the laboratory. A full description of the apparatus constructed at the Max Planck Institute for Dynamics and Self-Organisation and first results will be presented in an accompanying paper (Capelo et al., in prep.). Here, we restrict ourselves to highlight some aspects relevant to the understanding of the numerical results presented in this paper.

The experimental apparatus consists of a cylindrical vessel, housing a gas stream operating at pressures in the range of 0.5−10³ millibar. The axial component of the cylindrical flow is parallel with the direction of Earth’s gravity, similar to the sedimentation configuration in the simulations presented here. The upwards steady-state flow will be seeded with weakly inertial particles with typical sizes of 20–90μm. The range of operational pressures and temperatures, listed in Table 3, then enable us to span both the Stokes and Epstein drag regimes.

The particle entrainment happens upstream in the fluid flow. There, the system is in a brief transient state. The solids are transported and mixed with the gas by the time the flow reaches the steady-state conditions in which the measurements are to be made. This is performed to make a fair comparison with the nearly homogeneously mixed initial conditions of the two-fluid dust/gas models.

Table 3 summarizes the parameter region in which the experiment operates, including gas state variables, Mach and Reynolds numbers. The flow conditions are incompressible and laminar. The experiment will be the first in its kind probing the Epstein drag regime in a fluid with equal mass loading of gas and particles.

The experiment described here is somewhat analogous to particle suspension experiments in Newtonian fluids with low particle Reynolds number (Guazzelli & Hinch 2011). However, in these studies volume fractions, φ = n_pR³ (with n_p and R the particle number density and radius), are no lower than φ ≈ 0.01%. Our experiment operates at φ ≈ 10^-4%, when the dust-to-gas ratio is unity for solid spherical particles with densities ranging from that of vitreous carbon (ρ_• = 1.4 g cm^-3) to steel (ρ_• = 8 g cm^-3). The low particle Reynolds number, Re_p, in these suspension experiments comes from the use of a fluid with high dynamic viscosity. The particles are very buoyant and slowly creep through a thick liquid. Here, on the other hand, the low values of Re_p come from the fact that the kinematic viscosity becomes high when the gas density is low. It is encouraging that such experiments, even if in a regime that is different to the one studied here, show interesting particle dynamics (Batchelor 1972), such as particle Rayleigh-Taylor mushrooms and drafting particle trains (Pignatel et al. 2009; Matas et al. 2004).

Time-resolved data on the particle trajectories will be obtained from high-resolution cameras and 3D Lagrangian particle tracking (Xu 2008; Ouellette et al. 2006). This is a common technique to study both tracer and intertial particles in fluids. The typical measured and derived quantities are the probability density distributions of the particle velocities and accelerations, their statistical moments, and correlation and structure functions.

The obtained data will provide an interesting comparison to the results shown in this paper. The parameter regime is sufficiently similar to the numerical experiments that we expect the drafting instability to manifest itself. Particle tracking would not only enable the detection of particle swarms, but also the individual particle dynamics. For instance, Fig. 2 demonstrates that the growing maximum in particle velocity dispersion traces the increase in maximum particle density. These types of statistical measurements of the particles will allow a qualitative comparison between the numerical work and the experiments.

5.2. Numerical work

Here, we have presented several numeric experiments to demonstrate a drafting instability. Future work will refine the estimates made in this paper. For example, the numerical set up is currently limited to studying sedimentation on rather short timescales, set by the length of the simulation domain. This could be avoided in future work by implementing a form of periodic boundary conditions in the vertical direction, which would recycle particles.

To aid the interpretation of the experimental results, it will be necessary to specifically reproduce the parameter regime in which the apparatus operates. Additionally, refined boundary conditions will be needed to approximate the experimental set-up. This work is under progress, but evidently awaits the first experiments.

We have also argued that the drafting instability could be of relevance in a protoplanetary disc. To study this connection in more detail, it will be necessary to simulate numerically expensive larger domains that encompass the disc midplane. Additionally, the connection between the drafting and streaming instability could be studied in more detail. Ultimately, the results should be placed in the context of other sources of disc turbulence, such as the magnetorotational instability that operates in sufficiently ionized regions or in the penetration of vertical shear instability to the midplane (Turner et al. 2014). Additionally, the growth of particles through coagulation or condensation will need to be taken into account in a self-consistent matter. Future work is needed to understand the possibly constructive interplay of these mechanisms.

6. Summary

In this paper we have demonstrated the presence of a drafting instability when particles sediment through a fluid in hydrostatic balance. On timescales of tens of friction times, particles unmix out of a homogeneous mixture and particle concentrations increase by a factor 10.

The presence of such an instability was not expected because it evades detection in an analytic linear stability analysis. However, our numerical results demonstrate that the system is nonlinearly unstable. The exact nature of the instability is difficult to determine. We interpret the instability to be the result of an imbalance in the stratification, which locally disturbs the hydrostatic balance. We support this hypothesis with a simple toy model that captures some of the main characteristics of the instability. Growth is fastest at the smallest available scales, increases with the square root of the dust-to-gas ratio and a critical scale is identified at which viscosity overwhelms the instability.

By expressing our numerical results in a system of friction units, we can exploit our results by scaling them to either upcoming laboratory experiments or the protoplanetary disc. We argue that an experiment can achieve a dust-to-gas ratio of around unity and probe the regime of interest. In protoplanetary discs the drafting instability may take place in particle-rich layers above the midplane in the outer regions of the disc, on scales smaller than previously studied. In these regions, where the conditions for the drafting instability are met, we have shown that sedimenting particles spontaneously form dense clumps. Future work will be needed to investigate to what degree this clumping affects coagulation rates and whether the drafting instability can create the dense environment necessary for chondrule flash heating.

Acknowledgments

We would like to thank Haitao Xu, Wlad Lyra, Chao-Chin Yang, Karsten Dittrich, Hubert Klahr, Neal Turner, Mario Flock, and Andrew Youdin for valuable comments. The paper benefited greatly from insightful comments from two referees who helped clarify the stability analysis and the interpretation of the laboratory constraints on chondrule formation densities. Part of the computations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at LUNARC (Lund University). M.L. acknowledges funding from the Knut and Alice Wallenberg Foundation. A.J. acknowledges financial support from the European Research Council (ERC Starting Grant 278675- PEBBLE2PLANET), the Swedish Research Council (grant 2010-3710) and the Knut and Alice Wallenberg Foundation (grant Bottlenecks for particle growth in turbulent aerosols Dnr. KAW 2014.0048). H.C., J.B. and E.B. acknowledge support by the Deutsche Forschungsgemeinschaft under the grant INST 186/959-1 as part of the CRC 963 Astrophysical Flow Instabilities and Turbulence.

References

Appendix A: Linear stability analysis

We briefly rederive the stability analysis for a particle gas-mixture in a non-rotating flow (Youdin & Goodman 2005; Jacquet et al. 2011). We demonstrate the result in two spatial dimensions, but our conclusions remain valid when generalized to three dimensions.

Appendix A.1: Governing equations

We assume the gas to be incompressible, in line with Youdin & Goodman (2005), $\begin{array}{ccc}\frac{\mathit{\partial }{\mathit{u}}_{\mathit{x}}}{\mathit{\partial x}}\mathrm{+}\frac{\mathit{\partial }{\mathit{u}}_{\mathit{z}}}{\mathit{\partial z}}\mathrm{=}\mathrm{0}\mathit{,}& & \end{array}$(A.1)and use the standard momentum equations $\begin{array}{ccc}& & \frac{\mathit{\partial }{\mathit{u}}_{\mathit{x}}}{\mathit{\partial t}}\mathrm{+}{\mathit{u}}_{\mathit{x}}\frac{\mathit{\partial }{\mathit{u}}_{\mathit{x}}}{\mathit{\partial x}}\mathrm{+}{\mathit{u}}_{\mathit{z}}\frac{\mathit{\partial }{\mathit{u}}_{\mathit{x}}}{\mathit{\partial z}}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathit{\rho }}\frac{\mathit{\partial P}}{\mathit{\partial x}}\mathrm{+}\frac{\mathrm{1}}{{\mathit{t}}_{\mathrm{f}}}\frac{{\mathit{\rho }}_{\mathrm{p}}}{\mathit{\rho }}\mathrm{(}{\mathit{v}}_{\mathit{x}}\mathrm{-}{\mathit{u}}_{\mathit{x}}\mathrm{)}\mathit{,}\\ & & \end{array}$Similarly, for the particle fluid we make use of the continuity equation, $\begin{array}{ccc}\frac{\mathit{\partial }{\mathit{\rho }}_{\mathrm{p}}}{\mathit{\partial t}}\mathrm{+}\frac{\mathit{\partial }}{\mathit{\partial x}}\mathrm{(}{\mathit{\rho }}_{\mathrm{p}}{\mathit{v}}_{\mathit{x}}\mathrm{)}\mathrm{+}\frac{\mathit{\partial }}{\mathit{\partial z}}\mathrm{(}{\mathit{\rho }}_{\mathrm{p}}{\mathit{v}}_{\mathit{z}}\mathrm{)}\mathrm{=}\mathrm{0}\mathit{,}& & \end{array}$(A.4)and the set of momentum equations $\begin{array}{ccc}& & \frac{\mathit{\partial }{\mathit{v}}_{\mathit{x}}}{\mathit{\partial t}}\mathrm{+}{\mathit{v}}_{\mathit{x}}\frac{\mathit{\partial }{\mathit{v}}_{\mathit{x}}}{\mathit{\partial x}}\mathrm{+}{\mathit{v}}_{\mathit{z}}\frac{\mathit{\partial }{\mathit{v}}_{\mathit{x}}}{\mathit{\partial z}}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{{\mathit{t}}_{\mathrm{f}}}\mathrm{(}{\mathit{v}}_{\mathit{x}}\mathrm{-}{\mathit{u}}_{\mathit{x}}\mathrm{)}\mathit{,}\\ & & \end{array}$This completes the model with six parameters (ρ_p, ρ,v_x, v_z, u_x, u_z), and as many equations.

Appendix A.2: Equilibrium solution

In equilibrium, we have no vertical motion (u_x = 0,v_x = 0). Additionally, we assume the gas to be in the rest frame, u_z = 0, and particles to be uniformly spread (ρ_p = ρ_p,0 constant) initially. This leaves the particle continuity and z-momentum equations as the non-trivial equations that determine v_z and ρ,ρ_p, $\begin{array}{ccc}& & \frac{\mathit{\partial }}{\mathit{\partial z}}\mathrm{(}{\mathit{\rho }}_{\mathrm{p}}{\mathit{v}}_{\mathit{z}}\mathrm{)}\mathrm{=}\mathrm{0}\mathit{,}\\ & & {\mathit{v}}_{\mathit{z}}\frac{\mathit{\partial }{\mathit{v}}_{\mathit{z}}}{\mathit{\partial z}}\mathrm{=}\mathrm{-}\mathit{g}\mathrm{-}\frac{\mathrm{1}}{{\mathit{t}}_{\mathrm{f}}}{\mathit{v}}_{\mathit{z}}\mathit{,}\\ & & \end{array}$In the last equation, we have assumed an isothermal gas, $\mathit{P}\mathrm{=}\mathit{\rho }{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}$. The equilibrium solution then takes the form $\begin{array}{ccc}& & {\mathit{v}}_{\mathit{z}}\mathrm{=}\mathrm{-}\mathit{g}{\mathit{t}}_{\mathrm{f}}\mathit{,}\\ & & \end{array}$where ρ_b is the gas density at the z = 0 boundary.

Appendix A.3: Dispersion relation

We now consider a first order perturbation of this equilibrium state. For the gas we find $\begin{array}{ccc}& & \\ & & \\ & & \end{array}$For the particles we get $\begin{array}{ccc}& & \\ & & \frac{\mathit{\partial }{\mathit{v}}_{\mathit{x}}^{\mathrm{\prime }}}{\mathit{\partial t}}\mathrm{+}{\mathit{v}}_{\mathit{z,}\mathrm{0}}\frac{\mathit{\partial }{\mathit{v}}_{\mathit{x}}^{\mathrm{\prime }}}{\mathit{\partial z}}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{{\mathit{t}}_{\mathrm{f}}}\left({\mathit{v}}_{\mathit{x}}^{\mathrm{\prime }}\mathrm{-}{\mathit{u}}_{\mathit{x}}^{\mathrm{\prime }}\right)\mathit{,}\\ & & \end{array}$For modes of the form A′ ∝ exp(ωt−ik_x−ik_z), we find that the system only has non-zero solutions when the determinant is zero, $\begin{array}{ccc}& & \left(\mathit{\omega }\mathrm{-}\mathrm{i}{\mathit{v}}_{\mathit{z}}{\mathit{k}}_{\mathit{z}}\right)\mathrm{(}\mathit{\omega }\mathrm{-}\mathrm{i}{\mathit{v}}_{\mathit{z}}{\mathit{k}}_{\mathit{z}}\mathrm{+}{{\mathit{t}}_{\mathrm{f}}^{-1}}^{\mathrm{)}}\\ & & \end{array}$(A.18)The first term represents a pure particle mode falling at the terminal velocity, while the second term represents particle motion damped by gas drag. The last factor of this expression has the solutions $\begin{array}{ccc}\mathit{\omega }\mathrm{=}\frac{{\mathit{v}}_{\mathit{z}}{\mathit{k}}_{\mathit{z}}}{\mathrm{2}}\mathrm{i}\mathrm{+}\frac{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{ϵ}\mathrm{)}{\mathit{t}}_{\mathrm{f}}^{-1}}{\mathrm{2}}\left(\mathrm{-}\mathrm{1}\mathrm{\pm }\sqrt{\mathrm{1}\mathrm{-}\frac{{\mathit{k}}_{\mathit{z}}^{\mathrm{2}}{\mathit{v}}_{\mathit{z,}\mathrm{0}}^{\mathrm{2}}{\mathit{t}}_{\mathrm{f}}^{\mathrm{2}}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{ϵ}{\mathrm{)}}^{\mathrm{2}}}\mathrm{+}\mathrm{2}\frac{\mathrm{(}\mathit{ϵ}\mathrm{-}\mathrm{1}\mathrm{)}{\mathit{k}}_{\mathit{z}}{\mathit{v}}_{\mathit{z,}\mathrm{0}}{\mathit{t}}_{\mathrm{f}}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{ϵ}{\mathrm{)}}^{\mathrm{2}}}\mathit{i}}\right)\mathit{.}& & \end{array}$(A.19)The real part of the square root term is always below unity ($\mathrm{Re}{}^{\mathrm{\right(}}\sqrt{\mathrm{1}\mathrm{-}\mathrm{i}\mathit{x}}^{\mathrm{\left)}}\mathit{<}\mathrm{1}$ for any x), so the perturbation is damped.

In summary, this analysis shows that there are no growing modes under the assumption of incompressibility. Possibly, unstable modes could be found when relaxing the assumption of incompressibility, or by using more realistic equations of state or by exploring non-local perturbation techniques. We leave this for future work, given the complexity of such investigations.

Appendix B: Toy model dispersion relation

Fig. B.1 Linear behavior of the toy model. In black, the real solution to the dispersion relation is given (Eq. (10), with α = 0.01). The knee, the largest scale for fast growth, is located at λknee ≈ 0.08. The gray curve is the result for a dust-to-gas ratio of 0.1 as opposed to unity. The dashed red line gives the high k approximation (Eq. (11)) and the dashed blue line gives the low k approximation (Eq. (13)).

We start with making the ansatz that $\begin{array}{ccc}\mathit{u}\mathrm{=}\mathit{\alpha }\mathrm{(}\mathit{ϵ}\mathrm{-}{\mathit{ϵ}}_{\mathrm{0}}\mathrm{)}\mathit{v,}& & \end{array}$(B.1)which removes the explicit dependency on the equations for the gas component. Here α is a proportionality parameter that encapsulates the viscosity dependency and remains to be determined through numerical simulations³. For simplicity, we also assume a constant gas density, ρ = ρ₀. Subsequently, the drag term in the particle momentum equation takes the form $\begin{array}{ccc}\mathrm{-}\frac{\mathrm{1}}{{\mathit{t}}_{\mathrm{f}}}\mathrm{(}\mathit{v}\mathrm{-}{\mathit{u}}^{\mathrm{)}}\mathrm{=}\mathrm{-}\frac{\mathit{v}}{{\mathit{t}}_{\mathrm{f}}}\left[\mathrm{1}\mathrm{-}\mathit{\alpha }\mathrm{(}\mathit{ϵ}\mathrm{-}{\mathit{ϵ}}_{\mathrm{0}}\mathrm{)}\right]\mathit{.}& & \end{array}$(B.2)In equilibrium, the dust-to-gas ratio, ϵ = ρ_p/ρ, is constant. The momentum equation then shows that particles move, as desired, with terminal velocity, $\begin{array}{ccc}{\mathit{v}}_{\mathrm{0}}\mathrm{=}\mathrm{-}\mathit{g}{\mathit{t}}_{\mathrm{f}}\mathit{.}& & \end{array}$(B.3)We now linearize the system of particle equations $\begin{array}{ccc}& & {\mathit{\partial }}_{\mathit{t}}{\mathit{\rho }}_{\mathrm{p}}^{\mathrm{\prime }}\mathrm{+}{\mathit{\rho }}_{\mathrm{p}\mathit{,}\mathrm{0}}{\mathit{\partial }}_{\mathit{z}}{\mathit{v}}^{\mathrm{\prime }}\mathrm{+}{\mathit{v}}_{\mathrm{0}}{\mathit{\partial }}_{\mathit{z}}{\mathit{\rho }}_{\mathrm{p}}^{\mathrm{\prime }}\mathrm{=}\mathrm{0}\mathit{,}\\ & & {\mathit{\partial }}_{\mathit{t}}{\mathit{v}}^{\mathrm{\prime }}\mathrm{+}{\mathit{v}}_{\mathrm{0}}{\mathit{\partial }}_{\mathit{z}}{\mathit{v}}^{\mathrm{\prime }}\mathrm{=}\mathrm{-}\mathit{g}\mathrm{-}\frac{{\mathit{v}}_{\mathrm{0}}\mathrm{+}{\mathit{v}}^{\mathrm{\prime }}}{{\mathit{t}}_{\mathrm{f}}}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\alpha }{\mathit{ϵ}}^{\mathrm{\prime }}\mathrm{)}\mathit{.}\end{array}$The last two terms simplify⁴ to $\begin{array}{ccc}\mathrm{-}\mathit{g}\mathrm{-}\frac{{\mathit{v}}_{\mathrm{0}}\mathrm{+}{\mathit{v}}^{\mathrm{\prime }}}{{\mathit{t}}_{\mathrm{f}}}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\alpha }{\mathit{ϵ}}^{\mathrm{\prime }}\mathrm{)}\mathrm{=}\mathrm{-}\frac{{\mathit{v}}^{\mathrm{\prime }}}{{\mathit{t}}_{\mathrm{f}}}\mathrm{-}\mathit{g\alpha }\frac{{\mathit{\rho }}_{\mathrm{p}}^{\mathrm{\prime }}}{{\mathit{\rho }}_{\mathrm{0}}}\mathrm{·}& & \end{array}$(B.6)By taking now modes of the form A′ ∝ exp(ωt−ikz), we are left with the following system of equations: $\begin{array}{ccc}\left(\begin{array}{c}\\ & & \\ & & \end{array}\right)\left(\begin{array}{c}\\ \\ \end{array}\right)\mathrm{=}\left(\begin{array}{c}\\ \\ \end{array}\right)\mathit{.}& & \end{array}$(B.7)Non-zero solutions are found when $\begin{array}{ccc}{\mathit{\beta }}^{\mathrm{2}}\mathrm{+}\frac{\mathit{\beta }}{{\mathit{t}}_{\mathrm{f}}}\mathrm{+}\mathrm{i}\mathit{\alpha }{\mathit{ϵ}}_{\mathrm{0}}\mathit{gk}\mathrm{=}\mathrm{0}\mathit{,}& & \end{array}$(B.8)where β = ω−ikv₀. Thus we find $\begin{array}{ccc}\mathit{\beta }\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}{\mathit{t}}_{\mathrm{f}}}\left(\mathrm{-}\mathrm{1}\mathrm{\pm }\sqrt{\mathrm{1}\mathrm{-}\mathrm{4}\mathit{\alpha }{\mathit{ϵ}}_{\mathrm{0}}\mathit{g}{\mathit{t}}_{\mathrm{f}}^{\mathrm{2}}\mathit{k}\mathrm{i}}\right)\mathit{,}& & \end{array}$(B.9)where the last term has a positive real part larger than 1 for any product $\mathit{\alpha }{\mathit{ϵ}}_{\mathrm{0}}\mathit{g}{\mathit{t}}_{\mathrm{f}}^{\mathrm{2}}$ that is different from 0. This reproduces Eq. (10), allowing the approximation of the two limit cases, Eq. (11) at small scales and Eq. (13) at large scales. The shape of the dispersion relation and the two limit cases can be inspected in Fig. B.1.

Appendix C: Particle number test

Particle numbers of 16 superparticles per grid cell are sufficient to correctly capture the evolution of the particle-gas mixture. However, increased particle numbers decrease the noise amplitude that is initially injected. In Fig. C.1, we show the evolution of the maximal particle density and gas velocity dispersion, as a function of the particle number. Because of the nonlinear nature of the drafting instability, we can see that the decreased noise amplitude with increased particle number prolongs a dormant state before the instability comes fully into effect and growth rates decrease moderately. However, if we ignore the protracted dormant phase, growth rates between 16 (our nominal value) and 64 particles per grid cell are undistinguishable, although slower than the 4 particles per grid cell case.

Fig. C.1 Evolution of the maximal particle density (top) and the horizontal gas dispersion (bottom), for 4, 16, and 32 particles per grid cell (respectively red, black, and blue curves). The initial maximal particle overdensity is reduced from ϵ/ϵ0 = 2.4 to 1.6 and finally 1.3, when increasing the particle number by a factor of 4 each time. The lowest particle number simulation shows fastest growth. In higher particle number runs, the initial dormant phase persists longer and growth rates become lower. Results from run2.n4, run2, and run2.n64.

All Tables

All Figures

Fig. 1 Development of particle swarms by a drafting instability. Displayed is the evolution of the particle density in the two-dimensional simulation run2 (see Table 1). The left-most panel illustrates the initial conditions: a stratified gas column in the vertical directions with particles sedimenting at terminal velocity, placed randomly throughout the simulation domain (1 lf wide and 20 lf high, note that the figure aspect ratio is enlarged in the x-direction). In the other panels (time t = 10,15 tf), regions marginally overdense with particles locally break the stratification equilibrium and accelerate downwards while in particle-poor regions, a deceleration from terminal velocity occurs. This leads to continued particle pileups through drafting, originating from gas being dragged by the particles. The particle and fluid components unmix and remain in this state. At the end of the simulation, when most particles have sedimented out of the simulation domain, the maximal particle density has increased by a factor 10.

In the text

Fig. 2 Top panel: evolution of the particle density (normalized to the background dust-to-gas ratio ϵ0), for three different resolutions in 2D simulations and a 3D simulation (run3d.4, gray). Bottom panel: evolution of the instability based on the horizontal gas dispersion, $\sqrt{\mathrm{⟨}{\mathit{u}}_{\mathit{x}}^{\mathrm{2}}\mathrm{⟩}}$, (as an alternative tracer to the maximal particle density, which is an intrinsically noisy variable). We used a similar color coding as for the top panel. Note for the runs with ϵ0 = 4, we displaced the dashed curves for clarity with a factor 10 downwards.

In the text

Fig. 3 Correlation between the local dust-to-gas ratio ϵ and the deviations in the velocity of the gas uz (full line) and particles vz (dotted line). To aid interpretation, we subtracted the mean velocity of the sedimenting particles ⟨vz⟩ from the particle velocity and a small artificial upwards mean gas velocity ⟨uz⟩ from the gas velocity. Displayed are different times t = 0,5,10 tf, in respectively black, blue, and red, based on the region between z = 0–10 lf in run2. The standard deviations on the binned averages are relatively large, σ ≈ 0.03–0.04 vf for the gas and particle velocities between t = 5–10 tf. The black dashed line corresponds to α = 0.02 in the toy model (Sect. 3.5).

In the text

Fig. 4 Left panel: sedimenting particles with a kz = π/ 2-mode in the particle density. Represented is a snapshot of runRT at t = 5tf. No instabilities develop, suggesting that the drafting instability is not directly related to the Rayleigh-Taylor instabilities that should occur in this set up. Right panel: evolution of the sedimenting particle wave at an arbitrary placed slice at x = −0.23. Different curves represent different times (t = 0,1,2,...,6tf), that for clarity are offset by 0.2. The particle density ρp(z) is given in black, while the red dashed line gives the gas velocity uz(z). Particle ridges do not approach each other, as indicated by the blue line that tracks the position of a point advected with velocity vf. The pressure that is supporting the mass-loaded stratification adapts to compensate for small changes in the particle density. This allows uz(z) to remain zero while the particles sediment.

In the text

Fig. 5 Sedimenting particles with a kx = π/ 2-mode in the particle density. We show different details (left) from a larger simulated domain (right), at different times (t = 1,6,9,12tf), based on runKH. Initially, we added noise in the particle velocity to perturb the boundary between particle-rich and particle-poor regions. We believe that we see a “parasitic” Kelvin-Helmholtz instability appear over time.

In the text

	Fig. 6 The evolution of particle droplets, at different times (t = 0,5,10tf), resulting in the emergence of characteristic Rayleigh-Taylor mushroom clouds. The color bar gives the color scale for the shown particle density. Results from runEGG.
In the text

	Fig. 7 Slices through the evolution of the particle droplets, along the z-direction, similar to Fig. 4. The droplets sediment at terminal velocity, and the wave steepening is independent of the initial droplet scale. The left panel shows run runEGG with λx = 1,λz = 4 and the right panel runEGG2 with λx = 0.5,λz = 2.
In the text

	Fig. 8 The particle and gas evolution around particle droplets, at times t = 5,10tf. Streamlines show the gas velocity, while the particle density is color coded as in Fig. 6. Results from runEGG.
In the text

	Fig. 9 Viscosity dependency on the growth rate of the drafting instability. Dashed lines give the growth time scaling where ω ∝ ν− 1/3, as in the toy model. At high viscosities over ν′ ≈ 10-3, we no longer identify growth in the fluid velocity dispersion. Numerical results from runv2, runv3, run2, runv5, and runv6.
In the text

Fig. 10 Long term evolution of the dust-to-gas ratio max(ϵ) (black), combined with the horizontal $\sqrt{\mathrm{⟨}{\mathit{u}}_{\mathit{x}}^{\mathrm{2}}\mathrm{⟩}}$ (red) and vertical $\sqrt{\mathrm{⟨}{\mathit{u}}_{\mathit{z}}^{\mathrm{2}}\mathrm{⟩}}$ (blue). The grey dashed lines represent the growth rate scaling with $\sqrt{{\mathit{ϵ}}_{\mathrm{0}}}$, as found from the toy model (Eq. (11)). The growth rate decreases with dust-to-gas ratio, and we are unable to measure growth rates below ϵ0 = 0.05. The small oscillations seen in vertical velocity dispersion for the low metallicity runs are the result of vertical waves, which are due to a slight imbalance in the initial condition in elongated domains. The figure is based on the simulations run1.01, run2.01, run3.01, and run4.01 (in order of decreasing metallicity).

In the text

Fig. 11 Relevant length scales for the drafting instability in the minimum mass solar nebula: the friction length (lf), the upper scale for fast growth (λknee), and the scale at which viscosity dominates (λvisc). We consider particles located at a particle scale height above the midplane, with Hp/H = 0.1 and the midplane dust-to-gas ratio is 0.1. We take the toy model parameter to be α = 0.01. Particles are assumed to be 1 cm in radius, or 1 mm in a gas depleted disc with 10 times lower gas surface density (in which case the curves remain the same, but viscous scale λvisc,d is now the red dotted line). The instability likely does not operate in the inner (<5 AU) of the protoplanetary disc, unless particles are large or significant pre-concentration occurs.

In the text

Fig. B.1 Linear behavior of the toy model. In black, the real solution to the dispersion relation is given (Eq. (10), with α = 0.01). The knee, the largest scale for fast growth, is located at λknee ≈ 0.08. The gray curve is the result for a dust-to-gas ratio of 0.1 as opposed to unity. The dashed red line gives the high k approximation (Eq. (11)) and the dashed blue line gives the low k approximation (Eq. (13)).

In the text

Fig. C.1 Evolution of the maximal particle density (top) and the horizontal gas dispersion (bottom), for 4, 16, and 32 particles per grid cell (respectively red, black, and blue curves). The initial maximal particle overdensity is reduced from ϵ/ϵ0 = 2.4 to 1.6 and finally 1.3, when increasing the particle number by a factor of 4 each time. The lowest particle number simulation shows fastest growth. In higher particle number runs, the initial dormant phase persists longer and growth rates become lower. Results from run2.n4, run2, and run2.n64.

In the text