© ESO 2019

1. Introduction

Planet formation theories that are consistent with the architecture, elemental composition, and collisional record in our solar system require a population of approximately km-sized precursor bodies, termed planetesimals. Due to primarily gravitational interactions with one another and with the surrounding gas, dust, and pebbles, planetesimals, which undergo repeated collisions, grow into planetary cores and ultimately mature planets possessing atmospheres (Morbidelli et al. 2012). A similar process must be responsible for the planets known to exist in extra-solar systems (Seager & Lissauer 2010).

How first-generation planetesimals form is under debate, with fundamental agreement that the only available material to compose them are the solids (“particles”) observed in gaseous protoplanetary discs (PPDs; Testi et al. 2014) and likely to have been present in the early solar PPD. Studies employing isochrone age-determination techniques of star-cluster sources with near-infrared excess indicate that PPDs have statistically short lifetimes of ~2–4 Myr (Haisch et al. 2001; Wyatt 2008; Bell et al. 2013), providing upper limits on the timescale within which macroscopic solids can form from the disc material. Submillimetre observations place even tighter constraints of ~1–2 Myr for the survival of the large-grain population (Lee et al. 2011; Ansdell et al. 2017). The dust and ice grains of μm sizes represent ~1% of the total PPD mass, the rest of which is mainly hydrogen gas. Therefore, first-order PPD models are ideal-gas flows subject to orbital dynamics, with primitive state variables that depend upon distance from the star and height above the disc midplane (Weidenschilling 1977b; Hayashi et al. 1985; Armitage 2010).

Particle–gas interactions determine the movement of solids – with respect to each other and to the gas – at virtually all particle size scales (Turner et al. 2014). Depending upon local circumstances, the flow may be laminar or else turbulence may arise from a host of proposed fluid instabilities in the context of planetesimal formation (Johansen et al. 2014). Although hydrodynamic turbulence is known to create pressure and velocity gradients that sequester particles and facilitate collisions, it can, for the same reasons, serve to diffuse high particle-concentration regions or lead to high-energy particle collisions that are ultimately destructive (Blum & Wurm 2008; Youdin 2010).

The statistical outcomes of dust-aggregate collisions are well mapped from microgravity experiments conducted under vacuum conditions (Dominik & Tielens 1997; Blum et al. 1998, 1999; Blum & Wurm 2000; Güttler et al. 2009). While low-energy collisions result in fractal growth, aggregates accrue momentum with increasing mass, such that further collisions may result in erosion or fragmentation, instead of continued bottom-up growth (Zsom et al. 2010; Güttler et al. 2010; Bukhari Syed et al. 2017).

Drag dissipation from pressure-supported sub-Keplerian gas motions on particles with Keplerian velocities causes particle orbital decay (Whipple 1972). The drift rate of particles possessing typical solid densities peaks for particle diameters 10–100 cm, with associated depletion timescales in the 100s of orbits (Weidenschilling 1977a). Growth by collisional coagulation is thus stifled by the shortsurvival times of the pebbles that disappear due to systematic drift (Brauer et al. 2007; Birnstiel et al. 2012). Mechanisms to bypass the metre-size drift barrier must therefore be rapid, making strict collisional growth inefficacious.

A complementary approach to coagulation models is to follow the dust evolution as if it is a fluid without internal pressure (e.g. Nakagawa et al. 1981), where both the gas and solid phases obey continuum mass and momentum balance equations. Taking into account the collective drag-force coupling between the two phases, Youdin & Goodman (2005) found that such a flow possesses an oscillatory instability, referred to as the streaming instability (SI). The critical control parameters for the onset of instability are the spatially averaged dust-to-gas mass density ratio ϵ = ρ_p ∕ρ_g and the Stokes number τ_f = T_fΩ, defined in this context as the viscous coupling time T_f normalised by the particles’ Keplerian frequency Ω. For a particle with mass m_p and velocity |v|, the viscous coupling time derives from comparing its momentum to the drag force, F_d, that it experiences: $T_{\text{f}}=\frac{m_{\text{p}}\left|v\right|}{F_{\text{d}}}$. The growthof SI is strongest when particles are weakly decoupled from the gas (τ_f > ~0.1) and when ϵ is close to unity. When SI becomes fully developed, particles in regions with high particle density feel a reduced collective drag and therefore the mass density can become further enhanced by isolated particles that catch up to the slower drifting cluster. The amplitude of the resultant density wave could become sufficient for gravitational instability to occur and form planetesimals by direct collapse (Youdin & Johansen 2007; Johansen & Youdin 2007).

This effect has become integral to a number of subsequent theoretical studies on disc evolution and planet formation (e.g. Johansen et al. 2009, 2012; Bai & Stone 2010a, b, c; Miniati 2010; Kowalik et al. 2013; Blum et al. 2014; Yang & Johansen 2014; Carrera et al. 2015; Simon et al. 2017; Schäfer et al. 2017). The paradigm revisits the sequence originally proposed by Safronov (1969) and Goldreich & Ward (1973), whereby dust grains are first concentrated by vertical settling and collisional coagulation, followed by gravitational fragmentation of the resultant dust sub-disc.

A common finding in the aforementioned theoretical works is that the onset rate of the instability is appreciable for ϵ ≳ 1. This is consistent with the consensus from studies of dispersed two-phase flow, where two-way drag coupling is significant for ϵ = 1 and four-way coupling – fluid drag and particle back reaction, plus collisional momentum exchange – is significant when volumetric filling factors ϕ increase (Balachandar & Eaton 2010; Marchioli 2017). The work of Dra̧żkowska & Dullemond (2014) shows that the collisional coagulation process can produce agglomerates large enough to seed SI, with especially favourable conditions beyond the radial distance at which ices can form (Bitsch et al. 2015; Drążkowska & Alibert 2017; Schoonenberg & Ormel 2017).

Recent studies have sought to generalise the SI and understand its model dependencies. Jacquet et al. (2011) demonstrated that a laminar two-fluid shearing sheet is linearly stable to SI when Ω = 0. Lambrechts et al. (2016) numerically evolved the same model equations, that is, incompressible two-fluid drag-coupled mass and momentum balance, with terms relating to Keplerian rotation absent, and found the system to be unstable to non-linear perturbations. Together with the lack of requirement of rotation, the presence of an instability on small scales led the authors to suggest that the effect may occur in a range of circumstances where particles experience drag from a dilute gas, such as in the vertical sedimentation phase of dust sub-disc evolution, the dust in the atmosphere of an embedded giant planet, or in the formation of chondrules. Additional works find dust-drag induced fluid instabilities to occur generically and in various stages of PPD formation and evolution (e.g. Bate & Lorén-Aguilar 2017; Squire & Hopkins 2018a,b; Lin & Youdin 2017).

The SI is a robust effect that arises in a laminar disc model but also persists in the presence of turbulence. Independent of how the viscous evolution of PPDs proceeds, SI is considered a promising way to help overcome the long-standing metre-barrier problem. However, small-scale dust concentration is not yet possible to observe directly with astronomical observations, due to the optically thick nature of discs at short wavelengths, because spontaneous concentration of particles due to fluid instabilities is likely occurring over extremely rapid timescales in very young discs, and because the lengthscales involved are too small to resolve.

The topic of angular-momentum transport in PPDs by anomalous turbulent viscosity has been addressed experimentally with Taylor–Couette devices rotating at Keplerian frequencies (Dubrulle et al. 2005; Hersant et al. 2005). The study of dust-grain drag-induced aerodynamic instabilities, which could eventually lead to turbulence, in sufficiently dilute scaled-down laboratory flow has never been previously attempted and is the unique goal of this work.

We devise a very reduced system, starting with the premise that the instability arises principally from the differential motion between the gas (“carrier”) and solid (“disperse”) phases, and that this is simply analogous to the case where particles, initially homogeneously mixed within a fluid, settle against a pressure gradient under constant gravitational acceleration g and relax to their terminal velocity u_t. A minor difference between the setup in Lambrechts et al. (2016) and the present one is that they assumed a hydrostatic pressure gradient against which particles could settle and unmix, whereas the change in hydrostatic column density is negligible in our experimental volume and we instead drive an upward flow by applying a pressure differential. In the reference frame between the disperse and carrier phases, however, the effect of there being a net velocity difference is the same. As will be shown, the particle-tracking data analysis exhibits particle clustering and modified particle dynamics on small scales, which we interpret as a collective particle-drag induced instability.

Particle suspensions in the Stokes regime (continuum drag assumption) have been studied for centuries (see, e.g. Darcy 1856; Batchelor 1972 and references therein and thereto), but broaching into the Epstein-drag regime (Epstein 1924; free molecular flow or relaxed surface slip conditions) is a novelty of this work. In general, solids that participate in planetesimal formation can be in either regime. Similarly, the use of a rarefied gas as the carrier phase, where the value of ϵ ~ 1, results in volumetric filling factors ϕ several orders of magnitude (at least 3–4) lower than what has ever been studied before in the literature of fluidised beds and suspension chambers (see reviews such as Balachandar & Eaton 2010, or Guazzelli & Hinch 2011).

We describe the apparatus and our experimental techniques in Sect. 2. We then outline the parameters of the experiments in Sect. 3. The results of the particle tracking are shown in Sect. 4. We discuss our findings in Sect. 5 and suggest ideas for follow-up in Sect. 6. We summarise and conclude in Sect. 7.

2. Experimental techniques

2.1. Apparatus

The experimental facility is shown in Fig. 1. In a cylindrical plexiglas chamber, solid particles are transported by an upward stream of dilute gas under low pressure. The top of the chamber is connected to a vacuum pump and the base of the apparatus contains a pressure-reducing filter composed of porous material. The filter creates a pressure drop ΔP = 9.9 × 10⁵ Pa from the outside of the apparatus to the inside. The pressure differential across the whole system drives a continuous flow.

The facility is equipped to supply different gases to the apparatus, but all the experiments described here were conducted with dry air as the carrier phase. We choose spherical steel particles for their high mass density and high electrical conductivity to eliminate the effect of electrostatic forces. During each experiment, the steady-state pressure and particle-seeding density are controlled and we record the particle positions from three viewing angles in order to reconstruct their trajectories.

