Issue 
A&A
Volume 622, February 2019



Article Number  A151  
Number of page(s)  20  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201833702  
Published online  13 February 2019 
Observation of aerodynamic instability in the flow of a particle stream in a dilute gas
^{1}
Max Planck Institute for Dynamics and SelfOrganization (MPIDS),
Am Faβberg 17,
37077 Göttingen,
Germany
email: holly.capelo@space.unibe.ch
^{2}
Lund Observatory, Department of Astronomy and Theoretical Physics, Lund University,
22100 Lund,
Sweden
^{3}
Institut für Geophysik und extraterrestrische Physik, Technische Universität Braunschweig,
Mendelssohnstr. 3,
38106 Braunschweig, Germany
^{4}
Institute for Nonlinear Dynamics, University of Göttingen,
37073 Göttingen,
Germany
^{5}
Laboratory of Atomic and SolidState Physics and Sibley School of Mechanical and Aerospace Engineering, Cornell University,
Ithaca,
NY 14853, USA
Received:
22
June
2018
Accepted:
2
December
2018
Forming macroscopic solid bodies in circumstellar discs requires local dust concentration levels significantly higher than the mean. Interactions of the dust particles with the gas must serve to augment local particle densities, and facilitate growth past barriers in the metre size range. Amongst a number of mechanisms that can amplify the local density of solids, aerodynamic streaming instability (SI) is one of the most promising. This work tests the physical assumptions of models that lead to SI in protoplanetary discs (PPDs). We conduct laboratory experiments in which we track the threedimensional motion of spherical solid particles fluidised in a lowpressure, laminar, incompressible, gas stream. The particle sizes span the Stokes–Epstein drag regime transition and the overall dusttogas mass density ratio, ϵ, is close to unity. A recently published study establishes the similarity of the laboratory flow to a simplified PPD model flow. We study velocity statistics and perform timeseries analysis of the advected flow to obtain experimental results suggesting an instability due to particlegas interaction: (i) there exist variations in particle concentration in the direction of the mean relative motion between the gas and the particles, that is the direction of the mean drag forces; (ii) the particles have a tendency to “catch up” to one another when they are in proximity; (iii) particle clumping occurs on very small scales, which implies local enhancements above the background ϵ by factors of several tens; (iv) the presence of these density enhancements occurs for a mean ϵ approaching or greater than 1; (v) we find evidence for collective particle drag reduction when the local particle number density becomes high and when the background gas pressure is high so that the drag is in the continuum regime. The experiments presented here are precedentsetting for observing SI under controlled conditions and may lead to a deeper understanding of how it operates in nature.
Key words: hydrodynamics / instabilities / turbulence / planets and satellites: formation / protoplanetary disks
© 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 kmsized 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 extrasolar systems (Seager & Lissauer 2010).
How firstgeneration 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 agedetermination techniques of starcluster sources with nearinfrared 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 largegrain 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, firstorder PPD models are idealgas 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 particleconcentration regions or lead to highenergy particle collisions that are ultimately destructive (Blum & Wurm 2008; Youdin 2010).
The statistical outcomes of dustaggregate 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 lowenergy collisions result in fractal growth, aggregates accrue momentum with increasing mass, such that further collisions may result in erosion or fragmentation, instead of continued bottomup growth (Zsom et al. 2010; Güttler et al. 2010; Bukhari Syed et al. 2017).
Drag dissipation from pressuresupported subKeplerian 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 metresize 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 dragforce 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 dusttogas 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: . 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 subdisc.
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 twophase flow, where twoway drag coupling is significant for ϵ = 1 and fourway 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 twofluid shearing sheet is linearly stable to SI when Ω = 0. Lambrechts et al. (2016) numerically evolved the same model equations, that is, incompressible twofluid dragcoupled mass and momentum balance, with terms relating to Keplerian rotation absent, and found the system to be unstable to nonlinear 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 subdisc evolution, the dust in the atmosphere of an embedded giant planet, or in the formation of chondrules. Additional works find dustdrag induced fluid instabilities to occur generically and in various stages of PPD formation and evolution (e.g. Bate & LorénAguilar 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 longstanding metrebarrier problem. However, smallscale 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 angularmomentum 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 dustgrain draginduced aerodynamic instabilities, which could eventually lead to turbulence, in sufficiently dilute scaleddown 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 particletracking data analysis exhibits particle clustering and modified particle dynamics on small scales, which we interpret as a collective particledrag 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 Epsteindrag 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 followup 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 pressurereducing filter composed of porous material. The filter creates a pressure drop ΔP = 9.9 × 10^{5} 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 steadystate pressure and particleseeding density are controlled and we record the particle positions from three viewing angles in order to reconstruct their trajectories.
The particles are preloaded onto a particle seeding platform, which is composed of a fine wire mesh, through which air can pass – but particles cannot (see righthand 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 lefthand 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 particleseeding platform, where they reenter in the flow andcirculate. The result of the particle circulation is that the centreline of the flow has a persistent upwarddirected 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 steadystate pressure conditions.
The particle positions are tracked in time using threedimensional 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 timeresolved 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 pressurereducing 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 datalogging system. On the optical table are three highspeed 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. 

Open with DEXTER 
2.2. Flow profile
The carrierphase 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 timeseries 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 twodimensional measurements of a single slice through the centreline 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 μmsized 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 steadystate operational pressure of ~10 mbar.
The measurements were made using one highspeed camera placed perpendicularly to a laser sheet^{1} from a single cavity Qswitched 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^{2} to perform the image preprocessing and processing. The preprocessing consisted of the following steps: (i) mean background image subtraction; (ii) local nonlinear subtraction; (iii) slidingminimum filter over a scale of 64 pixels; (iv) linear 3 × 3 pixel Gaussian smoothing filter; (v) intensity renormalisation filter. The processing involved interrogating preprocessed image pairs with a 3pass 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.
Lineprofiles of the twodimensional vector fields in the PIV timeseries data demonstrate that the gas flow has a cylindrical Pouseille profile, with a centreline velocity of 1.45 m s^{−1}. 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 containerscale 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.
Cutaway views of components used to control the gas pressure and twophase 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 pressurereducing filter composed of sintered metal plates of porosity coefficient α = 10^{−12} m^{2}. The funnel insert holds a fine wire mesh that is used as a particle staging platform; lowpressure air streaming through the mesh fluidises the particles. 

Open with DEXTER 
2.3. Inertial particle motion
In the past decades, traditional Eulerian flowmeasurement techniques have been complemented by new methods to study individual flow elements in their Lagrangian frame. Using multiple synchronised highspeed 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 turbulentlike motions, if an instability develops. We therefore employ the same particletracking technique used in previous studies of turbulent flow, with customisation for the present setup.
We acquired simultaneous shadowgraph images with three highspeed 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^{−1}. 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 leastsquare minimisation to the triangulation errors on a mapping between sensor position and realworld coordinates.The laboratoryframe coordinates are known via acalibration procedure in which the three cameras shown in Fig. 1 were focussed on a gridpattern 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 sensortorealworld 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 twoslit camera model goodnessoffit by performing the reconstruction on the calibration images. We subtract the realworld 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 threedimensional 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 cm^{3} 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 timeseries PIV data. Bottom panel: calculation of the relative difference invelocity from the maximum centreline velocity. The shaded region represents the uncertainty range due to error bars in the top panel. 

Open with DEXTER 
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, V_{meas}, defined by the intersecting of lines of sight from all three camera sensors. 

Open with DEXTER 
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 6s LPT recording, processed in the manner previously stated. A single data set includes several repetitions of the experiment () at a given gas pressure. To perform a statistical analysis, we combine the threedimensional 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 lengthscales, 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 lowRe 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 roomtemperature in our vacuum system is μ_{air} = 1.8 × 10^{−5} kg m^{−1} s^{−1}, which gives a kinematic viscosity of η_{air} = 2.0–7.9 × 10^{−3} m^{2} s^{−1} for the pressure range 200–800 Pa. Therefore, for the pipe flow on the container scale, 15–45, and for flow at the scale of an individual particle with diameter d_{p}, Re_{p} = d_{p}v_{rms,z}∕η ≪ 1, where D_{pipe} is the pipe diameter, m s^{−1} is the average velocity in the pipe, d_{p} is the particle diameter, and v_{rms,z} is the rootmeansquare (rms) velocity of the particles in the zdirection. 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, (1)
where V_{meas} is the volume of the measurement region, and spatially averaged mass loading (2)
where the material density of the stainless steel particles used in the experiments is ρ_{m,p} = 8050 kg m^{−3}. We also calculate the momentum diffusion time across the mean interparticle distance, t_{d} = n^{−2∕3}∕η, to demonstratethat the experiments are in a tight particlecoupling 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 piecewise draglaw transition for Re_{p} < 1, (3)
where δu = u_{g} − u_{pz} is the difference between the zcomponent particle velocity, u_{pz}, and the gas velocity, u_{g} = 1.45 m s^{−1}, 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: (4)
Further subdivisions in flow regime use the following as approximate guidelines (see, e.g. Allen & Raabe 1985; Karniadakis et al. 2005):
Figure 5 shows λ as a function of pressure, with the Kndependent regimes colorcoded for reference. Particles with sizes used in our experiments are represented by vertical green bars spanning 15 < d_{p} < 65 μm, with orange subsegments 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 Kndependent Cunningham correction factor (Cunningham 1910) to Stokes drag: (5)
with . 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^{3}. 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 zdirection, denoted and also listed in Table 1, and demonstrate narrow dispersion in particle size. The quantity 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,
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 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 steadystate 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 timedependent term to be neglected according to Eq. (6).
Parameters of four data sets collected at different gas pressure (mbar), listed in the second column.
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 Kndependent 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 steadystate velocity equaling the measured particle velocity in the experiment: u_{p,ss} = ⟨u_{pz}⟩. 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. 

Open with DEXTER 
Fig. 6.
Left panel: height as a function of the ratio of particle velocity to the steadystate velocity u_{pz} ∕u_{p,ss}. Right panel: height as a function of time in units of the viscous relaxation timescale T_{f}. Particlesreach a steady state well before reaching the measurement window at ~ 50 cm. 

Open with DEXTER 
3.2. Summary statistics
The particle velocities, in each of the data sets listed in Table 1, are by far dominant in the verticalcomponent, with 0.3 m s^{−1} ≲ ⟨u_{pz}⟩ ≲ 0.5 m s^{−1}. By comparison, the horizontal components are 0.001 ≲ ⟨u_{px,y}⟩ ≲ 0.005 m s^{−1}. 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^{−1}.
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 , 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 gasvelocity fluctuations at the flow centreline of 0.0098 m s^{−1} 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.
Componentwise rms velocities in units of m s^{−1}.
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: (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 socalled Eulerian secondorder 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, (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 solidline curve is plotted to show the difference in gas velocity between a point on the vessel’s centreline 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.
Secondorder moment of the longitudinal relative velocity component, D_{ll}(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 centreline. The contribution to D_{ll} from the variation in shear velocity in the gas phase is thus negligible. 

Open with DEXTER 
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 zaxis is the same as in the laboratory frame and the rcoordinate is in the horizontal plane. The origin, however, is not fixed, but set on any given particle. We study the relative velocity (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. 

Open with DEXTER 
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. 

Open with DEXTER 
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 twoposition, twotime correlation function (10)
where and 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^{−1} 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 fixedprobe 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 singlelocation, twotime autocorrelation function of n′ at times t and t + τ (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 , and the other for it being less than 0.5 times the mean number density, denoted . 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 and , for the Kn ~ 1 data. The densecase ratio, , is always less than one for all z. The converse is true for the dilutecase ratio, . To gain intuition for this result, we consider a limiting case: if the particles were arranged in a perfectly homogeneous fashion with the interparticle 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 nonuniformly 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 decorrelate 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 decorrelate 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^{−1}. 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 crosscorrelation 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 advectiondiffusion process. The larger value of C_{nn}(τ, δz) at long timelag 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 crosscorrelation peaks corresponds to the time the particles take to move through the distance dz at the mean velocity. The broadening of the crosscorrelation 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 decorrelation of the particle number density observed in the bottom panel.
This broadening of the crosscorrelation 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^{−1}. Consequently, any initial inhomogeneities in particle concentration would have been smeared out, judging from the decorrelation 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. 

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

Open with DEXTER 
Fig. 12.
Temporal variation of the crosscorrelation 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 advectiondiffusion process. 

Open with DEXTER 
4.4. Evidence for collective continuumdrag 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 particledrag reduction. That we do not observe this behaviour for the lower gasdensity 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 zdirection velocity, conditional upon the number of particles in a 1mm 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. 

Open with DEXTER 
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. 

Open with DEXTER 
5. Discussion
5.1. Interparticle 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 , depending on the gas pressure of our experiment. Therefore, we can achieve relatively high average solidtogas ratio, ϵ ≈ 1, expressed by Eq. (2), even with extremely low average filling factors, 10^{−7} < ϕ < 10^{−6}, 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 materialdensity contrast between rock, with ρ_{rock,m} ~ 3 g cm^{−3}, and gas is . For a dusttogas ratio of unity in the midplane, we obtain a volume filling factor ϕ ≈ 4.3 × 10^{−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 selfgravity.
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: (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^{−1}, 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^{−1}, which we adopt as V_{rel}. We obtain Z ≈ 10^{−2} collisions s^{−1} 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 overestimate, 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 fluiddynamical in nature and unrelated to electrostatic attraction.
We use stainlesssteel 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 photoelectric 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 = σ_{e}A. It would require a very high surface charge density to overcome the potential energy of our small, heavy particles, , 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 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 (mmm), 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 , 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 theNavierStokes 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 particlepair and collectiveparticle 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 lowpressure 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 moderatetohigh Reynolds numbers, meaning 50 < Re_{p} < 850, which are decidedly nonStoksian values. Such relatively high values of Re_{p} are required because the mechanism quite literally involves a wellseparated 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 longrange 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 wellknown result, even though it may remind one of classical studies of strong flowsymmetry 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 nonnegligible. 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 crosscorrelation 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 lowdensity 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 “massloaded 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 preexisting framework for understanding the detailed physics across scales: to our knowledge, multiparticle 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 interparticle perturbations, namely longrange 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 particlefluid interface and at relatively long range. The boundarycondition dependence of fluid drag presents notorious challenges in studying multiparticle systems, because the presence of neighbouring particles in the flow alters the shearing rate at particlefluid 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 & ValdesParada 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 nonlinear instability. There is perhaps a seed of insight here regarding the longstanding question of why strong clumping in the nonlinear phase of SI is particlesize 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 semianalytical 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 wakerelated 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 lefthand side of the zdirection 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 particlepair 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 lowpressure 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 volumeaveraged, twofluid 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 freeflowing 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 dustgrain 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 lowpressure gas stream. We calibrated the fluid state parameters and gasvelocity 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 particlesedimentation 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 solidtogas mass density ratio of these two data sets was close to one.
Using a particletracking technique, we measured individual particle motions in three dimensions, within a cm from the centre of the chamber, within which distance the gasflow velocity varies by only a few percent. The particletracking statistics reveals a complex particlevelocity 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 particledensity waves.
In addition to the observation of convergent motions leading to particledense layers, we find enhanced mean settling velocities that depend upon the local particle number density. This effect is suggestive of collective particledrag 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 collectivedrag forces in SI (Sect. 5 of that article).
The two observations – onaverage approaching relative velocities and longlasting 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 particlegas drag instabilities in an astrophysical context.
The most appropriate theoretical framework for interpreting the observed instability comes from astrophysical flows where socalled “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 nonrotating system. We subsequently find a spontaneous concentration phenomenon to selforganise. Specifically, in a pressuredriven flow with orderone mass loading of particles in which coupled inertial particles experience linear drag forces, trailing particles tend to draw near to leading particles, localised dusttogas density enhancements occur, and – for a certain drag regime – so does collective particle drag reduction.
Since the dilutegas flow vessel used here is the first of its kind, we expect its existence to enable further novel laboratorybased investigations of particle aerodynamics in astrophysical flow. We have demonstrated the ability to achieve steadystate conditions for fixed dusttogas 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 lefthand panels the PDF of instantaneous velocity normalised by the rms velocity for the distribution,, for each of the i components. Similarly, the PDF of instantaneous acceleration, normalised by the rms acceleration, , is shown in the righthand 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 overplotted 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. 

Open with DEXTER 
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 , 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 secondorder and thirdorder 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 secondorder moment reveals a dependence of on radial separation, δr, for both data sets: the placement of the curves on the yaxis decreases with descending radial separation. The spread in values of is much larger than the flowprofile 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 pressurewave peak and those at larger separation respond to the pressurewave trough.
Fig. B.1.
Secondorder moment of δu_{z} vs. δz, conditioned upon horizontal separation δr. Top (bottom) panel: for Kn ~ 1 (Kn < 1) data. Conditioning and coloring same as for Fig. 9. 

Open with DEXTER 
There is a notable difference between the two data sets: for the Kn ~ 1 data shown in the top panel, the secondorder velocity difference in the zdirection 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 secondorder 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 thirdorder moment, 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 thirdorder 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 nonzero 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 wavelike nature of the clustering we report on in Sects. 4.2 and 4.3.
Fig. B.2.
Thirdorder moment of δu_{z} vs. δz, conditioned upon horizontal separation δr. Top (bottom) panel: for Kn ~ 1 (Kn < 1) data. Conditioning and coloring same as for Fig. 9. 

Open with DEXTER 
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 gasdust settling motion of ≈1m s^{−1}), the system should approach a new steady state, characterised by nonlinear waves (Cross & Greenside 2009). We illustrated in Fig. 6 that the equilibrium vertical velocities are not time dependent. The righthand 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 smallscale 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, thirdorder moments of the velocity difference in the radial direction, conditioned upon vertical separation, as a function of radius. 

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

Open with DEXTER 
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
 Allen, M. D., & Raabe, O. G. 1985, Aerosol Sci. Technol., 4, 269 [NASA ADS] [CrossRef] [Google Scholar]
 Ansdell, M., Williams, J. P., Manara, C. F., et al. 2017, AJ, 153, 240 [NASA ADS] [CrossRef] [Google Scholar]
 Armitage, P. J. 2010, Astrophysics of Planet Formation (Cambridge: Cambridge University Press) [Google Scholar]
 Bai, X.N., & Stone, J. M. 2010a, ApJ, 722, 1437 [NASA ADS] [CrossRef] [Google Scholar]
 Bai, X.N., & Stone, J. M. 2010b, ApJS, 190, 297 [NASA ADS] [CrossRef] [Google Scholar]
 Bai, X.N., & Stone, J. M. 2010c, ApJ, 722, L220 [NASA ADS] [CrossRef] [Google Scholar]
 Balachandar, S., & Eaton, J. K. 2010, Ann. Rev. Fluid Mech., 42, 111 [NASA ADS] [CrossRef] [Google Scholar]
 Batchelor, G. K. 1972, J. Fluid Mech., 52, 245 [NASA ADS] [CrossRef] [Google Scholar]
 Bate, M. R., & LorénAguilar, P. 2017, MNRAS, 465, 1089 [NASA ADS] [CrossRef] [Google Scholar]
 Bell, C. P. M., Naylor, T., Mayne, N. J., Jeffries, R. D., & Littlefair, S. P. 2013, MNRAS, 434, 806 [NASA ADS] [CrossRef] [Google Scholar]
 Birnstiel, T., Andrews, S. M., & Ercolano, B. 2012, A&A, 544, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bitsch, B., Johansen, A., Lambrechts, M., & Morbidelli, A. 2015, A&A, 575, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Blum, J., & Wurm, G. 2000, Icarus, 143, 138 [NASA ADS] [CrossRef] [Google Scholar]
 Blum, J., & Wurm, G. 2008, ARA&A, 46, 21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Blum, J., Wurm, G., Poppe, T., & Heim, L.O. 1998, Earth Moon Planets, 80, 285 [NASA ADS] [CrossRef] [Google Scholar]
 Blum, J., Wurm, G., Poppe, T., et al. 1999, Meas. Sci. Technol., 10, 836 [NASA ADS] [CrossRef] [Google Scholar]
 Blum, J., Gundlach, B., Mühle, S., & TrigoRodriguez, J. M. 2014, Icarus, 235, 156 [NASA ADS] [CrossRef] [Google Scholar]
 Bourgoin, M. 2006, Science, 311, 835 [NASA ADS] [CrossRef] [Google Scholar]
 Brauer, F., Dullemond, C. P., Johansen, A., et al. 2007, A&A, 469, 1169 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brennen, C. E. 2005, Fundamentals of Multiphase Flow (Cambridge: Cambridge University Press) [Google Scholar]
 Brinkman, H. C. 1949, Appl. Sci. Res., 1, 27 [Google Scholar]
 Bukhari Syed, M., Blum, J., Wahlberg Jansson, K., & Johansen, A. 2017, ApJ, 834, 145 [NASA ADS] [CrossRef] [Google Scholar]
 Capelo, H. L. 2018, Ph.D. Thesis, GeorgAugustUniversität Göttingen, https://ediss.unigoettingen.de/handle/11858/0017350000 002EE4020 [Google Scholar]
 Carman, P. 1937, Chem. Eng. Res. Des., 75, S32 [Google Scholar]
 Carrera, D., Johansen, A., & Davies, M. B. 2015, A&A, 579, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cross, M., & Greenside, H. 2009, Pattern Formation and Dynamics in Nonequilibrium Systems (Cambridge: Cambridge University Press) [Google Scholar]
 Cunningham, E. 1910, Proc. R. Soc. London Ser. A, 83, 357 [NASA ADS] [CrossRef] [Google Scholar]
 Darcy, H. 1856, Les Fontaines Publiques de la Ville de Dijon (Libraire des corps imperiaux des ponts et chaussées et des mines) [Google Scholar]
 Dominik, C., & Tielens, A. G. G. M. 1997, ApJ, 480, 647 [NASA ADS] [CrossRef] [Google Scholar]
 Drążkowska, J., & Alibert, Y. 2017, A&A, 608, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dra̧żkowska, J., & Dullemond, C. P. 2014, A&A, 572, A78 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dubrulle, B., Dauchot, O., Daviaud, F., et al. 2005, Phys. Fluids, 17, 095103 [NASA ADS] [CrossRef] [Google Scholar]
 Einstein, A. 1905, Ph.D. Thesis, Universität Zürich, Switzerland [Google Scholar]
 Epstein, P. S. 1924, Phys. Rev., 23, 710 [NASA ADS] [CrossRef] [Google Scholar]
 Feng, J., Hu, H. H., & Joseph, D. D. 1994, J. Fluid Mech., 261, 95 [NASA ADS] [CrossRef] [Google Scholar]
 Fissell, W. H., Conlisk, A. T., Datta, S., et al. 2011, Microfluid. Nanofluidics, 10, 425 [Google Scholar]
 Fornari, W., Picano, F., & Brandt, L. 2016, J. Fluid Mech., 788, 640–69 [NASA ADS] [CrossRef] [Google Scholar]
 Fortes, A. F., Joseph, D. D., & Lundgren, T. S. 1987, J. Fluid Mech., 177, 467 [NASA ADS] [CrossRef] [Google Scholar]
 Gibert, M., Xu, H., & Bodenschatz, E. 2010, EPL, 90, 64005 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Goldreich, P., & Ward, W. R. 1973, ApJ, 183, 1051 [NASA ADS] [CrossRef] [Google Scholar]
 Guazzelli, E., & Hinch, J. 2011, Ann. Rev. Fluid Mech., 43, 97 [NASA ADS] [CrossRef] [Google Scholar]
 Guazzelli, É., & Hinchdoychev, J. 2011, Ann. Rev. Fluid Mech., 43, 97 [NASA ADS] [CrossRef] [Google Scholar]
 Guazzelli, É., Morris, J., & Pic, S. 2011, A Physical Introduction to Suspension Dynamics, Cambridge Texts in Applied Mathematics (Cambridge: Cambridge University Press) [Google Scholar]
 Güttler, C., Krause, M., Geretshauser, R., Speith, R., & Blum, J. 2009, AIP Conf. Ser., 1145, 941 [NASA ADS] [CrossRef] [Google Scholar]
 Güttler, C., Blum, J., Zsom, A., Ormel, C. W., & Dullemond, C. P. 2010, A&A, 513, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Haisch, Jr. K. E., Lada, E. A., & Lada, C. J. 2001, ApJ, 553, L153 [NASA ADS] [CrossRef] [Google Scholar]
 Hamdan, M. 2016, Int. J. Eng. Res. Appl., 6, 41 [Google Scholar]
 Hayashi, C., Nakazawa, K., & Nakagawa, Y. 1985, in Protostars and Planets II, eds. D. C. Black & M. S. Matthews (Tucson, AZ: University of Arizona Press), 1100 [Google Scholar]
 Herminghaus, S., & Mazza, M. G. 2017, Soft Matter, 12, 1 [Google Scholar]
 Hersant, F. B., & Hure, J. M. 2005, A&A, 429, 531 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jacquet, E., Balbus, S., & Latter, H. 2011, MNRAS, 415, 3591 [NASA ADS] [CrossRef] [Google Scholar]
 Johansen, A., & Youdin, A. 2007, ApJ, 662, 627 [NASA ADS] [CrossRef] [Google Scholar]
 Johansen, A., Youdin, A., & Mac Low M.M. 2009, ApJ, 704, L75 [NASA ADS] [CrossRef] [Google Scholar]
 Johansen, A., Youdin, A. N., & Lithwick, Y. 2012, A&A, 537, A125 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Johansen, A., Blum, J., Tanaka, H., et al. 2014, Protostars and Planets VI (Tucson, AZ: University of Arizona Press), 547 [Google Scholar]
 Kajishima, T., & Takiguchi, S. 2002, Int. J. Heat Fluid Flow, 23, 639 [Google Scholar]
 Karniadakis, G., Beskok, A., & NR, A. 2005, MicroFlows and Nanoflows  Fundamentals and Simulation (New York: Springer) [Google Scholar]
 Kowalik, K., Hanasz, M., Wóltański, D., & Gawryszczak, A. 2013, MNRAS, 434, 1460 [NASA ADS] [CrossRef] [Google Scholar]
 La Porta, A., Voth, G. A., Crawford, A. M., Alexander, J., & Bodenschatz, E. 2001, Nature, 401, 1017 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Lambrechts, M., Johansen, A., Capelo, H. L., Blum, J., & Bodenschatz, E. 2016, A&A, 591, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Landau, L., Lifschitz, E., & Pitaevskii, L. 1982, Electrodynamics of Continuous Media: Volume 8 Course of Theoretical Physics (Burlington, MA: Elsevier Butterworth) [Google Scholar]
 Lasseux, D., & ValdesParada, F. 2017, Comptes Rendus Mécanique, 345, 660 [Google Scholar]
 Lee, N., Williams, J. P., & Cieza, L. A. 2011, ApJ, 736, 135 [NASA ADS] [CrossRef] [Google Scholar]
 Lin, M.K. & Youdin, A. N. 2017, ApJ, 849, 129 [NASA ADS] [CrossRef] [Google Scholar]
 Marchioli, C. 2017, Collective Dynamics of Particles: From Viscous to Turbulent Flows, Book CISM International Centre for Mechanical Sciences, 576 (Cham, Switzerland: Springer) [Google Scholar]
 Miniati, F. 2010, J. Comput. Phys., 229, 3916 [NASA ADS] [CrossRef] [Google Scholar]
 Morbidelli, A., Lunine, J. I., O’Brien, D. P., Raymond, S. N., & Walsh, K. J. 2012, Ann. Rev. Earth Planet. Sci., 40, 251 [NASA ADS] [CrossRef] [Google Scholar]
 Mordant, N., Crawford, A. M., & Bodenschatz, E. 2004, Physica D Nonlinear Phenomena, 193, 245 [NASA ADS] [CrossRef] [Google Scholar]
 Nakagawa, Y., Nakazawa, K., & Hayashi, C. 1981, Icarus, 45, 517 [NASA ADS] [CrossRef] [Google Scholar]
 Ouellette, N. T., Xu, H., & Bodenschatz, E. 2006, Exp. Fluids, 40, 301 [Google Scholar]
 Pope, S. B. 2000, Turbulent Flows (Cambridge, England: Cambridge University Press) [Google Scholar]
 Safronov, V. S. 1969, Evolution of the Protoplanetary Cloud and Formation of the Earth and Planets (Moscow: Nauka), English translation: NASA TTF677 (1972) [Google Scholar]
 Saw, E.W., Shaw, R. A., Salazar, J. P. L. C., & Collins, L. R. 2012, New J. Phys., 14, 105031 [NASA ADS] [CrossRef] [Google Scholar]
 Saw, E.W., Bewley, G. P., Bodenschatz, E., Sankar Ray, S., & Bec, J. 2014, Phys. Fluids, 26, 111702 [NASA ADS] [CrossRef] [Google Scholar]
 Schäfer, U., Yang, C.C., & Johansen, A. 2017, A&A, 597, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schaffer, N., Yang, C.C., & Johansen, A. 2018, A&A, 618, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schoonenberg,D. & Ormel, C. W. 2017, A&A, 602, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Seager, S., & Lissauer, J. J. 2010, Introduction to Exoplanets (Tucson, AZ: University of Arizona Press), 3 [Google Scholar]
 Sharipov, F. 2011, J. Phys. Chem. Ref. Data, 40, 023101 [NASA ADS] [CrossRef] [Google Scholar]
 Simon, J. B., Armitage, P. J., Youdin, A. N., & Li, R. 2017, ApJ, 847, L12 [NASA ADS] [CrossRef] [Google Scholar]
 Squire, J., & Hopkins, P. F. 2018a, MNRAS, 477, 5011 [NASA ADS] [CrossRef] [Google Scholar]
 Squire, J., & Hopkins, P. F. 2018b, ApJ, 856, L15 [NASA ADS] [CrossRef] [Google Scholar]
 Tee, S.Y., Mucha, P. J., Brenner, M. P., & Weitz, D. A. 2008, J. Fluid Mech., 596, 467 [NASA ADS] [CrossRef] [Google Scholar]
 Tennekes, T., & Lumley, J. L. 1972, A First Course in Turbulence (Cambridge, USA: The MIT Press) [Google Scholar]
 Testi, L., Birnstiel, T., Ricci, L., et al. 2014, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning (Tucson, AZ: University of Arizona Press), 339 [Google Scholar]
 Tietjens, O., & Prandtl, L. 1957, Applied Hydro and Aeromechanics: Based on Lectures of L. Prandtl No. v. 2 (New York: Dover Publications) [Google Scholar]
 Turner, N. J., Fromang, S., Gammie, C., et al. 2014, Protostars and Planets VI (Tucson, AZ: University of Arizona Press), 411 [Google Scholar]
 Uhlmann, M., & Doychev, T. 2014, J. Fluid Mech., 752, 310 [NASA ADS] [CrossRef] [Google Scholar]
 van Dyke, M. 1982, An Album of Fluid Motion (Stanford, CA: Parabolic Press) [Google Scholar]
 Weidenschilling, S. J. 1977a, MNRAS, 180, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Weidenschilling, S. J. 1977b, Ap&SS, 51, 153 [NASA ADS] [CrossRef] [Google Scholar]
 Whipple, F. L. 1972, in From Plasma to Planet, ed. A. Elvius (New York: Wiley Interscience Division), 211 [Google Scholar]
 Wyatt, M. C. 2008, ARA&A, 46, 339 [NASA ADS] [CrossRef] [Google Scholar]
 Xu, H. 2008, Meas. Sci. Technol., 19, 075105 [NASA ADS] [CrossRef] [Google Scholar]
 Xu, H., & Bodenschatz, E. 2008, Physica D: Nonlinear Phenomena, 237, 2095 [NASA ADS] [CrossRef] [Google Scholar]
 Xu, H., Bourgoin, M., Ouellette, N. T., & Bodenschatz, E. 2006, Phys. Rev. Lett., 96, 024503 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Xu, H., Ouellette, N. T., & Bodenschatz, E. 2007, Phys. Rev. Lett., 98, 050201 [NASA ADS] [CrossRef] [Google Scholar]
 Yang, C.C., & Johansen, A. 2014, ApJ, 792, 86 [NASA ADS] [CrossRef] [Google Scholar]
 Youdin, A. N. 2010, EAS Pub. Ser., 41, 187 [Google Scholar]
 Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459 [NASA ADS] [CrossRef] [Google Scholar]
 Youdin, A., & Johansen, A. 2007, ApJ, 662, 613 [NASA ADS] [CrossRef] [Google Scholar]
 Zsom, A., Ormel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, A&A, 513, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
This camera and illumination scheme is different than the one shown in Fig. 1, since the two measurement techniques, PIV and LPT, require a different setup. The tracer (smoke, 1 μm) particles used in the carrier phase flow calibration should not be confused with the inertial (steel, 25–65 μm) particles used in the LPT experiments.
Assuming a drag law with abrupt transition might lead to slightly overestimating the particle size and hence deriving somewhat different values for the other parameters that depend on d_{p}. However, the two approaches to calculating d_{p} do not yield qualitatively different results in terms of the flow or massloading regimes that result.
All Tables
Parameters of four data sets collected at different gas pressure (mbar), listed in the second column.
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 pressurereducing 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 datalogging system. On the optical table are three highspeed 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. 

Open with DEXTER  
In the text 
Fig. 2.
Cutaway views of components used to control the gas pressure and twophase 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 pressurereducing filter composed of sintered metal plates of porosity coefficient α = 10^{−12} m^{2}. The funnel insert holds a fine wire mesh that is used as a particle staging platform; lowpressure air streaming through the mesh fluidises the particles. 

Open with DEXTER  
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 cm^{3} 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 timeseries PIV data. Bottom panel: calculation of the relative difference invelocity from the maximum centreline velocity. The shaded region represents the uncertainty range due to error bars in the top panel. 

Open with DEXTER  
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, V_{meas}, defined by the intersecting of lines of sight from all three camera sensors. 

Open with DEXTER  
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 Kndependent 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 steadystate velocity equaling the measured particle velocity in the experiment: u_{p,ss} = ⟨u_{pz}⟩. 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. 

Open with DEXTER  
In the text 
Fig. 6.
Left panel: height as a function of the ratio of particle velocity to the steadystate velocity u_{pz} ∕u_{p,ss}. Right panel: height as a function of time in units of the viscous relaxation timescale T_{f}. Particlesreach a steady state well before reaching the measurement window at ~ 50 cm. 

Open with DEXTER  
In the text 
Fig. 7.
Secondorder moment of the longitudinal relative velocity component, D_{ll}(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 centreline. The contribution to D_{ll} from the variation in shear velocity in the gas phase is thus negligible. 

Open with DEXTER  
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. 

Open with DEXTER  
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. 

Open with DEXTER  
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. 

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

Open with DEXTER  
In the text 
Fig. 12.
Temporal variation of the crosscorrelation 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 advectiondiffusion process. 

Open with DEXTER  
In the text 
Fig. 13.
PDFs of zdirection velocity, conditional upon the number of particles in a 1mm 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. 

Open with DEXTER  
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. 

Open with DEXTER  
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. 

Open with DEXTER  
In the text 
Fig. B.1.
Secondorder moment of δu_{z} vs. δz, conditioned upon horizontal separation δr. Top (bottom) panel: for Kn ~ 1 (Kn < 1) data. Conditioning and coloring same as for Fig. 9. 

Open with DEXTER  
In the text 
Fig. B.2.
Thirdorder moment of δu_{z} vs. δz, conditioned upon horizontal separation δr. Top (bottom) panel: for Kn ~ 1 (Kn < 1) data. Conditioning and coloring same as for Fig. 9. 

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

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

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.