The particles are pre-loaded onto a particle seeding platform, which is composed of a fine wire mesh, through which air can pass – but particles cannot (see right-hand panel of Fig. 2). Particles that reach the top of the container enter a hollow expansion where the gas slows down (due to mass conservation), preventing particles from travelling further upstream (see left-hand panel of Fig. 2 for a view of the expansion chamber). Because the gas velocity close to the wall is smaller, particles that migrate towards the walls fall down, then they funnel towards the particle-seeding platform, where they re-enter in the flow andcirculate. The result of the particle circulation is that the centre-line of the flow has a persistent upward-directed particle stream that is constantly replenished with particles. The particle stream persists indefinitely, therefore we wait until transients in the pressure disappear and take recordings of the particle motions during steady-state pressure conditions.

The particle positions are tracked in time using three-dimensional Lagrangian particle tracking (LPT), described in Sect. 2.3. The stereoscopic camera arrangement is shown in Fig. 1, and we used LED backlighting to image the solid particles using shadowgraphy.

The measurements presented in this section will demonstrate that, before adding particles, the underlying gas flow is incompressible and is laminar. We made time-resolved pressure (temperature) measurements with Pirani pressure heads (thermocouples) secured at three heights in the measurement volume. We measured the radial profile of the upward gas velocity using particle image velocimetry (PIV), described in Sect. 2.2. Many more details of the apparatus and its tests can be found in Capelo (2018); see also the summary section of Lambrechts et al. (2016). Below we describe our measurement techniques.

Fig. 1. Particle fluidisation and sedimentation chamber, shown towards the right of the image, is supported on a heavy aluminium plate and fixed to an extruded aluminium profile. A gas line leads into the pressure-reducing filter at the bottom and vacuumtubing leading to a vacuum pump extends from the hollow gas expander at the top. Cables extending from the device connect pressure transducers to a data-logging system. On the optical table are three high-speed phantom V10 cameras, positioned to have overlapping field of view inside of the apparatus, and the signal generator used to synchronise the three cameras via external triggering. LED spotlights are mounted near the apparatus for use as backlighting.

2.2. Flow profile

The carrier-phase flow was calibrated using PIV, with smoke aerosols as the seeding material. This technique involves correlating the positions of tracer particles in subsequent image pairs to extract the instantaneous vector flow field. An average of time-series PIV data was used to determine the mean field and its fluctuations. Due to the axisymmetry of the flow in question, we determined the radial flow profile from two-dimensional measurements of a single slice through the centre-line of the apparatus. Particles with short viscous coupling time, T_f, are desirable to approach tracer particles and, due to the high contrast in density between the gas and solid particles, it is necessary to use the smallest possible tracers. We therefore generated μm-sized smoke particles by combustion at atmospheric pressure outside the apparatus and entrained the fluidised aerosols into the bottom of the sedimentation chamber, where the flow was maintained at a steady-state operational pressure of ~10 mbar.

The measurements were made using one high-speed camera placed perpendicularly to a laser sheet¹ from a single cavity Q-switched Nd:YAG laser at 532 nm (IB Laser Chronos 400 mm IC SHG). The laser was operating at a repetition rate of 3.3 kHz, with a pulse duration of 140 ns and pulse energy of 6.4 mJ. A signal generator was used to synchronise the rising edge of the laser pulses with the camera. The image pairs consisted of 20 μs exposures recorded at a temporal separation of 300 μs. We measured for a 6 s duration, with an image pair acquired every 450 μs.

We used the commercially available DaViS 7.3 software² to perform the image pre-processing and processing. The pre-processing consisted of the following steps: (i) mean background image subtraction; (ii) local non-linear subtraction; (iii) sliding-minimum filter over a scale of 64 pixels; (iv) linear 3 × 3 pixel Gaussian smoothing filter; (v) intensity renormalisation filter. The processing involved interrogating pre-processed image pairs with a 3-pass PIV scheme, with the first pass using 1.5 × 1.5 mm windows at 50% overlap and subsequent passes using 0.75 × 0.75 mm windows at 50% overlap.

Line-profiles of the two-dimensional vector fields in the PIV time-series data demonstrate that the gas flow has a cylindrical Pouseille profile, with a centre-line velocity of 1.45 m s⁻¹. Figure 3 zooms into the central 2 cm of the gas velocity profile. The measured mean profile is shown with discrete points, indicating the spatial resolution of the vector field. The error bars represent the fluctuations (standard deviation) on the mean flow. A parabolic fit to the mean field, shown as a solid purple line in the figure, goes to zero at radii of ± 45 mm (not shown), as expected for a laminar flow with small container-scale Knudsen number (Sharipov 2011). Although there is radial variation in the flow velocity, it has a difference of no more than 5% over the central 2 cm, which is about twice the size of the region where we conduct LPT measurements.

Fig. 2. Cut-away views of components used to control the gas pressure and two-phase flow velocity. Left panel: hollow gas expansion unit where particles enter and slow when they reach the top of the sedimentation chamber. Right panel: apparatus base with pressure-reducing filter composed of sintered metal plates of porosity coefficient α = 10−12 m2. The funnel insert holds a fine wire mesh that is used as a particle staging platform; low-pressure air streaming through the mesh fluidises the particles.

2.3. Inertial particle motion

In the past decades, traditional Eulerian flow-measurement techniques have been complemented by new methods to study individual flow elements in their Lagrangian frame. Using multiple synchronised high-speed cameras with overlapping fields of view, fluid motions of tracer particles have been used to further understand the energy transfer across scales in turbulent flows (La Porta et al. 2001; Ouellette et al. 2006; Xu et al. 2006, 2007; Xu 2008), and the behaviour of inertial particles in turbulent fluids has been studied as well (Bourgoin 2006; Xu & Bodenschatz 2008; Gibert et al. 2010; Saw et al. 2012, 2014). Although we have devised a laminar gas flow, one expects that the presence of the disperse particle phase will result in complex turbulent-like motions, if an instability develops. We therefore employ the same particle-tracking technique used in previous studies of turbulent flow, with customisation for the present setup.

We acquired simultaneous shadowgraph images with three high-speed Phantom 10 cameras at a frame rate of 2000 Hz and frame size of 768 × 1024 pixels. The spatial scale of the recordings is ~12 μm pixel⁻¹. After the image acquisition has taken place, the particle trajectory reconstruction consisted of the following steps in the order stated: (i) particle image finding; (ii) stereoscopic position reconstruction; (iii) tracking in time.

In step (i), we performed a background image subtraction and then determined the sensor position of each particle in all movie frames by performing a Gaussian fit to pixel intensity maxima. In step (ii), we applied a least-square minimisation to the triangulation errors on a mapping between sensor position and real-world coordinates.The laboratory-frame coordinates are known via acalibration procedure in which the three cameras shown in Fig. 1 were focussed on a grid-pattern calibration mask inserted in the centre of the tube. The mask was translated in the laboratory frame, with a calibration image being obtained at 2 mm separations.

It is common to construct the sensor-to-real-world mapping using a “pinhole” model, which assumes that all points inthe volume can be traced back to the imaging plane via a mutual intersection point. However, we must also account for aberrations in projected particle position due to the cylindrical walls of the experimental vessel, and we therefore implemented a generalised pinhole model that assumes that the lines of sight can pass through twomutually perpendicular slits. We also correct for smaller scale distortions due to variations in the thickness of the experimental vessel wall by fitting a distortion profile along the horizontal axis. We verified the two-slit camera model goodness-of-fit by performing the reconstruction on the calibration images. We subtract the real-world coordinates from the reconstructed positions and find a standard deviation in offset of particle position of ~2 μm, which we take to be the error on particle position.

In step (iii), the remaining task is to link the positions of particles in subsequent frames, thus producing trajectories in time. In each frame, likely particle positions are estimated based upon previous particle positions and a maximum velocity in the previous frame. We note that the particles’ displacement, due to their velocity is hundreds of μm, which is much greater than the noise due to the position reconstruction, stated above. At each frame in which a particle appears, we use the minimum variance linear extrapolation from the previous trajectory positions to determine which trajectory the particle belongs to.

From the three-dimensional particle trajectories, the instantaneous velocity and acceleration at all times is calculated by applying a finite difference scheme, smoothed with a Gaussian filter (Mordant et al. 2004). The filter lengths for the velocity and acceleration were 10 and 40 frames, respectively. We also calculate the moments and distributions of the velocity and acceleration, in addition to other explicit quantities that will be reported in the following sections.

Fig. 3. Flow profile calibration measurement using PIV, limited to the region of possible overlap with the LPT measurements, which are conducted within a ~1 cm3 region centred within the vessel with an accuracy of ~0.5 cm. Top panel: measured mean velocity (blue dots), and parabolic fit (purple line), limited to the central 2 cm of the flow vessel. Errorbars represent the standard deviation of the time-series PIV data. Bottom panel: calculation of the relative difference invelocity from the maximum centre-line velocity. The shaded region represents the uncertainty range due to error bars in the top panel.

Fig. 4. The neighbourhood of the origin of the LPT measurements, where the black, red, and blue lines correspond to the lines of sight and centre point for cameras 0, 1, and 2, respectively. The thick purple lines delineate the measurement volume, Vmeas, defined by the intersecting of lines of sight from all three camera sensors.

3. Properties of the data

3.1. Parameters and their definitions

We conducted experiments at fixed gas pressures in the range 200–800 Pa (2–8 mbar). Each experiment consists of a 6-s LPT recording, processed in the manner previously stated. A single data set includes several repetitions of the experiment ($\mathcal{O}(10)$) at a given gas pressure. To perform a statistical analysis, we combine the three-dimensional trajectories from all experiments of the data set. Table 1 summarises the parameters of the experiments, including the pressure at which the facility was operating and the total number of trajectories, N_traj, comprising the data set.

Flows with vastly different length-scales, L, will exhibit statistically similar velocity fields, if their Reynolds numbers, Re = uL∕η, are the same. In the numerator of Re is the flow’s characteristic velocity u and in the denominator is the kinematic viscosity: the dynamic viscosity μ divided by the fluid density η = μ∕ρ_g. When the Reynolds number is high, the flow will become unstable and even turbulent. Conversely, in low-Re flows the viscous damping is significant and the flows remain laminar. For gases with a normal temperature range, the dynamic viscosity is roughly constant, especially when the pressure is low. Therefore, low gas density ρ_g means high η, consequently low Re, which implies a laminar flow. The dynamic viscosity of dry air at room-temperature in our vacuum system is μ_air = 1.8 × 10⁻⁵ kg m⁻¹ s⁻¹, which gives a kinematic viscosity of η_air = 2.0–7.9 × 10⁻³ m² s⁻¹ for the pressure range 200–800 Pa. Therefore, for the pipe flow on the container scale, $R{e}_{\text{C}}=D_{\text{pipe}}\overline{u}/\eta \approx$ 15–45, and for flow at the scale of an individual particle with diameter d_p, Re_p = d_pv_rms,z∕η ≪ 1, where D_pipe is the pipe diameter, $\overline{u}=0.95$ m s⁻¹ is the average velocity in the pipe, d_p is the particle diameter, and v_rms,z is the root-mean-square (rms) velocity of the particles in the z-direction. These values are listed in Table 1.

The mean number of particles in the measurement volume over all experiments in a data set is n, where the error reported is the standard deviation in n across individual iterations of the experiment. From n, we derive a mean volumetric filling factor, $$\varphi =\frac{\pi }{6}n{d}_{\text{p}}^3/V_{\text{meas}},$$(1)

where V_meas is the volume of the measurement region, and spatially averaged mass loading $$ϵ={\rho }_{\text{p}}/{\rho }_{\text{g}}=\varphi {\rho }_{\text{m},\text{p}}/{\rho }_{\text{g}},$$(2)

where the material density of the stainless steel particles used in the experiments is ρ_m,p = 8050 kg m⁻³. We also calculate the momentum diffusion time across the mean inter-particle distance, t_d = n^−2∕3∕η, to demonstratethat the experiments are in a tight particle-coupling regime, as indicated by this timescale being shorter than the viscous relaxation time (see also Eq. (4) of Lambrechts et al. 2016).

All of the experiments listed in Table 1 were conducted with the same distribution of particle diameters: 15–65 μm initially seeded into the apparatus. On the other hand, because the experiments are conducted in a flow regime where the velocity of particles depends upon ρ_g, for each experiment under a given gas pressure, the flow naturally selects only particles with diameters in a subrange of this size distribution by settling and levitation.

We determined the particle diameters by using the fact that at steady state, the measured relative velocity between the particles and the gas stream should equal the terminal velocity of a particle of size d_p. It is common to assume a piece-wise drag-law transition for Re_p < 1, $$F_{\text{D}}=\frac{\pi }{6}{d}_{\text{p}}^3{\rho }_{\text{m},\text{p}}g=\{\begin{array}{ll}3\pi d_{\text{p}}{\mu }_{\text{air}}\delta u\hfill & \text{Stokes},\hfill \\ \frac{\pi }{3}{\rho }_{\text{g}}{d}_{\text{p}}^2{v}_{\text{th}}\delta u\hfill & \text{Epstein},\hfill \end{array}$$(3)

where δu = u_g − u_pz is the difference between the z-component particle velocity, u_pz, and the gas velocity, u_g = 1.45 m s⁻¹, $v_{\text{th}}=\sqrt{8R_{\text{g}}T/\pi M}$ is the thermal speed of molecules of molecular weight M, and R_g is the universal gas constant.

The resistance that a sphere experiences is thus broadly classified into whether continuum viscous forcing, or rather impingement from individual molecules, is responsible for the drag. Between these two extremes, the boundary conditions of the flow around the particle will be characterised by increasing slip on particle boundaries. The appropriate drag regime (Whipple 1972; Weidenschilling 1977a) is given in general by the Knudsen number: $$Kn\equiv \lambda /d_{\text{p}}.$$(4)

Further sub-divisions in flow regime use the following as approximate guidelines (see, e.g. Allen & Raabe 1985; Karniadakis et al. 2005): $$\begin{array}{lll}\hfill & \hfill & Kn<0.01\text{\hspace{0.17em}continuum\hspace{0.17em}}(\text{Stokes}),\hfill \\ \hfill & \hfill & 0.01<Kn<0.1\text{\hspace{0.17em}continuum\hspace{0.17em}}(\text{with\hspace{0.17em}free\hspace{0.17em}slip\hspace{0.17em}correction}),\hfill \\ \hfill & \hfill & 0.1<Kn<10\text{\hspace{0.17em}transition},\hfill \\ \hfill & \hfill & 10<Kn\text{\hspace{0.17em}free\hspace{0.17em}molecular\hspace{0.17em}}(\text{Epstein}).\hfill \end{array}$$

Figure 5 shows λ as a function of pressure, with the Kn-dependent regimes color-coded for reference. Particles with sizes used in our experiments are represented by vertical green bars spanning 15 < d_p < 65 μm, with orange sub-segments indicating d_p selected by the flow.

Equating the expressions for the Stokes and Epstein drag forces gives an algebraic dividing point for particle diameter of (9∕2)λ, which is shown as a dashed line and divides our data sets in Fig. 5.

While our data is representative of the two different drag regimes delineated by Eq. (3), Fig. 5 also demonstrates that even the data collected at higher pressure can be considered transitional flow and that a purely continuum assumption for the gas flow would require significantly higher pressure or much larger particles. In observance of this issue, we further constrain d_p, as is appropriate for transitional flow, using the empirically derived Kn-dependent Cunningham correction factor (Cunningham 1910) to Stokes drag: $$F_D=F_{D,\text{Stokes}}/C_{Kn}$$(5)

with $C_{Kn}=1.0+2Kn\left(\alpha +\beta \exp (-\gamma /2Kn)\right)$. We adopt the values from Allen & Raabe (1985): α = 1.142, β = 0.558 and γ = 0.99. The values of d_p obtained using this correction are shown in Table 1 and are those used to calculate the other listed parameters³. For comparison, we list d_p,St−Ep obtained using Eq. (3). In both cases, the errors are calculated by propagating the variance on the mean settling velocity in the z-direction, denoted ${\sigma }_{\langle u_{\text{z}}\rangle }$ and also listed in Table 1, and demonstrate narrow dispersion in particle size. The quantity ${\sigma }_{\langle u_{\text{z}}\rangle }$ is distinct from and much smaller than the fluctuation velocities reported in Table 2 below.

We verify that the particles have had time to reach their terminal velocity before reaching the measurement window, which is h_meas = 5 × D_pipe ~ 50 cm downstream of the bottom mesh where the particles settle when there is no flow, by solving the dynamic equation for particle velocity, $$\frac{\text{d}{u}_{\text{pz}}}{\text{d}t}=\frac{u_{\text{g}}}{T_{\text{f}}}-\frac{u_{\text{pz}}}{T_{\text{f}}}+g,$$

which has the solution $$u_{\text{pz}}=(u_{\text{t}}-u_{\text{g}})\left(e^{\frac{-t}{T_{\text{f}}}}-1\right)=(g{T}_{\text{f}}-u_{\text{g}})\left(e^{\frac{-t}{T_{\text{f}}}}-1\right),$$(6)

where u_t = gT_f is the terminal velocity. At long times, such as when t ≫ T_f, u_pz reaches its steady state value u_p,ss = u_g − u_t. Equation (6) is used to obtain the relation between the height $z_{\text{p}}={\displaystyle {\int }_0^t{u}_{\text{pz}}}\text{d}t$ above the bottom mesh and the normalised particle velocity u_pz∕u_p,ss, and between z_p and the normalised time t∕T_f. The results are plotted in Fig. 6, which shows that particles closely approach steady-state velocity in no more than 20 cm above the bottom mesh, which corresponds to a time approximately 4–6 T_f and that is indeed justified for the time-dependent term to be neglected according to Eq. (6).

Fig. 5. Flow conditions at the scale of particle size in the experiments. Black curve: the mean free path of gas molecules as a function of pressure λ(P). Coloured swaths denote Kn-dependent drag regimes. Grey dashed curve: Stokes–Epstein transition. The four vertical green bars represent each of the four data sets; the position of the bars denotes the pressure, their height denotes the spread in particle size for all experiments, and their width denotes the offset error in the pressure measurement. The orange subsegments indicate the particle sizes that have the steady-state velocity equaling the measured particle velocity in the experiment: up,ss = ⟨upz⟩. All the data sets can be considered transition flow, with the two at higher pressure further towards continuum and those at lower pressure approaching free molecular flow.

Fig. 6. Left panel: height as a function of the ratio of particle velocity to the steady-state velocity upz ∕up,ss. Right panel: height as a function of time in units of the viscous relaxation timescale Tf. Particlesreach a steady state well before reaching the measurement window at ~ 50 cm.

3.2. Summary statistics

The particle velocities, in each of the data sets listed in Table 1, are by far dominant in the vertical-component, with 0.3 m s⁻¹ ≲ ⟨u_pz⟩ ≲ 0.5 m s⁻¹. By comparison, the horizontal components are 0.001 ≲ ⟨u_px,y⟩ ≲ 0.005 m s⁻¹. We note that in the laboratory frame, the particles travel vertically up during the experiments, against gravity. However, because the mean gas velocity is much larger, the particles are actually falling downward relative to the gas, with a relative velocity δu in the range of ~0.9–1.2 m s⁻¹.

We also find comparably larger velocity spread in the vertical than in the horizontal components, as quantified by the rms velocity for each component i, defined as $u_{\text{rms},i}\equiv {\langle u_i^2\rangle }^{1/2}$, which we also refer to as the velocity dispersion or fluctuation. Table 2 shows u_rms,i at the gas pressures we used in the experiments (data set).

We measured gas-velocity fluctuations at the flow centre-line of 0.0098 m s⁻¹ when the pressure is at 10 mbar, corresponding to Re_C = 56, which is larger than the Reynolds number at any of the lower pressures at which our experiments were conducted (see Fig. 3 and Table 1). Such fluctuations in the gas phase velocity may contribute to the dispersion in vertical particle velocity. Its effect is expected to decrease at lower gas pressures as the Reynolds number is smaller, which is consistent with the observed smaller particle velocity fluctuations at lower gas pressures. Given that the fluctuations in the gas are smaller thanthose in the particle phase, the trend of increasing particle velocity fluctuations with increasing gas pressure is most likely related to the different dynamics associated with the drag regime to which the particles belong, which will become evident due to the results in Sect. 4.4. We note that the dispersion in particle size can also contribute to the rms velocities, but, the fluctuations are much too large to be accounted for by this effect. While we have found that there is a small variance on the mean settling velocity, this does not effect the rms values in Table 2 because the global mean of the data set is not subtracted as an ensemble. In the appendix, we plot the probability density functions (PDFs) of instantaneous particle velocity and acceleration, which show tails much wider than Gaussian.

The possibility of a fluid instability complicates arguments based on simple terminal velocity assumptions of the particles, since one would expect a complex velocity field. However, the low level of gas velocity fluctuations compared to the mean gas flow (i.e. see shaded region in bottom panel of Fig. 3) and compared to the values of u_rms,z, together with the fact that the effect of the mean shear in the gas flow is negligible, indicates that the gas flow itself is not the direct reason to cause inhomogeneities in the particle phase, if there are any.

In the following, we therefore present our analysis of the particle motions with the focus on particle–pair relative velocity statistics. We will use u and u_i to denote particle velocity and its components and drop the subscript “_p” for simplicity, since we henceforth only discuss particle velocity and its statistics.

4. Results

In this section, we present the results of the experiments conducted at the conditions given in Table 1. We present an overview of the magnitudes of the relative motion between particles in Sect. 4.1, which shows that the relative motion between particles is not determined by their terminal velocities caused by their size differences, nor by the gas velocity profile. In Sect. 4.2, we consider the components of the relative velocity, which demonstrates the occurrence of particle aggregation and clustering on small scales.

4.1. Magnitude of the relative motion between particles

If particles of the same size in a uniform gas stream are moving at their terminal velocities relative to the gas as if they are isolated, then the relative velocity among particles would be zero. Even if we consider the dispersion of particle size (and hence their terminal velocity), the statistics of the relative velocity between particles would not depend on the distance between particles, provided that particles with different sizes are uniformly distributed in space. In this section, we first study the variation of the magnitude of the relative velocity between particles with the distance between particles.

At any instance t, for two particles located at positions x and x +l, with velocities u(x, t) and u(x +l, t), respectively, their relative velocity is δu = u(x +l) −u(x). In order to see how the statistics of δu depend on l, we consider the following quantity: $$D_{jk}(l)\equiv \langle [u_j(x+l,t)-u_j(x,t)][u_k(x+l,t)-u_k(x,t)]\rangle ,$$(7)

where “⟨⋅⟩” means ensemble average, that is the average over all particle pairs that are separated by l, for all x and t. If the particles are fluid tracers, thenthe quantity D_jk defined by Eq. (7) is the so-called Eulerian second-order velocity structure function of the flow field. For heavy particles, like those we study here, it gives us a direct measure of the magnitude of the relative velocities between particles. As we are more interested in how two particles approach or separate from each other, we consider only the second moment of the relative velocity component in the direction along their separation vector, $$D_{ll}(l)\equiv \langle {\left[(\delta u)\cdot \frac{l}{l}\right]}^2\rangle .$$(8)

In Fig. 7, D_ll(l) for all data sets are shown as open symbols. The dependence on scale l is clearly visible. A solid-line curve is plotted to show the difference in gas velocity between a point on the vessel’s centre-line and a point at radial distance l, due to the Poiseuille flow profile. If particles simply adjust to the spatially varying gas velocity, this would be the maximum possible velocity differences they could have. The curve is well below all measured values of D_ll, which means that the variation of D_ll with l cannot be attributed to the variation in the gas flow profile. This is another evidence that the mean gas flow does not contribute much to the particle velocity variations.

For all four data sets, D_ll are rather small at separations less than ~1 mm, indicating that particle velocities may be correlated at very small separations. In all cases, D_ll increases with l, showing that statistically the relative velocity increases with separation. The result from data set kn1_2mb [DS1] is noisy, which is due to the low number density in the experiment. As shown in Table 1, there is on average one particle in the measurement volume at any given time, which results in poor statistics of particle relative velocity. Similarly, data set kng1_5mb also has relatively low number density and the statistics suffers from convergence problem as well, at least at small separations.

The second moment of the relative velocity shows how intense the relative motion is. While some hints of collective motion can be implied from that statistics, it is more informative to study directly the average relative velocity when checking for particle aggregation, as we will do next. Since kn1_2mb [DS1] is not a large enough sample for relative velocity analysis, we do not consider this data set further. Of the other three data sets which exhibit the potential for collective behaviour, we select kn1_3mb [DS2] as representative of the Kn ~ 1 regime and kng1_8mb [DS3] as representative of the Kn < 1 regime. We only show kng1_5mb [DS4] when the results are not obscured by the wide range of particle number density in the experimentscomprising the data set; note from Table 1 that the fluctuation in the number density n is greater than the mean, hence although there are some high number density experiments, the mass loading is on average low.

Fig. 7. Second-order moment of the longitudinal relative velocity component, Dll(l), versus l for all data sets. Purple solid curve corresponds to the magnitude of the velocity difference of the gas flow profile as a function of radial distance from the flow centre-line. The contribution to Dll from the variation in shear velocity in the gas phase is thus negligible.

4.2. Particle aggregation revealed by the relative velocity

Now we explore the details of the relative motion between particles. In the last section, we presented statistics only as a function of distance, without concerning the directional variations. The particle motion is actually rather axisymmetric, as indicated by the fluctuation velocities shown in Table 2. Therefore, we present here the average of relative motion in cylindrical coordinates, where the z-axis is the same as in the laboratory frame and the r-coordinate is in the horizontal plane. The origin, however, is not fixed, but set on any given particle. We study the relative velocity $$\delta u_i(\delta r,\delta z)\equiv u_i(r+\delta r,z+\delta z)-u_i(r,z),$$(9)

which is a function of δr and δz but does not depend on r and z for homogeneous flow fields, which is an accepted approximation for our cases as we perform the measurements in the central region of the apparatus, where the motion of particles is hardly influenced by the wall. Moreover, because of the symmetry of our setup, the relative velocity does not depend on the direction of δr. Thus, the relative velocity varies only with the magnitudes of δr and δz.

For the two data sets of Kn ~ 1 and Kn < 1, the mean of δu_i (δr < 1 mm, δz), as a function of δz is shown in Fig. 8, in which the statisticsof δu_i is limited to the cases with two particles radially very close, δr < 1mm, that is the two particles are nearly one on top of each other, but with a vertical distance δz.

The means of δu_x and δu_y are nearly the same, thus supporting our axisymmetry arguments. Both ⟨δu_x⟩ and ⟨δu_y⟩ are negligibly small, which is consistent with the very small fluctuation velocities in the horizontal plane (see Table 2).

The relative velocity in the vertical direction, however, is interesting. Here a negative ⟨δu_z ⟩ means that on average the particle that is below is moving faster and hence would “catch up” to the one that is above, which implies aggregation. Particles clearly aggregate towards one another at an increasing rate for diminishing separation in δz. This behaviour occurs for both the Kn ~ 1 and Kn < 1 cases. The separation at which this effect is most pronounced is for δz ≲ 4 mm. In neither data set does the relative velocity return to zero for increasing separation, discussed below.

In Fig. 9, we plot ⟨δu_z(δz)⟩ for other ranges of δr, such that the second particle lies in a cylindrical shell with different radius around the first particle. The statistics of the horizontal velocity components ⟨δu_x⟩ and ⟨δu_y⟩ is not shown, since they are nearly vanishing in all cases. Results in Fig. 9 show very weak dependence on δr, which implies that the aggregation occurs in thevertical direction, at least at the scale of the measurement volume. If this is due to the instabilities in the motion of the particle phase, the wavelength of the unstable modes should be at least of the size of the measurement volume. Theoretical calculations presented in Lambrechts et al. (2016) predict that the smallest lengthscale, λ_knee, of the particle flow should be around λ_knee = 0.08L_f. Using valuesof T_f from Table 1, λ_knee ≈ 1 cm for the two data sets discussed here, which is consistent with our observation.

Figures 8 and 9 show that δu_z does not vanish even at the largest measured separation δz ≈ 10 mm, consistent with our estimate above that the smallest lengthscale of the unstable modes in particle motion is ~1 cm, such that our entire measurement volume is always in one of such “contraction regions” when we observe particles. Please see also details of the higher order behaviour of the same quantity in Appendix B. We acknowledge that the apparent limiting behaviour may be consistent with a developing phase lag, of the type proposed in Lin & Youdin (2017). Further identification of the mechanism requires simultaneous measurements of the local gas speed or pressure, together with inertial particles, which is a very challenging task and beyond the capability of our apparatus.

Fig. 8. Mean relative velocity components between particle pairs in units of m s−1, as a functionof vertical separation. The Kn ~ 1 data is shown in the top panel and the Kn < 1 data in the bottom panel. Particles aggregate for separations of 4 mm or less in both data sets.

Fig. 9. Mean vertical relative velocity between particle pairs in units of m s−1, as a functionof vertical separation. The Kn ~ 1 data is shown in the top panel and the Kn < 1 data in the bottom panel. The different curves refer to increasingly large values in δr. The aggregation dynamics is the same at all radial (horizontal) separations, suggesting the formation of a particle layer.

4.3. Clustering of particles seen from concentration variation

The statistics of the relative velocity between particles indicates that particles tend to aggregate, which implies that the spatial distribution of particles cannot be uniform, nor stationary. Now we explore how the local concentration of particles varies with time. To that end, we divide the measurement volume into horizontal slabs, each of height 1 mm, and count the number of particles in each slab, which is then a function of time: we measure n′ (t, z), where n′ is the number density in the 1 mm slab and z is in steps of 1 mm. In general, the spatial and temporal variation of n′ (t, z) can be studied by the two-position, two-time correlation function $$C_{nn}(\tau ,\delta z)\equiv \frac{\langle n^{\prime }(t,z)n^{\prime }(t+\tau ,z+\delta z)\rangle }{{\sigma }_n^2(z){\sigma }_n^2(z+\delta z)},$$(10)

where ${\sigma }_n^2(z)$ and ${\sigma }_n^2(z+\delta z)$ are the variances of n′ at position z and z + δz, respectively.

Because the mean particle velocity in the laboratory frame ⟨u_pz⟩, which is roughly 0.4 m s⁻¹ for all cases, is much larger compared to the fluctuation velocities, the change in particle relative position is very small when particles fly pass the measurement volume. Hence the particle configuration might be regarded as “frozen” in the measurement and the number concentration n′ measured atfixed z but different time could be used as a surrogate of the instantaneous number concentration in the vertical direction. This is a standard approach in turbulence research when performing fixed-probe measurement in a wind tunnel, where it is called “Taylor’s frozen turbulence hypothesis” and it is valid as long as the fluctuation velocity is much smaller than the mean (Tennekes & Lumley 1972; Pope 2000). Under this approximation, the spatial variation of particle concentration can be probed by the single-location, two-time autocorrelation function of n′ at times t and t + τ $$C_{nn}(\tau ,0)\equiv \frac{\langle n^{\prime }(t,z)n^{\prime }(t+\tau ,z)\rangle }{{\sigma }_n^2(z)}.$$(11)

While C_nn(τ, 0) can provide information on the average size of the local concentration variation, the more revealing statistics is actually C_nn conditioned on the local values of n′(z). We perform two conditional statistics, one for the local number density being greater than 1.5 times the mean number density, denoted $C_{nn}(\tau ,0)|n^{\prime }(t,z)>1.5\overline{n}(t;z)$, and the other for it being less than 0.5 times the mean number density, denoted $C_{nn}(\tau ,0)|n^{\prime }(t,z)<0.5\overline{n}(t;z)$. We refer to the former conditioning as the dense case and the latter as the dilute case.

In Fig. 10, we plot the conditional number density concentration correlation normalised by the unconditional ones, that is $C_{nn}(\tau ,0)|n^{\prime }(t,z)>1.5\overline{n}(t,z)/C_{nn}(\tau ,0)$ and $C_{nn}(\tau ,0)|n^{\prime }(t,z)<0.5\overline{n}(t,z)/C_{nn}(\tau ,0)$, for the Kn ~ 1 data. The dense-case ratio, $C_{nn}(\tau ,0)|n^{\prime }(t,z)>1.5\overline{n}(t,z)/C_{nn}(\tau ,0)$, is always less than one for all z. The converse is true for the dilute-case ratio, $C_{nn}(\tau ,0)|n^{\prime }(t,z)<0.5\overline{n}(t,z)/C_{nn}(\tau ,0)$. To gain intuition for this result, we consider a limiting case: if the particles were arranged in a perfectly homogeneous fashion with the inter-particle separation much smaller than the volume under consideration, then the correlation would remain close to one. Consequently, there would be no distinction between the dense and dilute cases and so the ratio represented by Fig. 10 would also remain constant at one. On the other hand, if the particles are non-uniformly distributed, when we start thestatistics from a locally dense region, then the correlation would decrease as we move away from the dense region. Hence, the normalised conditional correlation is smaller than one, and similarly, the normalised correlation conditioned on the dilute case is larger than one. Furthermore, the extent of the deviation from one is an indication of the size of the inhomogeneity in number density. One expects for the local number density to de-correlate faster, the more localised that inhomogeneities of the mixture are. At larger distances, the effect of local concentration should vanish and the ratio of either the dense or the dilute case to the unconditional value of C_nn should eventually return to 1, which is indeed what we observe in Fig. 10.

In Fig. 11, we plot the conditional number density concentration correlation normalised by the unconditional ones for the Kn < 1 case, which are basically the same as for the Kn ~ 1 case, except that the magnitude of the ratios is slightly smaller, which could be due to the relatively dense overall seeding.

From Figs. 10 and 11 we see that the inhomogeneities in particle number density de-correlate in approximately 0.03–0.05 s, which can be converted to approximately 1–2 cm when we use the mean particle velocity in the laboratory frame ⟨u_pz⟩ ≈ 0.4 m s⁻¹. This lengthscale is actually larger than our measurement volume size and thus cannot be directly observed. The frozen hypothesis, on the otherhand, provides a tool for us to probe this spatial inhomogeneity beyond the size of the measurement volume. We also note that this lengthscale is consistent with the wavelength of the most unstable modes that we estimated in Sect. 4.2.

In Fig. 12, we plot the cross-correlation of particle number densities at two different locations as defined by Eq. (10). For both the kn1_3mb [DS2] (top) and kng1_8mb [DS3] (bottom) data, C_nn (τ, δz) decreases in magnitude and broadens as the distance between the two locations increases, which is consistent with the behaviour of an advection-diffusion process. The larger value of C_nn(τ, δz) at long time-lag in the bottom panel is simply due to the higher average particle number density of the kng1_8mb [DS3] run. The time at which the cross-correlation peaks corresponds to the time the particles take to move through the distance dz at the mean velocity. The broadening of the cross-correlation curves is due to the velocity fluctuation of the particles. For the kng1_8mb [DS3] (bottom) case, the rms velocity is much larger than the kn1_3mb [DS2] (top) case, as shown in Table 2, which explains the faster de-correlation of the particle number density observed in the bottom panel.

This broadening of the cross-correlation curve implies that the clustering of particles that we observed in the measurement volume is a dynamic effect rather than some initial concentration variations introduced at the bottom of the apparatus; our measurement volume is ~50 cm above the bottom where the particles are entrained in the air stream, which converts to ~1.25 s of travel time at a mean particle velocity of ~0.4 m s⁻¹. Consequently, any initial inhomogeneities in particle concentration would have been smeared out, judging from the de-correlation time shown in Fig. 12.

Fig. 10. Conditional autocorrelation function of particle number density, normalised by the unconditional correlation, for the Kn ~ 1 case. Top panel: conditional on the dense case; bottom panel: conditional on the dilute case. Heterogeneity in the particle density is evident from the ratios deviating from unity. The timescale for the ratio to return to unity corresponds to a characteristic lengthscale of the clustering.

Fig. 11. Same as Fig. 10, for the Kn < 1 case, showing generally similar clustering behaviour.

Fig. 12. Temporal variation of the cross-correlation of particle number densities at different heights in the measurement volume. Top panel: the case of kn1_3mb [DS2]; bottom panel: kng1_8mb[DS3]. Shifting of the peaks with time and broadening of the shape indicate that the concentration evolves like an advection-diffusion process.

4.4. Evidence for collective continuum-drag reduction

Having demonstrated particle aerodynamic focussing and spatial clustering, let us now consider how the vertical particle velocities depend upon the local particle density. To do so, we measure the PDF of the vertical component of the instantaneous particle velocities u_z, conditioned upon the number of particles in the 1 mm high slab which the particle under investigation belongs to. The resultant PDFs, P(u_z|N_p(z)), are shown in Fig. 13, where the two panels are for the two data sets with gas pressures in ascending order.

The sets of PDFs corresponding to the two Kn < 1 data sets (middle and bottom panels) show that, as N_p increases, the peak of the distribution both narrows and shifts to lower positive velocity. For the Kn ~ 1 data, there are some differences in the PDFs depending upon N_p, but with no clear trend; perhaps there is some broadening of the distribution with increasing number density but no obvious shift in the mean.

These observations are summarised and clarified in Fig. 14, which shows the mean of all instantaneous vertical particle velocities as a function of the number of particles found in a 1 mm slab. The two curves representing data with Kn < 1 decrease as N_p increases, in contrast to the Kn ~ 1 data set.

While the net velocity of the particles is always positive, a decreasing positive value represents a tendency for the particles to resist being pushed upwards by drag forces. The dependence on number density for the two Kn < 1 data sets strongly suggests collective particle-drag reduction. That we do not observe this behaviour for the lower gas-density data set could be attributable to the considerable slip one expects on the surface of the particles. Indeed, how the mean u_p (z) is always slower for this data set than for the other two is consistent with the fact that, for the same gas velocity, u_t is higher for particles in the Epstein drag regime and so one expects greater δu between thecarrier and disperse phases. Similarly, the momentum diffusivity, expressed by the dynamic viscosity, μ =η∕ρ_g, is higher in a low density medium and so perhaps the local gas accelerations caused by the particles decay more readily than when the momentum is allowed to diffuse more slowly.

Figures 13 and 14, which show different outcomes with regard to collective drag reduction for each of the two representative data sets, provide more insight into the slightly dissimilar behaviour observed between the two data sets in Sects. 4.2 and 4.3. This will be discussed in Sect. 5.4.

Fig. 13. PDFs of z-direction velocity, conditional upon the number of particles in a 1-mm slab, for kn1_3mb [DS2] (top panel) and kng1_8mb [DS3] (bottom panel). Curves shift towards lower velocities with increasing number density in the lower panel.

Fig. 14. Mean settling velocity in m s−1 for kn1_3mb [DS2] in blue, kng1_8mb[DS3] in red, and kng1_5mb[DS4] in green, conditional upon the number of particles in a 1 mm cylindrical slab. Drag dissipation is reduced for high density regions in the Kn < 1 data, yet notin the Kn ~ 1 data.

5. Discussion

5.1. Inter-particle interactions and the consequences of low ϕ

The experimental apparatus we use for this study bears some resemblance to a fluidised bed, which is a type of system known to show various signs of collective particle phenomena. An important distinction between this and previous studies is the drastic contrast in material density of the two phases, which is within ${10}^7<\frac{{\rho }_{\text{p},\text{m}}}{{\rho }_{gas}}<{10}^{10}$, depending on the gas pressure of our experiment. Therefore, we can achieve relatively high average solid-to-gas ratio, ϵ ≈ 1, expressed by Eq. (2), even with extremely low average filling factors, 10⁻⁷ < ϕ < 10⁻⁶, expressed by Eq. (1). The parameter range of the system in question can be compared directly to the situation in a PPD. Assuming a MMSN disc, at 1 AU, the material-density contrast between rock, with ρ_rock,m ~ 3 g cm⁻³, and gas is $\frac{{\rho }_{\text{rock},\text{m}}}{{\rho }_{gas}}=2\times {10}^9$. For a dust-to-gas ratio of unity in the midplane, we obtain a volume filling factor ϕ ≈ 4.3 × 10⁻¹⁰. Such extremely diffuse mass distribution is in some sense the essence of why rapid particle concentration mechanisms such as SI are required to facilitate self-gravity.

It might be imprudent to try to study this extremely diffuse limit, either of the carrier or disperse phase, in a laboratory flow, since the total volume under consideration can never approach interplanetary lengthscales and so continuum fluid assumptions would not be met. Still, this study mimics the extremely diffuse regime, due to the relatively large values of Kn on the particle scale and particularly in the sense that the system meets the criterion for the solid phase to be considered a pressureless fluid. In a PPD, the solid component is pressureless because the particles do not – in their laminar equilibrium drift configuration – directly interact, either collisionally or viscously. The dusty fluid therefore has no effective “temperature” in the kinetic sense. The laboratory flow we produced is the same, which is quantified by the criterion T_f > t_d, demonstrated to be true in Table 1. This criterion is demonstrated to be true in the MMSN by Lambrechts et al. (2016).

Collisions were not observed in our experimental system. We estimate the collision frequency in the measurement volume between any two particles, per unit time, Z, as is done in the kinetic theory of gases: $$Z=\frac{1}{2}{V}_{\text{rel}}\pi d_{\text{p}}^2{n}^2/V_{\text{meas}}.$$(12)

Using the values of d_p and n reported in Table 1. Figure 7 shows that the square root of D_ll for particles with separations l ≲ 1 mm is approximately 0.05 ms⁻¹, consistent with the findings in Figs. 8 and 9, that the relative velocity at separations δr < 1 mm and δz < 1 mm approaches 0.05 ms⁻¹, which we adopt as V_rel. We obtain Z ≈ 10⁻² collisions s⁻¹ in the measurement volume, and therefore we could observe no more than one collision in each of our experimental data sets, which are each, cumulatively, ~30 s. We note that the collision frequency calculated this way is an over-estimate, because the relative velocity is dominated by the vertical component, and so a direct collision will also depend upon two particles having a favourable vertical alignment.

5.2. On the negligibility of electrostatic forces

The relative velocity statistics in Sect. 4.2 bear direct witness to particles in the act of spontaneously aggregating. We recall that the particles are always on average approaching one another. However, notably, they do so only in the direction coincident with the mean direction of drag force. We stated in Sect. 2.1 that we ignore electrostatic forces between our particles. We reinforce this point to explain why the clustering phenomemon we observe is fluid-dynamical in nature and unrelated to electrostatic attraction.

We use stainless-steel particles, which cannot become charged by friction. We recollect the fundamental properties of conductive metals: all excess charge resides on the surface, with charge density σ_e; electrons freely flow within, to, and from them – because there is no stable equilibrium point at which the potential energy has a minimum. Furthermore, in a uniform electrical field, the total force on an uncharged conductor is identically zero (see, e.g. Landau et al. 1982, for formal proofs).

There is no external energy source in this system to induce or sustain a net charge on the particles (e.g., the energy of our light source is too low for the photo-electric effect to occur on the metallic surfaces, and even if it were not, the result would be net repulsion, not attraction). Still, supposing that the particles are not completely neutral, the field around each particle of surface area A and charge density σ_e would have everywhere a potential energy, in the direction n normal to the surface E_n = 4πσ_e = − ∂Φ∕∂n, where Φ is the potential. The total charge Q on the sphere, inscribed within a Gaussian surface, is Q = σ_eA. It would require a very high surface charge density to overcome the potential energy of our small, heavy particles, $E_{\varphi }={\rho }_{\text{m},\text{p}}(g{h}_{\text{meas}}-\frac{1}{12}\pi d_{\text{p}}^3{u}_{\text{pz}})$, where h_meas is the height at which measurements are made which, along with all other variables, is defined in Sect. 2 of this paper.

The attractive force given by Coulomb’s law $F\propto \frac{Q_i{Q}_j}{d^2}$ between particles i and j is directed along the separation vector d between particles. By stark contrast, we find that the direction of the attractive phenomenon we observe in Figs. 8 and 9 is not along the separation vector, but strictly along the vertical component of it. There is indeed a quickening of the attractive behaviour in the vertical direction with decreasing vertical distance, however over rather great distances of up to 4 mm ~ 80−90 d_p, where the Coulomb force should be extremely weak, and we do not find any particle pairs at separations close to d_p (even though the positional accuracy of the measurements is much smaller than this scale, see Sect. 2), where the Coulomb forces actually, in general, operate.

5.3. Drafting and the consequence of Re_p

Particles that systematically drift in a PPD correspond to sizes of pebbles or boulders (mm-m), depending upon location or disc model assumptions. Such particles are very close to the transition given by Eq. (3) and in the MMSN, this transition is given by 9∕2λ = 6.4 cm ${\left(\frac{r}{\text{AU}}\right)}^{2.75}$, that is, 6 cm particles at 1 AU. In Fig. 5, we demonstrated that our experiments bracket this transition.

In either the Stokes or Epstein drag regime, an important property of flow around individual spheres is that it is symmetric, both in space and time (Guazzelli et al. 2011). As a pertinent example, Youdin & Goodman (2005) exploited this property to linearise theNavier-Stokes equations and derive an analytical dispersion relation for the SI, and refer to this simplification as a “terminal velocity approximation”. When removing the time dependence of the particle velocity, it is sometimes customary to also retract Navier’s name from the momentum balance equations and refer to them simply as “Stokes equations” or “Stoksian” flow, and in practical terms this mainly means that the drag law (be it “Stokes” or “Epstein”) is linear. A flow can be considered Stoksian when the shearing rate is small compared to the viscous stress and Re_p ≪ 1; in Sect. 2 we demonstrated that this is true in all of our experiments.

There are several highly studied regimes in particle-pair and collective-particle sedimentation phenomena that are due to flow symmetry breaking by increasing Re_p much higher than 1. A notable example is “wake attraction”, which occurs when the flow separates from the downstream side of a sphere to generate a wake, where the low-pressure conditions essentially create a vacuum and pull a closely aligned and relatively nearby (never more than 10 d_p downstream) particle into the wake. Various studies differ regarding to precisely which Reynolds number is required for this effect to result in particle clustering (depending upon whether it is a 2D or 3D system, whether the particles are tightly confined by walls, how high is ϕ, and whether the wake sheds vortexes or not). All studies (see e.g. Fortes et al. 1987; Feng et al. 1994; Kajishima & Takiguchi 2002; Uhlmann & Doychev 2014) are unanimous in stating that it requires moderate-to-high Reynolds numbers, meaning 50 < Re_p < 850, which are decidedly non-Stoksian values. Such relatively high values of Re_p are required because the mechanism quite literally involves a well-separated or even turbulent wake, for example of the types beautifully visualised in Tietjens & Prandtl (1957) and van Dyke (1982).

In addition to our experiments being in the incorrect regime to expect wake attraction, the phenomenology as it is formally studied is inconsistent with the details of our data. Wake attraction requires a very favourable alignment, and in Fig. 9 the aggregation effect occurs for both vertically aligned and misaligned particles – as expressed by how all of the curves appear the same for every value of horizontal separation. While wake attraction is a relatively long-range effect, it is not so long as what is found here; as mentioned above in Sect. 5.2, the separation at which the particle layers aggregate in this system is much greater than 10 d_p downstream.

Due to all of these considerations, we resist any speculation that the aggregation in this system is simply a manifestation of this well-known result, even though it may remind one of classical studies of strong flow-symmetry breaking.

5.4. Dragdissipation and consequence of drag regime (Kn)

We remind the reader that the two “representative” data sets we have presented are those with values of ϵ for which inertial coupling between the carrier and disperse phases is non-negligible. This criterion is considered important for reaching critical growth rates both in the formulation of a parallel shear flow with Keplerian gravitational potential (Youdin & Goodman 2005), and as in the present case of sedimentation with constant gravitational acceleration (Lambrechts et al. 2016). We achieve this criterion in two different ways: in kn1_3mb[DS2] (Kn ~ 1), we use fewer particles and lower gas pressure; in kng1_8mb[DS3] (Kn < 1), we use moreand slightly larger particles, with the gas pressure over twice the gas pressure of kn1_3mb[DS2].

In Sect. 4.4, we find that even for the same local number density conditions, drag is reduced only for the Kn < 1 data. We do not equate the aggregation behaviour we find in Figs. 8 and 9 with the lowered drag dissipation rate we report in Fig. 14. Our data unambiguously show that spontaneous mass concentration occurs even in the absence of a correlation between particle number density and mean settling rate. We believe this effect, when it happens, plays an auxiliary role, such as by increasing the magnitude of the baseline value of δu_z in the bottom panels in Figs. 8 and 9 and by shortening the timescale for the cross-correlation of number density to decay in the bottom panel of Fig. 12. Because of this contribution, particles find it easier to “catch up” to one another and are exchanged more rapidly between high and low-density regions.

As noted above, due to the high contrast in material densities between the two phases, high ϵ does not imply high ϕ. So, while increases in the local particle number density by several orders of magnitude proportionately augment ϵ, the volume filling factor ϕ remains very low. Therefore, the collective drag reduction we observe (only for Kn < 1, not for Kn ~ 1) does not come from augmenting the total surface area presented to the flow.

It is much more likely that here, as in Youdin & Johansen (2007) and Johansen & Youdin (2007), the collective momentum of the particles is slowing the gas locally, which in turn reduces the drag the particles feel. Unfortunately, we cannot visualise or measure the local gas flow simultaneously with the inertial particles. However, the dataset for which collective drag reduction does not occur, while active aggregation does, serves as proxy and supports this hypothesis.

6. The state of the art

We interpret the spontaneous aggregation behaviour we observe to fulfil the predictions of the “mass-loaded rain” simulations in Lambrechts et al. (2016). We emphasise, however, that the interpretation of the Lambrechts et al. (2016) simulations – along with attempts to grasp the physical meaning of the mathematical and numerical calculations that produce SI – is being constantly revisited by theory (see e.g. Squire & Hopkins 2018a).

While the interpretation of our results is guided by fundamental principles, there is at the same time very little pre-existing framework for understanding the detailed physics across scales: to our knowledge, multi-particle physics in the transition drag region, additionally with such low particle packing fractions (ϕ), has never been studied before.

We believe that, by virtue of the very low values of ϕ, we have removed the primary cause of inter-particle perturbations, namely long-range hydrodynamic interactions, which are deemed responsible for strong velocity and density fluctuations in fluidised beds or, equivalently, particle suspensions (Batchelor 1972; Guazzelli & Hinchdoychev 2011; Tee et al. 2008). Yet, the current state of the art in understanding the decay of the perturbation field around particles comes from investigations of low Kn systems where pure continuum fluid assumptions apply. We therefore highlight the question of how such perturbations decay in rarefied gas of Kn ≈ 1 or greater as a potentially impactful future direction for investigation.

Since the measurements of Cunningham (1910) were made, there has been significant progress in developing an experimentally verified unified theory for drag on individual objects in transitional flow (Karniadakis et al. 2005; Fissell et al. 2011). However, the results are quite dependent on the boundary conditions both at the particle-fluid interface and at relatively long range. The boundary-condition dependence of fluid drag presents notorious challenges in studying multi-particle systems, because the presence of neighbouring particles in the flow alters the shearing rate at particle-fluid interfaces. There is a storied history of attempts to understand how viscous drag is modified due to this fact, building primarily on the seminal works of Darcy (1856), Einstein (1905), Carman (1937), and Brinkman (1949). A drag force correction can be derived, and many authors have proposed variations on this theme (see, e.g. Hamdan 2016; Lasseux & Valdes-Parada 2017). However, all such proposed corrections universally take the form F_d = F_d,Stokes(1 + f(ϕ)), where F_d,Stokes and ϕ are defined the same as in this work (although often written in terms of the porosity, which is P = 1 − ϕ). This is necessarily the case, because drag by definition depends on the total surface area presented to the flow. Due to the vanishing values of ϕ in these experiments, such theories cannot explain the apparent drag reduction shown in the top panel of Fig. 13.

Whether collective drag reduction occurs, or whether it does not, depends upon Kn. The experimental result is nearly in agreement with the simulations presented in Lambrechts et al. (2016), who only considered particles in the Epstein regime and find a very weak effect of enhanced particle settling by particle swarms formed in a non-linear instability. There is perhaps a seed of insight here regarding the long-standing question of why strong clumping in the non-linear phase of SI is particle-size dependent (see e.g. Schaffer et al. 2018), and so future work to derive a precise explanation for the Kn dependence should be worthwhile. We suggest that an important factor to take into consideration when pursuing this line of inquiry is that the kinematic viscosity is equivalent to a momentum diffusivity, which increases inversely with decreasing gas density. See also the semi-analytical analysis in Lambrechts et al. (2016), in which it is shown that increasing viscosity has a damping effect on the toy model that is studied there. Squire & Hopkins (2018a) also address how differences in drag regime have an impact on the growth of resonant drag instability, which is a family of instabilities to which, they argue, SI belongs.

We point out that wake formation and associated attractive forces in rarefied gas is not a topic of previous investigation. Therefore, the possibility that we are the first to observe an effect related to wake attraction in a completely new regime cannot be completely dismissed. However, even if some form of wake-related behavior is operating in this system, it has a very different characteristic than traditional sedimentation. Drafting events in dilute 3D sedimentation generally do not lead to clustering, but to the renowned “draft, kiss, tumble” (DKT) effect which, as indicated by the name, involves a momentary but fleeting contact between particles (Fornari et al. 2016). As expressed above, we do not observe contact. Moreover, simulations of sedimentation show that isolated DTK events do not dominate the velocity statistics, but can be identified as a bump in the tail in the negative end of the settling velocity PDF; they do not dominate the mean because a very unique alignment is required. We would also be able to detect this bump on the left-hand side of the z-direction velocity PDFs in Appendix A. As a caveat, we do not have enough data here to study extremely rare events in the statistical tails.

We urge caution in making direct analogies to wake attraction as it has previously been studied primarily as a particle-pair process. We find tentative evidence in Sect. 4.2 (see also Appendix B) for a developing pressure phase lag with heuristic similarity to the mechanism proposed in Lin & Youdin (2017). It is well known that flow through densely packed porous beds results in a pressure drop (this is the operational principle in the right hand panel of Fig. 2, in which the packing fraction of the porous filter is maximum in order to create a steady low-pressure flow). However, in the experimental test chamber, low ϕ implies that there is no flow impedance to induce the local pressure gradient, in the classical sense of Darcy (1856) and Carman (1937).

We find more grounds for interpreting our results in the physical principles inherent to SI, than we do from classical sedimentation experiments and theory of suspension dynamics. We cannot find any major contradiction in the predictions of the SI and the results of the present study. However, reconciling volume-averaged, two-fluid models with discrete systems is in general an outstanding field of inquiry in sedimentation (Brennen 2005). In a real system, one must grapple with the fact that particles are discrete entities that also interact individually with free-flowing molecules. Therefore, our ongoing collaborative efforts include code development starting from gas and solid kinetics (Herminghaus & Mazza 2017), in order to test the constitutive equations of the macroscopic scale.

7. Summary

We constructed an experimental facility uniquely suited for the study of dust-grain aerodynamics under conditions similar to those in a PPD. The apparatus contains a suspension of solid spherical particles (of size 15–65 μm) settling under constant gravitational acceleration against a low-pressure gas stream. We calibrated the fluid state parameters and gas-velocity field in the absence of solid particles to confirm that the bulk flow is laminar and incompressible and that the conditions are steady, as discussed in Sects. 2 and 3.

We repeated particle-sedimentation experiments at four fixed pressures in the mbar range and featured the results from two experimental data sets, each of which representative of being on either side of the Stokes–Epstein drag law transition. The typical volumetrically averaged solid-to-gas mass density ratio of these two data sets was close to one.

Using a particle-tracking technique, we measured individual particle motions in three dimensions, within a cm from the centre of the chamber, within which distance the gas-flow velocity varies by only a few percent. The particle-tracking statistics reveals a complex particle-velocity structure that varies as a function of scale.

The statistics of relative velocity between particle pairs shows that, on average, particles approach one another when their vertical separation is small. While this effect depends strongly on vertical separation, it is independent of horizontal separation, at least within the range of the measurement volume. It is interesting to note that the vertical direction is the direction of the mean relative motion between the gas and the solid phases. This dynamic phenomenon operates in concert with variations in the mean particle density, measured at fixed position while particles advect through the measurement window, and therefore corresponds to the mean particle density variation in the vertical direction. Such a variation thus can be interpreted as particle-density waves.

In addition to the observation of convergent motions leading to particle-dense layers, we find enhanced mean settling velocities that depend upon the local particle number density. This effect is suggestive of collective particle-drag reduction. However, it is only clearly seen when the Knudsen number corresponds to nearly continuum flow conditions around the particles, and does not occur for the data set with Kn approaching free molecular flow. Since it is only present in the Kn < 1 experiments, yet both the Kn ~ 1 and Kn < 1 data exhibit complex convergent particle motions, the collective drag reduction might be neither exactly the same as, nor the cause of the particle aggregation behaviour; rather, it may be a secondary effect that occurs in systems where particle clustering is already emerging, but this secondary effect helps accelerate the clustering. See also Youdin & Johansen (2007) for an explanation of the role of energy dissipation by collective-drag forces in SI (Sect. 5 of that article).

The two observations – on-average approaching relative velocities and long-lasting particle density variations – indicate that the disperse phase of the system in question is unstable. The particle density variations are however extremely localised in scale and also illustrate a diffusive tendency, and so we find no indication that they are static structures. We interpret the strong directional dependence of the observed inhomogeneities to be a wavenumber selection process, such that the flow is sensitive to disturbances in the K_z direction, and either K_x and K_y are 0, or slightly larger than 2π∕w_meas, with w_meas the width of the measurement volume V_meas. Because we measured the gas flow to be laminar in the absence of particles, we have strong reason to suspect that the observed highly localised behaviour arises from the dynamic interaction between the two phases. We explore an appropriate parameter range to represent the diffuse particle distribution in PPDs, and thus our experiments provide evidence for the presence of particle-gas drag instabilities in an astrophysical context.

The most appropriate theoretical framework for interpreting the observed instability comes from astrophysical flows where so-called “streaming” or “drafting” instabilities (or “resonant drag”, more broadly) are theoretically predicted to occur and are considered a promising route to planetesimal formation. Such theories have not previously been directly tested. This work is the first to generate a real physical system that meets the essential physical assumptions of these models, most closely that presented in Lambrechts et al. (2016) of particle sedimentation in a non-rotating system. We subsequently find a spontaneous concentration phenomenon to self-organise. Specifically, in a pressure-driven flow with order-one mass loading of particles in which coupled inertial particles experience linear drag forces, trailing particles tend to draw near to leading particles, localised dust-to-gas density enhancements occur, and – for a certain drag regime – so does collective particle drag reduction.

Since the dilute-gas flow vessel used here is the first of its kind, we expect its existence to enable further novel laboratory-based investigations of particle aerodynamics in astrophysical flow. We have demonstrated the ability to achieve steady-state conditions for fixed dust-to-gas density ratio and for fixed Knudsen number; further experiments to systematically vary these two parameters independently will be an important future direction within the existing capabilities of the facility. Additional development of a moving measurement system will be required to assess how and when the instability saturates and if the growth rate depends upon varying the system parameters. Developing the capability to visualise the gas response to the inertial particles would also be a useful innovation.

Acknowledgements

We thank Hubert Klahr and Philip Armitage for helpful input and also acknowledge inciteful conversations with Greg Voth and Andrew Youdin. This work was supported by the German Deutsche Forschungsgemeinschaft, DFG collaborative research centre 963 Astrophysical Flow Instabilities and Turbulence, and by the grant “Bottlenecks for particle growth in turbulent aerosols” from the Knut and Alice Wallenberg Foundation (Dnr. KAW 2014.0048). H.L.C. also received support within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation. The anonymous reviewer’s constructive critique improved the manuscript. The mechanical drawings of the experimental apparatus are credited to Dr. Artur Kubitzek.

Appendix A: Intermittent statistical distributions

We present an overview of the particle behaviour by considering the particle velocity and acceleration PDFs. For each of the two representative data sets, Fig. A.1 shows in the left-hand panels the PDF of instantaneous velocity normalised by the rms velocity for the distribution,$u_i^{\star }\equiv u_i/{\langle u_i^2\rangle }^{1/2}$, for each of the i components. Similarly, the PDF of instantaneous acceleration, normalised by the rms acceleration, $a_i^{\star }\equiv a_i/{\langle a_i^2\rangle }^{1/2}$, is shown in the right-hand panels of the same figure. The colour and symbol coding in Fig. A.1 refers to the three different components,where z is vertical and x and y are horizontal. An equivalent Gaussian curve is over-plotted for comparison. We notice in the Kn ~ 1 case, that both the velocity and acceleration PDFs deviate from a normal distribution. The same is not true of the Kn < 1 case, since the velocities are nearly Gaussian, but the accelerations are clearly not. There are a number of possible reasons for the different appearance of these sets of distributions. In the Kn < 1 case, the number density is relatively high, and so it is possible that coordinated motions of the particles serve to reign in the excursions, even while the accelerations remain intermittent. The rms velocity, which can be read off of Table 2, is much higher, and so excursions may not be as apparent due to the normalisation, whereas the opposite is true in the horizontal component PDFs of velocity.

Fig. A.1. Velocity (left panels) and acceleration (right panels) for two representative data sets. Blue x’s, green plus marks, and red circles correspond to the x, y, and z components,respectively. All PDFs are normalised by their rms value. A Gaussian curve is overplotted for comparison.

Appendix B: Higher order moments of the relative velocity increment in axisymmetric form

In Sect. 4.2, the mean relative velocity is presented in component form and in an axisymmetric cylindrical coordinate system, with $r=\sqrt{x^2+y^2}$, and z the vertical component. Here we present higher order moments of the relative velocity, using the same geometry and conditioning. Figures B.1 and B.2 show the second-order and third-order moments of the relative vertical velocity, respectively, conditioned upon separation δr, plotted against vertical separation, δz, for each of the two representative data sets.

The second-order moment reveals a dependence of $\delta u_z^2$ on radial separation, δr, for both data sets: the placement of the curves on the y-axis decreases with descending radial separation. The spread in values of $\delta u_z^2$ is much larger than the flow-profile velocity variation shown in Fig. 7. Therefore, we take the trend in both panels of Fig. B.1 to indicate that the correlation in particle motion decays with radial separation. This result may also be consistent with a developing pressure phase lag, for example, particles at small radial separation are close to the lagging pressure-wave peak and those at larger separation respond to the pressure-wave trough.

Fig. B.1. Second-order moment of δuz vs. δz, conditioned upon horizontal separation δr. Top (bottom) panel: for Kn ~ 1 (Kn < 1) data. Conditioning and coloring same as for Fig. 9.

There is a notable difference between the two data sets: for the Kn ~ 1 data shown in the top panel, the second-order velocity difference in the z-direction turns towards smaller values for increasing δz. This indicates that, for separations larger than the wavelength of the unstable wave mode, the relative velocity will eventually return to zero. The same does not occur in the Kn < 1 data shown in the bottom panel.

The major difference between the two data sets is illustrated in Sect. 4.4: the Kn < 1 data shows a dependence of settling velocity on particle concentration but the Kn ~ 1 data does not. Only when this effect is present can it serve to increase the velocity dispersion and hence prevent the particle relative velocity to return to zero for large distance.

The second-order moment, being even, is unsigned and so it does not indicate direction. The apparent influence of the collective drag effect can therefore be better seen in the third-order moment, $\delta u_z^3$ vs δz, shown in Fig. B.2. The difference between the two data sets is striking. In the top panel, which shows the result for the Kn ~ 1 data, the third-order vertical velocity difference approaches zero relative velocity. In the bottom panel, which shows the result for the Kn < 1 data, there is an apparent turnover that depends upon δr and occurs around δz of 4–5 mm. The turnover obviously cannot be due to the gas flow profile because it would then be apparent in both data sets. Instead, it seems that we slightly pick up a non-zero horizontal wave mode. The implication is that the collective drag reduction effect is likely to increase the wave numbers k_x and k_y, and supports our interpretation of the wave-like nature of the clustering we report on in Sects. 4.2 and 4.3.

Fig. B.2. Third-order moment of δuz vs. δz, conditioned upon horizontal separation δr. Top (bottom) panel: for Kn ~ 1 (Kn < 1) data. Conditioning and coloring same as for Fig. 9.

For completeness, we also show the moments of the velocity difference in the horizontal direction in Figs. B.3 and B.4. We see that, despite the complex velocity field we observe, the relative velocity remains homogeneous in the horizontal spatial directions, up to third order. Therefore, this system meets the essential criteria for fluid instability: the base state must be statistically homogeneous in at least one spatial direction and the governing equations in equilibrium must be deterministic, and not time dependent; given enough time to develop, and a free energy source (here the differential gas-dust settling motion of ≈1m s⁻¹), the system should approach a new steady state, characterised by non-linear waves (Cross & Greenside 2009). We illustrated in Fig. 6 that the equilibrium vertical velocities are not time dependent. The right-hand panel of Fig. 6 shows that the fluid instability predicted by Lambrechts et al. (2016) to occur specifically in this experimental system, has had time to nearly fully develop. According to Fig. 2 of that publication, the growth reaches its saturation state in no more than 5–10 T_f, which is comparable to the travel time of particles in this system before reaching the measurement location. Figure 9 of this paper shows that the particles organise into density waves. Figures 10 and 11 show that, while structures (and voids) survive long enough for the particle density to remain correlated during their passage through the measurement window, the number density eventually shows the same correlation as average seeding density, implying that small-scale heterogeneities form due to growth of perturbations upon a steady base state. Finally, in this appendix, Figs. B.3 and B.4 illustrate that the deterministic equations of the new, unstable state remain homogeneous in two spatial directions.

Fig. B.3. Shown for Kn ~ 1 data. Top, middle, bottom panels: first-, second-, third-order moments of the velocity difference in the radial direction, conditioned upon vertical separation, as a function of radius.

Fig. B.4. Shown for Kn < 1 data. Top, middle, bottom panels: first-, second-, third-order moments of the velocity difference in the radial direction, conditioned upon vertical separation, as a function of radius.

Figures B.3 and B.4 also reinforce why the aggregation effect cannot be due to Coulomb forces, which would not cause zero radial component. Furthermore, it is abundantly clear from these figures that the particle flux is not controlled by any hypothetical anomaly outside the view of the measurement window.

References

All Tables

All Figures

Fig. 1. Particle fluidisation and sedimentation chamber, shown towards the right of the image, is supported on a heavy aluminium plate and fixed to an extruded aluminium profile. A gas line leads into the pressure-reducing filter at the bottom and vacuumtubing leading to a vacuum pump extends from the hollow gas expander at the top. Cables extending from the device connect pressure transducers to a data-logging system. On the optical table are three high-speed phantom V10 cameras, positioned to have overlapping field of view inside of the apparatus, and the signal generator used to synchronise the three cameras via external triggering. LED spotlights are mounted near the apparatus for use as backlighting.

In the text

Fig. 2. Cut-away views of components used to control the gas pressure and two-phase flow velocity. Left panel: hollow gas expansion unit where particles enter and slow when they reach the top of the sedimentation chamber. Right panel: apparatus base with pressure-reducing filter composed of sintered metal plates of porosity coefficient α = 10−12 m2. The funnel insert holds a fine wire mesh that is used as a particle staging platform; low-pressure air streaming through the mesh fluidises the particles.

In the text

Fig. 3. Flow profile calibration measurement using PIV, limited to the region of possible overlap with the LPT measurements, which are conducted within a ~1 cm3 region centred within the vessel with an accuracy of ~0.5 cm. Top panel: measured mean velocity (blue dots), and parabolic fit (purple line), limited to the central 2 cm of the flow vessel. Errorbars represent the standard deviation of the time-series PIV data. Bottom panel: calculation of the relative difference invelocity from the maximum centre-line velocity. The shaded region represents the uncertainty range due to error bars in the top panel.

In the text

	Fig. 4. The neighbourhood of the origin of the LPT measurements, where the black, red, and blue lines correspond to the lines of sight and centre point for cameras 0, 1, and 2, respectively. The thick purple lines delineate the measurement volume, Vmeas, defined by the intersecting of lines of sight from all three camera sensors.
In the text

Fig. 5. Flow conditions at the scale of particle size in the experiments. Black curve: the mean free path of gas molecules as a function of pressure λ(P). Coloured swaths denote Kn-dependent drag regimes. Grey dashed curve: Stokes–Epstein transition. The four vertical green bars represent each of the four data sets; the position of the bars denotes the pressure, their height denotes the spread in particle size for all experiments, and their width denotes the offset error in the pressure measurement. The orange subsegments indicate the particle sizes that have the steady-state velocity equaling the measured particle velocity in the experiment: up,ss = ⟨upz⟩. All the data sets can be considered transition flow, with the two at higher pressure further towards continuum and those at lower pressure approaching free molecular flow.

In the text

	Fig. 6. Left panel: height as a function of the ratio of particle velocity to the steady-state velocity upz ∕up,ss. Right panel: height as a function of time in units of the viscous relaxation timescale Tf. Particlesreach a steady state well before reaching the measurement window at ~ 50 cm.
In the text

	Fig. 7. Second-order moment of the longitudinal relative velocity component, Dll(l), versus l for all data sets. Purple solid curve corresponds to the magnitude of the velocity difference of the gas flow profile as a function of radial distance from the flow centre-line. The contribution to Dll from the variation in shear velocity in the gas phase is thus negligible.
In the text

	Fig. 8. Mean relative velocity components between particle pairs in units of m s−1, as a functionof vertical separation. The Kn ~ 1 data is shown in the top panel and the Kn < 1 data in the bottom panel. Particles aggregate for separations of 4 mm or less in both data sets.
In the text

	Fig. 9. Mean vertical relative velocity between particle pairs in units of m s−1, as a functionof vertical separation. The Kn ~ 1 data is shown in the top panel and the Kn < 1 data in the bottom panel. The different curves refer to increasingly large values in δr. The aggregation dynamics is the same at all radial (horizontal) separations, suggesting the formation of a particle layer.
In the text

Fig. 10. Conditional autocorrelation function of particle number density, normalised by the unconditional correlation, for the Kn ~ 1 case. Top panel: conditional on the dense case; bottom panel: conditional on the dilute case. Heterogeneity in the particle density is evident from the ratios deviating from unity. The timescale for the ratio to return to unity corresponds to a characteristic lengthscale of the clustering.

In the text

	Fig. 11. Same as Fig. 10, for the Kn < 1 case, showing generally similar clustering behaviour.
In the text

	Fig. 12. Temporal variation of the cross-correlation of particle number densities at different heights in the measurement volume. Top panel: the case of kn1_3mb [DS2]; bottom panel: kng1_8mb[DS3]. Shifting of the peaks with time and broadening of the shape indicate that the concentration evolves like an advection-diffusion process.
In the text

	Fig. 13. PDFs of z-direction velocity, conditional upon the number of particles in a 1-mm slab, for kn1_3mb [DS2] (top panel) and kng1_8mb [DS3] (bottom panel). Curves shift towards lower velocities with increasing number density in the lower panel.
In the text

	Fig. 14. Mean settling velocity in m s−1 for kn1_3mb [DS2] in blue, kng1_8mb[DS3] in red, and kng1_5mb[DS4] in green, conditional upon the number of particles in a 1 mm cylindrical slab. Drag dissipation is reduced for high density regions in the Kn < 1 data, yet notin the Kn ~ 1 data.
In the text

	Fig. A.1. Velocity (left panels) and acceleration (right panels) for two representative data sets. Blue x’s, green plus marks, and red circles correspond to the x, y, and z components,respectively. All PDFs are normalised by their rms value. A Gaussian curve is overplotted for comparison.
In the text

	Fig. B.1. Second-order moment of δuz vs. δz, conditioned upon horizontal separation δr. Top (bottom) panel: for Kn ~ 1 (Kn < 1) data. Conditioning and coloring same as for Fig. 9.
In the text

	Fig. B.2. Third-order moment of δuz vs. δz, conditioned upon horizontal separation δr. Top (bottom) panel: for Kn ~ 1 (Kn < 1) data. Conditioning and coloring same as for Fig. 9.
In the text

	Fig. B.3. Shown for Kn ~ 1 data. Top, middle, bottom panels: first-, second-, third-order moments of the velocity difference in the radial direction, conditioned upon vertical separation, as a function of radius.
In the text

	Fig. B.4. Shown for Kn < 1 data. Top, middle, bottom panels: first-, second-, third-order moments of the velocity difference in the radial direction, conditioned upon vertical separation, as a function of radius.
In the text