Inertia induces strong orientation fluctuations of non-spherical atmospheric particles

The orientation of non-spherical particles in the atmosphere, such as volcanic ash and ice crystals, influences their residence times, and the radiative properties of the atmosphere. Here, we demonstrate experimentally that the orientation of heavy submillimeter spheroids settling in still air exhibits decaying oscillations, whereas it relaxes monotonically in liquids. Theoretical analysis shows that these oscillations are due to particle inertia, caused by the large particle-fluid mass-density ratio. This effect must be accounted for to model solid particles in the atmosphere.

Introduction.-Thetransport, dispersion, and settling of volcanic ash [1,2], microplastic particles [3,4], and ice crystals in cold atmospheric clouds [5][6][7][8] have significant environmental impact.These non-spherical particles are subject to gravity, as well as viscous and inertial hydrodynamic forces and torques.An essential parameter characterising the latter is the particle Reynolds number, defined by Re p = av g /ν, where a is the linear size of the particle, v g is its settling speed, and ν is the kinematic viscosity of air.In general, the transport of non-spherical particles depends strongly on their angular dynamics [9][10][11], which directly affects the settling speed [1,6,12,13].This, in turn, determines its residence times and dispersion ranges in the atmosphere.The settling speed influences, for instance, how far microplastic, dust and volcanic ash can be transported away from a source, or how much time an ice crystal spends in a cloud [1][2][3]14].In addition, the angular dynamics determines the volume swept out.Together with the settling speed, this volume is a key parameter determining particle-particle collision rates [15], e.g.relevant for the formation of aggregates of ice particles in clouds [16][17][18][19] or volcanic ash [2].The particle orientation also has a direct impact on the absorption and scattering of radiation by the atmosphere, which affects albedo of atmospheric clouds [20][21][22], an effect still not understood quantitatively, despite its importance [23][24][25][26][27].
Many studies investigated the drag and stable orientation of non-spherical particles settling in viscous liquids at rest [13,[28][29][30][31][32], where the ratio R = ρ p /ρ f between the particle-mass density ρ p and the fluid-mass density ρ f is close to unity.When the particle Reynolds number Re p is larger than unity, the particle orientation aligns rapidly and monotonically, causing the particle to settle with its broad side down [11,33].The angular dynamics becomes unsteady at larger Re p .Settling particles in liquids, for example, exhibit a rich variety of motion patterns at Re p ∼ 100 [34].In air with R ∼ 1000, the angular dynamics of thin settling disks exhibits bistabil-ity [35] at Re p ∼ 50.One expects that particle inertia plays an important role in explaining this qualitative difference, since R is very large.
We remark that at much higher Re p , the angular dynamics of non-spherical particles settling in still air becomes unstable [1].Vortex shedding causes the characteristic fluttering first considered by Maxwell, see Ref. [36] and references cited there.
Experimentally, it remains a challenge to study the inertial angular dynamics of particles that settle very rapidly.One needs precise particle tracking over long periods of time, high-magnification imaging, a particlerelease mechanism that does not set the fluid in motion, and the container must be large enough to avoid spurious interaction with the wall [30].We developed a new experimental apparatus that overcomes these challenges [Fig.1(a)].Here we report on the first measurements, showing that the particle orientation exhibits characteristic oscillations with time scales comparable to those of atmospheric turbulence.The particles -with Reynolds numbers between 2 and 35 -do eventually align so that they settle with their broad side down.We explain that the oscillations are induced by particle inertia, and conclude that particle inertia can significantly enhance the extent to which turbulence randomises particle orientation.
Experiments.-Thesetup consists of an air-filled settling chamber (SC) with a novel particle injector, and two high-speed camera pairs (TX, TY) and (BX, BY), synchronised with a high-intensity pulsed LED array [Fig.1(a)].Each camera images a fall distance of ∼ 30 mm at a nominal resolution of 6.75 µm px −1 , The ap- Experimental setup.(a) Optical table with top cameras (TX and TY), and bottom cameras (BX and BY) named after the shown coordinate system (Z is the direction of gravity g), the settling chamber (SC), the synchronized pulsed LED unit (LED), controlled with a waveform generator (WG), and the oscilloscope (OSC).(b) Snapshots of a settling prolate spheroid recorded at 2932 frames per second.The tilt-angle -the angle between the particle symmetry axis and gravity -is shown in 5.1ms intervals in units of degrees.See Supplemental Material (SM) [45] for details.
paratus allows us to image individual solid 0.1 − 5 mmparticles settling in quiescent air.See Supplemental Material (SM) [45] for a complete description of the setup.
We used the Photonic Professional GT 3D printer (Nanoscribe GmbH) with sub-micrometer printing resolution to print submillimetre-sized spheroids with mass density ρ p = 1200 kg m −3 [46].Three size groups of spheroids were produced (Table I), with different aspect ratios λ = a ∥ /a ⊥ (2a ∥ is the length of the particle symmetry-axis, and 2a ⊥ is the perpendicular diameter).Three-dimensional scans of the particle shape using a 3D laser-scanning microscope (VK-X200K, Keyence)  I. Size groups of particles studied.Parameters: aspect ratio λ = a ∥ /a ⊥ (2a ∥ is the particle length along its symmetry axis, and 2a ⊥ is its diameter); volume Vp; Reynolds number Rep (using a = max{a ∥ , a ⊥ } and the observed settling speed); Stokes time τp = (2ρp/9ρ f )a ⊥ a ∥ /ν, where ρp and ρ f are the mass densities of the particle and the fluid, and ν = 1.5 × 10 −5 m 2 /s is the kinematic viscosity of air.show that the unevenness in the surface features of the spheroids is negligible (Fig. S2 in the SM [45]).In total, we carried out between 9 and 22 measurements per particle shape and size, resulting in a total of 170 successful experimental runs where the particle was in sharp focus for all four cameras.Fig. 1(b) shows recorded images of a prolate spheroid (2a ∥ = 410 µm, 2a ⊥ = 82 µm) as it falls in the settling chamber.Figs.S3-S4 [45] contain recorded images of all particles in Table I.Fig. 1(b) indicates that the orientation of the particle oscillates as it settles, in sharp contrast with previous experiments in liquids, where the alignment is monotonic [13,33,37,38,40].To explain the oscillations observed here, we developed a theoretical model that includes the effect of particle inertia.
Model.-Since the particle-to-fluid mass-density ratio R is large, one suspects that the observed oscillations in the angular dynamics are due to particle inertia.This expectation is consistent with a recent theoretical study [27] of tilt-angle fluctuations of small particles in turbulence.However, we cannot directly apply the model from Ref. [27], because that model was derived assuming Re p ≪ 1, whereas our particles have Re p of order 10.Therefore, we developed a theoretical description that extends the validity of the Re p ≪ 1model.Specifically, we modified the fluid-inertia contributions to the hydrodynamic force and torque by introducing two scalar functions C F and C T that depend on the settling speed, determined using ab-initio simulations of spheroids fixed in a uniform flow [40,[47][48][49][50][51].Details are given in Appendix A. As shown below, the new model accurately describes the experimental results, and highlights the key significance of particle inertia.The model has three non-dimensional parameters (Appendix C): the aspect ratio λ = a ∥ /a ⊥ , the non-dimensional particle vol-

FIG. 2. Comparison between experiments and theory. (a)
Time evolution of tilt-angle φ from experiment (blue) and its model prediction (red) for spheroids with aspect ratios λ = 0.2 and 5 from group I (Table I).(b) Steady-state settling speed v * g , frequency, and decay rate of the tilt-angle oscillations as functions of the aspect ratio λ.Markers show averages obtained for all experiments with error bars indicating 95% confidence bounds for groups I (red/•), II (green/□), and III (blue/⋄) (Table I).Solid lines show large-time asymptotes from a linear-stability analysis of the model described in Appendix A. Shaded regions indicate how much the theoretical predictions change as the measured settling speed varies along the particle trajectory.See Appendices A and B for details.Dashed lines show results of a linear-stability analysis of a harmonic-oscillator approximation, Eq. ( 6) in Appendix D .
ume, V = gV p /ν 2 , where g is the gravitational acceleration, V p = (4π/3)a 2 ⊥ a ∥ is the particle volume, and the mass-density ratio R = ρ p /ρ f .These parameters determine the steady-state settling speed v * g , and thus Re p .For small settling speeds, Re p ≈ (1/6π)RV , up to a λdependent factor of order unity.For oblate disks, particle inertia can be parameterised by a non-dimensional inertia ratio, J * .It is related to our non-dimensional parameters by J * = (π/64)Rλ (Appendix C).Our J * -values are at least two orders of magnitude larger than for disks settling in water [28,34,52], and larger than in [35] by a factor ∼ 2.
Results.-The experimental results are compared with model predictions in Fig. 2. Panel (a) shows how the angle between the particle-symmetry vector and gravity (the tilt-angle φ) decays to the steady value φ * = 0 for disks and φ * = π 2 for rods.The decay is oscillatory as opposed to the behavior in water [13,32] where the dynamics is overdamped (particle inertia is negligible) and the decay is monotonic.This is the case when the damping time τ ω of the angular velocity is smaller than the decay time τ φ ∼ ν/[v * g ] 2 of the tilt angle [33], which happens when R 3 V 2 ≪ 1 (Appendix C).In the opposite limit, τ ω ≫ τ φ .So particle inertia lengthens the transient.
Fig. 2(b) demonstrates that the model captures the observed settling dynamics very well.The largest disagreement is in the decay rate, which is hard to measure especially for nearly spherical particles.The white markers in the bottom panel show selected experiments where the decay was best fitted by an exponential.The smaller scatter of these data points, closer to the theory, suggests that systematic errors in extracting the data provide the most likely explanation of differences between theory and experiments.Fig. 2 reveals good agreement over the whole range of Re p and λ covered by our experiments.We mention that differences between model and experiment are expected to grow at larger Re p , because the determination of the inertial torque becomes less reliable beyond Re p ≈ 30 [51].
To develop a qualitative understanding of the oscillations, we simplified the model further, assuming that δφ = φ − φ * remains small and that the settling speed is large.In this limit we obtain a harmonic-oscillator equation for δφ, namely δ φ+δ φ+(V * g ) 2 C T |h(λ)|R 3 V 2 δφ = 0 in non-dimensional units (Appendix C).In particular, V * g is a non-dimensional settling speed ∝ v * g /(gτ p ), and h(λ) is a shape-dependent function, shown in Fig. S6.For Re p ≪ 1 and |δφ| ≪ 1, this equation simplifies to the form given in Refs.[27,33,53].Linear stability analysis of the harmonic-oscillator approximation shows that the particles approach alignment with an exponential decay, with decay rate The non-dimensional steady-state settling speed V * g is of order unity.For all particles in our experiments, the values of |h(λ)|R 3 V 2 were large enough to ensure that ∆ < 0, as shown in Fig. 3(a).A qualitative change occurs when ∆ becomes positive: then the particle orientation relaxes without oscillation.The expression for ∆ shows that this bifurcation cannot occur in the overdamped limit R 3 V 2 ≪ 1.We conclude: the bifurcation is due to particle inertia.Now consider the very slender fibres measured in Ref. [54].The asymptotic form of |h(λ)| for large λ im- ⊥ , disregarding factors of log λ.It follows that, for R ∼ 1000, only fibres with a ⊥ larger than ∼ 20µm can oscillate.This explains why the fibres with diameters ∼ 10µm used in Ref. [54], represented by the gray region in Fig. 3(a), did not oscillate.We conclude that the angular dynamics of slender fibres in the atmosphere can be very different from that of particles of moderate aspect ratios.In case of very slender disks, |h(λ)| ≈ 7 × 10 −5 λ, so the parameter combination |h(λ)|R 3 V 2 depends on particle geometry as (a ⊥ √ λ) 6 .Our estimates indicate that oscillations are observable for thin disks when a ⊥ √ λ is larger than ∼ 20µm.This condition is very well fulfilled for the oblate particles in Table I, as well as for a large class of ice crystals [6].A potential shortcoming of this analysis is that the forces and torques acting on thin disks have not been thoroughly tested for values significantly smaller than λ ∼ 0.1 [49,51].
In turbulence, tilt-angle fluctuations are determined by balancing the inertial relaxation described above with the effect of turbulent fluid-velocity gradients which upset alignment.Typical Kolmogorov times for weak atmospheric turbulence, τ K ∼ (ν/ε) 1/2 ∼ 0.4s to 0.04s for the dissipation rate per unit mass ε from 10 −4 m 2 /s 3 to 10 −2 m 2 /s 3 [6] are of the same order as the time scales in Fig. 2, indicating that particle inertia can significantly increase the randomising effect of turbulence, just as for small spherical particles, where particle inertia matters most for Stokes numbers of order unity and larger [15,[55][56][57].
To describe the effect of turbulence, we added a stochastic forcing [57] to our model (Appendix E). Figure 3(b) summarises the results.It shows the standard deviation σ φ of the tilt-angle fluctuations for small turbulent dissipation rate, obtained by simulations of the model from Appendix E. Also shown is the overdamped approximation obtained by neglecting translational and rotational accelerations (Eq.(S6) in the SM [45]).In Fig. 3(b), Re p was varied by changing the particle volume, keeping all other parameters the same.For small Re p , simulations and overdamped theory agree, so particle inertia has no effect.As Re p grows, σ φ decreases at first, because the fluid-inertia torque aligns the particle more strongly as the settling speed increases.At the same time, a difference between simulations and overdamped theory develops: particle inertia enhances the tilt-angle fluctuations.
At still larger Re p , σ φ starts to increase again, forming a characteristic minimum.This can be understood in terms of the harmonic-oscillator approximation: the minimum in Fig. 3(b) occurs at a critical value of Re p (details in SM [45]) where ∆ becomes negative, indicating that the increase is due to the bifurcation described above, causing transient oscillations that result in larger angular fluctuations.The bifurcation occurs when the Green's function of the harmonic-oscillator kernel becomes oscillatory, in the same way as without turbulence.The critical Re p where this happens is much smaller than the particle Reynolds numbers where bistability [35] I are shown as red/•, group II as green/□, group III as blue/⋄.Particles from Ref. [32] as +, and fibers from Ref. [54] as a gray region.We approximated cylindrical fibers as slender prolate spheroids and estimated ∆ by setting CT = CF = 1, since the corresponding Rep are very small.(b) Standard deviation of tilt-angle fluctuations for spheroids settling in weakly turbulent air with dissipation rate ε, as a function of Rep. Shown are simulation results of the model described in Appendix E in symbols, for spheroids with λ = 0.2 but different volumes, Vp.Solid lines correspond to the overdamped approximation, neglecting particle inertia (Eq.(S6) in SM [45]).The harmonic-oscillator bifurcation (∆ = 0) is shown as a dashed vertical line.
tering [34,36] is observed.We remark that a preliminary analysis suggests that our model may also explain the bistable angular dynamics observed in Ref [35].Also, dealing with more complex shapes is an important challenge.Ice crystals, for example, may be hollow or lack fore-aft symmetry [58].This may give rise to additional torques [13,59].
In summary, our experiments and model calculations show that particle inertia has a strong effect on the angular dynamics of atmospheric particles, generally enhancing the orientation fluctuations of settling atmospheric particles, not only in still but also in turbulent air.This causes increased settling velocities and lateral drift, in contrast to the drift-free pattern observed for steadily settling particles in liquids.Orientation fluctuations also affect the rate at which non-spherical particles collide [15] or fragment, a process important for secondary ice particle production [60,61].In addition, fluctuations in the orientation of ice crystals affect the radiative properties of ice-laden clouds, for example, by reducing cloud albedo when solar radiation is parallel to gravity.
Conclusions.-Weidentified the key importance of particle inertia for the motion of non-spherical atmospheric particles.Our results, made possible by the concurrent development of a unique experimental setup and by a reliable modelling strategy, show that heavy spheroids settling in air at Re p ∼ 1 − 30 -values typical for atmospheric particles -approach their stable orientation through decaying oscillations.We demonstrated that this behavior is a consequence of particle inertia.This physical effect must therefore be accounted for in models of important atmospheric processes, such as the radiative properties and evolution of ice-laden clouds, and residence times and dispersion ranges of volcanic ash or microplastics in the atmosphere.
TB  Appendix A Model for Re p up to ∼ 30 The dynamics of a settling particle is determined by Newton's equations for translation and rotation: Here x is the particle position, v its velocity, n is a unit vector parallel to the symmetry-axis of the particle, and ω is the angular velocity of the particle.Its mass is m p , and J p is the particle-inertia tensor [62].The gravitational acceleration is denoted by g.The main difficulty is to determine the hydrodynamic force F h and torque T h .For Re p ≪ 1, they can be determined in perturbation theory [9,11,63,64].For larger Re p -as in the experiment -one can parameterise forces and torques on a spheroid in uniform flow using ab-initio computer simulations [49,50].The conclusion is generally that force and torque can be parameterised by introducing empirically determined correction factors in the perturbative equations of motion.Here we use with correction factors C F (Re p , λ) and C T (Re p , λ).F and resistance tensors A(n, λ) and C(n, λ) (Eqs.(4,7) in [33]).The Re pcorrections are [9,11], with a = max(a ∥ , a ⊥ ).For spheroids, the shape factor F (λ) is given in Eqs.(4.1,4.2) of Ref. [64].For C F = C T = 1, Eq. ( 2) simplify to known expressions for Re p ≪ 1 [27].
Eq. ( 2) yields good results for F h up to Re p ∼ 100, for prolate particles with moderate λ [50], and there is similar qualitative agreement for oblate particles [49].We determined C F (Re p , λ) from interpolations of ab-initio simulation results for fixed Re p and λ for oblate [49] and prolate [50] spheroids as follows.The full model, Eqs.(1,2) has 11 dimensions.Since x is slaved to the other variables, it is sufficient to analyse the eight-dimensional system for v, ω, and n.In the experiments, the particles settled with speeds close to the steady-state settling speed.Therefore it suffices to evaluate C F at the steady state Rotational symmetry dictates that the polar angle θ * can take any value.Solving Eqs.(1,2) with a given (not yet known) value of C F , we find Here A (g) is the component of the translational resistance tensor A in the direction of gravity, for a particle falling with its steady-state orientation φ * .We used Eq. ( 4) to evaluate the Reynolds number, Re p = av * g (C F )/ν, with a = max{a ⊥ , a ∥ }.Since C F depends on Re p , and Re p in turn depends on C F through v * g (C F ), we solved the resulting implicit equation numerically to find the desired value of C F .The results are shown in Fig. 4(a).
We determined C T (Re p , λ) by interpolating the data in Refs.[49,50], using Re p determined as described above.The results are shown in Fig. 4(b), they are consistent with the ab-initio simulations from Ref. [51].
It is not guaranteed that the model works outside the tested parameter range, for example for very thin disks, or for very nearly spherical particles.Therefore we do not report numerical values for λ close to unity in Figs. 2  and 4. As mentioned above, we determined the functions C F (Re p , λ) and C T (Re p , λ) only near the steady state, for small δφ = φ − φ * (where φ is the tilt angle, and φ * = 0, π 2 is its steady-state value).This is sufficient as long as the dynamics does not depart too far from the steady state.Figs.2(a) and S5 show that this works very well.The small drift of the angular dynamics in Fig. 2(a) may be due to inaccuracies in C F or C T .The accuracy of the model could be improved by introducing correction matrices C F and C T in Eq. ( 2), instead of scalars, with elements that depend on φ, in addition to Re p and λ.

Appendix B Fitting the model to experimental data
This appendix contains the details needed to reproduce the theoretical fits in Fig. 2. First, Eq. ( 1) and ( 2) can be solved numerically for any initial condition v 0 , n 0 and ω 0 (eight parameters).To reduce the number of parameters, we fitted only the initial tilt-angle, φ, its angular velocity, φ, the initial settling speed v g , and the velocity component v ⊥ perpendicular to gravity.We assumed steady-state values for the remaining parameters.As a result, the dynamics resides in a plane determined by gravity and the direction of v ⊥ .The red lines in Fig. 2(a) and Fig. S5 in the SM [45] were obtained in this way.We see that the approximation works very well.
Second, to determine the shaded regions in Fig. 2(b), we perturbed the initial angular velocity and settling speed away from the above initial conditions, using typical experimental values for the particles from Table I.
Third, the solid lines in the top panel of Fig. 2(b) were obtained using Eq. ( 4).The solid lines in the middle and bottom panels of Fig. 2(b) were determined from linear stability analysis of the eight-dimensional dynamics of v, ω, and n.Linearising the dynamics around the fixed point (3), we obtained the eigenvalues describing exponential relaxation to the steady-state.Two eigenvalues form a complex pair.The real part gives the decay rate and the imaginary part gives 2π times the frequency.

Appendix C Non-dimensional parameters
We non-dimensionalise velocities with ṽ ≡ gτ p /A (g) , obtained from Eq. ( 4) in the limit of small settling speed, time with the angular-velocity relaxation time τ ω ≡τ p J ⊥ /(m p C ⊥ ), force with m p ṽ/τ ω , and torque with 5 a 2 ⊥ is the moment of inertia of a spheroid around an axis perpendicular to its symmetry axis, and C ⊥ is the rotational resistance coefficient around this axis [33].In particular, the steadystate settling speed is non-dimensionalised as V * g = v * g /ṽ.In these non-dimensional units, all terms in Eq. ( 1) are of order unity, except for ⊥ a ∥ is the volume of the particle.So there are two non-dimensional parameters in addition to λ = a ∥ /a ⊥ : the mass-density ratio R = ρ p /ρ f , and the non-dimensional particle volume V = gV p /ν 2 .In the experiment, R= 996, and V ≈ 0.06 (group I in Table I), V ≈ 0.1 (group II) and V ≈ 1 (group III).In the limit of small settling speeds, the parameters V and R are connected to Re p by Re p ∼ 1 6π RV , up to a λ-dependent prefactor.The overdamped limit of the angular dynamics is obtained when τ ω ≪ τ φ , where τ φ ∼ ν/[v * g ] 2 is the relaxation time of the tilt angle in this limit.So the overdamped limit corresponds to Willmarth et al. [28] quantified particle inertia for settling disks by the phenomenological parameter ⊥ is the moment of inertia of a short cylinder around its large axis, and 5 is proportional to the moment of inertia of a fluid sphere with diameter 2a ⊥ .For oblate spheroids, this expression reduces to J * ∝ λR.

Appendix D Harmonic-oscillator approximation
The planar dynamics described in Appendix B can be further simplified if δφ is so small that its feedback upon the settling speed v g can be neglected, and if the transversal velocity is much smaller than the settling velocity.Then the tilt-angle obeys a damped-pendulum equation.In the dimensionless units introduced in Appendix C it reads: The first term on the r.h.s. of the equation for ω is the rotational damping due to particle inertia.The time scale is chosen so that the prefactor is unity.The dependence of the function h(λ) on λ is shown in Fig. S6.For thin disks, |h(λ)| ∝ λ.In this limit, the prefactor in the last term on the r.h.s. of Eq. ( 5) evaluates to |h(λ)|R 3 V 2 ∝ J * Re 2 p .This rationalises the use of J * to describe the effect of fluid inertia on settling disks [28].
Although Eq. ( 5) is approximate, it captures the main features of the experimentally observed dynamics.The approach to the steady-state is conveniently analysed by linear stability analysis.Linearisation of Eq. ( 5) yields the harmonic-oscillator approximation In the limit of Re p ≪ 1, a corresponding equation was considered earlier, see for example Eq. ( 45) in Ref. [33].
Eq. ( 6) implies that the tilt-angle fluctuation δφ relaxes to zero exponentially, δφ ∼ a + exp(µ + t) + a − exp(µ − t), with eigenvalues , and purely imaginary otherwise.In the former case, relaxation towards the steady-state is monotonic.In the latter case, it involves oscillations.We remark that |h(λ)| takes small values [Fig.S6].This is compensated by R 3 V 2 ≳ 10 5 for the particles in Table I, with R ∼ 10 3 .For these particles, ∆ is negative.For much smaller ratios of particle-to-fluid mass densities, the condition ∆ < 0 is harder to satisfy.
The results of the harmonic-oscillator analysis are shown in Fig. 2 as dashed lines.In the top and middle panels they are indistinguishable from the solid lines.For the exponential decay rate, by contrast, the harmonicoscillator approximation differs from the solid lines.The bifurcation predicted by the harmonic-oscillator analysis is shown in Fig. 3 (dashed lines).

Appendix E Effect of turbulence
The effect of turbulence can be modeled by adding a stochastic forcing to the model described in Appendix A, representing the turbulent fluid velocity by a Gaussian random function u(x, t) (Eq.( 5) in Ref. [57]).Since C F and C T in Eq. ( 2) are approximated assuming small tilt angles (Appendix A), we must assume that the turbulent dissipation rate is sufficiently small.The particle equations of motion (1,2) change in the presence of the turbulent flow u(x, t).In the expressions for F (0,1) h and T (1) h , the particle velocity v is replaced by the slip velocity v − u(x, t) at the particle position x.The second change is that there is an additional torque due to the gradients of the imposed flow [65]: ] is half the turbulent vorticity, and S(x, t) is the strain-rate matrix, the symmetric part of the matrix of fluid-velocity gradients.The tensors C and H depend on particle shape.For spheroids, they are given in Eq. ( 7) in Ref. [33].The colon symbol represents a double contraction of indices [10].
here, the maximum cross-sectional area is only ∼ 0.3 mm 2 , which is less than 0.004% of the SC cross-section.In other words, the particles were far enough from the walls of the container to neglect additional torques due to wall interactions [2].The SC has four glass windows, a detachable top cover with a square opening to mount the PI, and a detachable bottom drawer to collect the dropped particles.With the exception of the opening in the removable top cover, the remaining five edges of the SC are sealed airtight to minimise potential airflow in the SC caused by the motion of the PI, as well as airflow in the laboratory.The high-precision XY Z-stage is attached beneath the SC, and allows to move the SC with 10 µm spatial resolution in all three directions.This resolution is a crucial requirement for precise camera calibration.

I.1.2. Particle injector (PI)
Since the aim of this study was to investigate the unsteady settling motion of single non-spherical particles, it was necessary to release the particles individually, with sufficient control over their translational and angular velocities.In particular, one must control the initial direction of motion so that the particle stays in the focus of all cameras.It was also important to ensure that the particles detach from the needle when desired, and that the initial release speed is neither too small nor too large, so that steady-state motion is reached in the field of view of the cameras.Finally, the flow disturbance in the SC due to particle injection should be minimised.We designed a particle injection mechanism that meets these requirements, as much as possible: it releases single particles of different sizes with minimal flow disturbance, and with a high success rate in imaging the particles with all cameras.Fig. 1(a) in the main text and Fig. S1 present the design of the PI.
The experiment was performed by taking a particle from the glass substrate on which it was printed, and placing the particle on the tip of the PI needle (Fig. S1).The needle is a cylindrical rod with a length of 200 mm, a diameter of 12 mm, and with a conical tip with a end diameter of 1 mm or 1.5 mm for attaching the particles.The needle was inserted into a hollow needle guide, and locked in place with a release key inserted into one of the needle grooves (Fig. S1).When the release key was removed, the needle dropped vertically until stopped by a needle block, where the sudden impulse caused the particle to detach.Depending on which groove was used, the needle reached a maximal speed of 0.42 − 1.5 m s −1 , and the particle therefore reached a similar initial speed upon detachment (the locking and eventual release of the needle from higher grooves, at smaller distances between the needle and the needle block, resulted in lower initial particle speeds).Upon release, the cameras were activated by an external photoelectric trigger placed in front of the needle block.
The PI was designed so that the particles had negligible lateral velocities, resulting in a horizontal displacement angle of at most 0.3 • .The initial angular velocity was harder to control, because it is more sensitive to how the particle is attached to the needle tip, as well as the detailed detachment mechanisms.

I.1.3. Cameras
Four Phantom VEO4K 990L cameras (Vision Research Inc.) were used to observe the SC.The cameras were positioned at different heights (Z).The optical axes of the top cameras (TX and TY) were in the X-and Y -directions, and the same for the bottom cameras (BX and BY), see Fig. 1(a) in the main text.Only the TX camera received direct light from the LED, the camera BX received light reflected by a mirror.The light from the LED was reflected through the SC in the Y -direction using a mirror of 204 × 254 mm 2 surface area at a 45 • -angle facing the LED.The TY camera received light from this mirror, while the BY camera was illuminated after a 90 • reflection by another mirror.
All cameras used a Nikon ED AF Micro Nikkor 200 mm 1 : 4D lens with a focal length of 200 mm configured at F-stop 22, and a magnification of 1 : 1, which corresponds to a image-pixel-size of 6.75 µm px −1 .The effective exposure time was 2.5 µs.For the range of measured settling speeds, 0.2 − 1.5 m s −1 , the pixel shift (defined as the exposure time multiplied by the particle velocity, essentially the image streak) was 0.08 − 0.56 px during the effective exposure time.The cameras were synchronised as follows.The camera TX was the leader camera.Its strobe signal was synced to the time codes of the three remaining cameras.The Phantom camera control (PCC) software was used to configure the role of the cameras and the trigger functions.
With a fixed data throughput of 9 Gpx s −1 [Phantom VEO-Series Datasheet, Vision Research, AMETEK (2021)], increasing the spatial resolution comes at the cost of temporal resolution.We used the maximum possible resolution (4096 px) in the Z-direction (the observation volume measured 27.648 mm in Z, at 1 : 1 magnification, corresponding to images with 6.75 µm px −1 ).The horizontal extent of the observation volume was covered with either 520 px (3.51 mm) at 4000 frames per second (for spheroids with λ = [0.8,1] in group I, and for all spheroids in groups II and III in Table 1 in the main text, or with 720 px (4.86 mm) at 2932 frames per second (for spheroids with large aspect ratios λ = [0.2,0.5, 1.25, 2, 5] in group I).Therefore, the horizontal extent of the observation volume was only four times the largest dimension of the λ = 4 particle from Group III, which indicates how small the observation volume had to be in order to successfully track the particle.The lateral displacement of the particle increases when the aspect ratio is much smaller or larger than unity, which affects the probability of observing the particle in focus in all cameras, i.e. the experiment success rate.The success rate was highest for almost spherical spheroids at around 34 % and lowest for spheroids with the largest aspect ratios at around 14 %.

I.1.4. Illumination
The LED contains 72 white LEDs, summing to an average power consumption of 300 W. It is operated in the 'pulseoverdrive' mode which can emit light at a pulse rate of 0 − 20 kHz [Product Manual, LED-Flashlight 300 white / blue / blue, IP54, LaVision GmbH (2021)].The brightness of this LED was larger than 1 000 000 lx at a distance of 1 m in the pulse-overdrive mode, and the duty cycle was between 0 − 10%.
A waveform generator (WG) was used to pulse the LED.The WG also controlled the pulse rate, amplitude, the offset of the waveform, and the duty cycle of the LED.The strobe signal from the TX leader camera was fed into the WG, causing it to generate a synchronisation signal for exposure times of all cameras.An oscilloscope (OSC) was used to monitor the strobe signals of the cameras, indicating the camera exposure times and phases.

I.1.5. Image analysis
The setup was calibrated immediately before or after each day of experiment.This involved taking images of a transparent and planar calibration target covered with 47 250 black circular dots simultaneously by all cameras, while the target was placed inside the SC.The circular dots on the calibration target had diameter 250 µm, printed with a feature tolerance of less than 0.8 µm, and a dimensional tolerance over the whole calibration target of less than 6 µm.The calibration target was placed inside the SC at an angle of 45 • and moved by the XY Z-stage in the XY -plane (see Fig. S1) at 0.5 mm increments up to 1 mm on either side of the focus plane.
Each camera was first independently calibrated with a camera model [3], to minimise the difference between the world coordinates given by the target and their projection recorded during the calibration procedure.Then, the camera parameters for each camera pair were further optimised by minimising the error between the mean of the three-dimensional world coordinates of the calibration target viewed by both cameras, and the coordinates obtained by stereo triangulation.In all cases, the root-mean-square error in estimating the X, Y and Z world coordinates of dots on the calibration targets was smaller than a pixel, and in most cases below one-third of a pixel.
For each camera, a median image background was calculated using all available frames.The background was then removed from the images.A small region around the particle in each frame was then cropped, and the resulting image was independently thresholded.The centroid position, the longest length, and a fitted ellipse were extracted from the thresholded image using the OpenCV 2 library in Python.Particle centroids were used to determine the three-dimensional positions of the particles using the stereo-optimised camera model.The particle velocities were then extracted using a sixth-order finite-difference scheme.For each spheroid image in a given frame, an ellipse was fitted (shape and orientation).For prolate particles, the angle between the major axis of the fitted ellipse and the gravity defines the tilt-angle.

I.1.6. Sensitivity analysis
Temperature sensitivity.To ensure that the air column inside the SC was not affected by the illumination, the change in air temperature in the SC must remain low during operation of the LED.We monitored the temperature inside the SC with a Humicap ® HMP7 humidity-and-temperature probe from Vaisala.The temperature sensitivity of the SC was tested in two LED illumination settings: 1) 15 s of illumination following a 105 s of no illumination, and 2) continuous illumination for a total duration of 30 min, at the highest LED illumination setting of 200 000 lx for both the cases.In the first setting, similar to the way we conducted the experiments, the temperature rise was 0.1 • C during the illumination time and dropped close to the ambient temperature during the non-illumination time, resulting in a total rise of 0.2 • C after 30 minutes.In the second setting, a temperature increase of 1.5 • C was observed after 30 min.The change in air mass density and viscosity due to the 1.5 • C-temperature rise is negligible, so no buoyancy-induced flow should affect the settling particle.
Impact of particle injector on the flow inside the SC and setup calibration.The movement of the needle before particle injection may result in undesired air currents in the SC.To quantify this airflow, oil-based smoke was filled into the SC.We released the PI needle from different locking grooves corresponding to different particle injection speeds, and tracked the smoke movement inside the SC.We found that the smoke remained stationary when the PI needle was released from most heights, except the highest ones.When the needle was positioned at the highest available height (second-highest height), it reached a velocity of 1.5 m s −1 (1.32 m s −1 ) before hitting the needle block.In this case, smoke movement was recorded with a vertical velocity of about 0.07 m s −1 (0.01 m s −1 ) near the tip of the needle.
A detailed study of the settling speed of spheres of different volumes V p = 1.44×10 −3 , 1.70×10 −3 , 17.97×10 −3 mm 3 was carried out to verify the setup.Table S1 compares the observed steady-state settling speed to the empirical model of Clift & Gauvin [4].The relative error is smaller than 5%.We note that the stated relative error of the empirical model is about 6% at Reynolds numbers below 3 × 10 5 [5].

I.2. Printed particles
The particles were printed using a 'Photonic Professional GT' 3D printer from NanoScribe.This printer polymerises the photoresist liquid with two-photon polymerization (2PP) technique, using 780 nm photons from a 3D Galvo laser scanning system [Photonic Professional (GT) User Manual, Nanoscribe GmbH (2015)].Two printer objectives are available in this setup.The 63×objective was used for high-resolution small-feature printing with an IP-dip 2PP resin on a fused silica substrate with a 200 × 200 µm 2 horizontal print field.The standard hatching distance (horizontal width of each polymerised dot) is 0.2 µm, while the standard slicing distance (vertical thickness of each polymerised dot) is 0.3 µm.The 25×objective was used for medium-resolution printing with an IP-S 2PP resin on an ITO(Indium Tin Oxide)-coated substrate with a 400 × 400 µm 2 horizontal print field.The standard hatching distance for this case was 0.5 µm, while the standard slicing distance was 1 µm.
Once the particles were printed onto the substrate, the un-polymerised portion of the 2PP resin was removed, which is necessary for the particles to be developed.In this development process, the substrate was immersed together with the particles and the remaining 2PP resin in a SU-8 developer, consisting mainly of 1-methoxy-2-propanol acetate.This process took about 30 min, when the un-polymerised resin dissolved in the SU-8.The SU-8 developer on the substrate and printed particles was then removed by dipping it in isopropanol alcohol for 2 min.
In order to reproduce the increasingly curved outlines of the slender particles, they were printed using the 25×objective in horizontal slices of 0.6 µm thickness.The printer polymerises the photoresist liquid IP-S to a specific voxel size (of 0.3 µm thickness in our case) using two photons on a conductive ITO-coated substrate.When the largest axis of the spheroid exceeds 280 µm, the particle was printed in two separate halves, which were then stitched together layer by layer during a later stage of the printing process.The quality of the resulting particle surface was checked with a Keyence VK-X200K laser microscope, see Fig. S2.The resulting unevenness of the particles is in the sub-µm range.We note that it is important to control the particle shape very precisely, because small shape irregularities can result in additional hydrodynamic torques that change the angular dynamics.Weak breaking of fore-aft symmetry, for example, changes the steady alignment angle [6,7].How it affects the transient is not known.

I.3. Supplemental experimental figures
Figs. S3-S4 show the recorded images of one experiment for each particle type listed in Table 1 in the main text.Fig. S5 shows in blue solid lines the time series of tilt-angle and settling speed post-processed from the experimental data of Figs.S3 and S4, together with the model predictions in red solid lines.The noise seen in the experimental time series for the settling speed is due to experimental uncertainties of the particle position in the sub-pixel range, caused by mirror distortion, out-of focus imaging, and uneven backlighting between adjacent LED lamps.In turbulence, fluid-velocity gradients give rise to an additional torque, Jeffery's torque [8] in the Stokes approximation.Refs.[9,10] used this model to describe the angular dynamics of small spheroids settling in a turbulent flow.They accounted for the inertial torque T (1) h , because its contribution is small in the limit of small Re p .The results indicate that particle inertia -the acceleration terms in Eq. ( 1) in Appendix A -can have a significant effect on the settling dynamics.This model can at best give qualitative insight into the angular dynamics of atmospheric particles, because it was derived for Re p ≪ 1, while atmospheric particles tend to have much larger particle Reynolds numbers.
Starting from the improved model described in Appendix A, the effect of weak turbulence can be modeled as described in Appendix E, by adding a stochastic forcing.To this end, one represents the turbulent fluid velocity by a Gaussian random function u(x, t) with root-mean square velocity u f , correlation time τ f , and correlation length ℓ f [11].The particle equations of motion -Eqs.( 1) and (2) in Appendix A -change in the presence of an imposed external flow u(x, t).In the expressions for F (0,1) h and T (1) h , the particle velocity v is replaced by the slip velocity v − u(x, t) at the particle position x.The second change is that there is an additional torque due to the gradients of the imposed flow, computed first by Jeffery [8].In dimensional form it reads: Here Ω(x, t) = 1 2 [∇ × u(x, t)] is 1 2 times the vorticity of the undisturbed flow, and S(x, t) is the strain-rate matrix, the symmetric part of the matrix of fluid-velocity gradients of the undisturbed flow.The tensors C and H depend on particle shape.For spheroids, they are given in Eq. ( 7) in Ref. [9].The double contraction represented by the colon symbol is explained in Ref. [12].Fig. 3(b) in the main text shows results for the r.m.s. of the tilt-angle fluctuations in weak turbulence that were obtained by simulations of the resulting model (symbols).
Also shown in Fig. 3(b) is the r.m.s.tilt-angle in the overdamped limit (solid lines).In this limit, the relaxation time of the angular velocity is much smaller than that of the tilt angle, τ ω ≪ τ φ .Here τ ω = τ p J ⊥ /(m p C ⊥ ) and τ φ ∼ν/[v * g ] 2 [9].In these expressions, v * g is the steady-state settling velocity, ν is the kinematic viscosity of air, m p is the mass of the spheroidal particle, C ⊥ is the component of C perpendicular to the symmetry axis, and J ⊥ is the moment of inertia around an axis perpendicular to the symmetry axis.
The overdamped limit is obtained by neglecting translational and angular accelerations.This corresponds to solving for particle orientation n, velocity v, and angular velocity ω by setting force and torque in Eq. ( 1) in the main text to zero.In a quiescent fluid, the solution is given by the stable fixed point in Eq. ( 3) of the main text, and the tilt angle δφ vanishes at this fixed point.
In turbulence, by contrast, turbulent flow-gradients modify the steady-state orientation, resulting in a small, but non-zero fixed point δφ * that depends on the local flow gradients [9,10].As a consequence, components transversal to gravity of the steady-state settling velocity obtain non-zero values.The steady-state angular velocity, however, remains zero [9,10].An expression for the tilt-angle variance in turbulence in the overdamped approximation was derived in Ref. [10] (see also references cited there) in the limit of small particle Reynolds number, Re p ≪ 1, and disregarding fluid-inertia correction to translation (the term F (1) h in Eq. ( 2) in Appendix A).Now we explain how to generalise this calculation to include the fluid-inertia correction to translation, building on the results of Ref. [10].
When the tilt angle remains small, its angular dynamics is governed by Eq. ( 18) in Ref. [10].This dynamics is driven by both the flow gradients and the settling through a matrix Y with components where W = v − u is the slip velocity, O and S are the antisymmetric and symmetric parts of the flow-gradient matrix, and Λ = (λ 2 − 1)/(λ 2 + 1) is a shape factor that depends on the particle aspect ratio λ.Moreover, |A | = Aligning gravity with the e 1 direction, the overdamped angular dynamics of a disk-like particle approaches a non-zero stable fixed point (Eq.(S6) in Ref. [10]) By solving Eq. (1a) of the main text with zero acceleration, with the orientation given by Eq. (S2), gives the transversal components of the steady-state settling velocity of disk-like particles.In the limit of large enough settling velocity, so that δφ is small, and the vertical component W g is much larger than the transversal components as well as the fluid gradients times the Kolmogorov length, we have for k = 2, 3, and where v * g is given by Eq. ( 4) in the main text.Using these expressions in Eq. (S1) allows to evaluate δφ * in Eq. (S2) The steady-state tilt-angle for rod-like particles can be derived in a similar manner, following the analysis in Ref. [10].
Averaging over a Gaussian distribution of fluid gradients gives Here A (g) = A ∥ and A (p) = A ⊥ for disk-like particles (λ < 1) and vice versa for rod-like particles.Further, ε is the turbulent dissipation rate per unit mass.From Eq. (S6), one obtains Eq. ( 16) in Ref. [10] by letting C F = 0 and C T = 1.In this limit, the steady-state settling velocity in a quiescent flow, Eq. ( 4) in Appendix A in the main text, simplifies to v * g = gτ p /A (g) .In Fig. 3(b) in the main text, the results obtained from Eq. (S6) are shown as solid lines.The tilt-angle variance σ 2 φ given by Eq. (S6) depends linearly on the energy-dissipation rate per unit mass ε.The dependence upon the kinematic viscosity ν, by contrast, is non-linear because the viscosity is also contained in the steady-state settling speed v * g , and in the particle Reynolds number Re p = v * g a/ν.The parameters ε and ν must be mapped to the statistical-model parameters.As described in Ref. [11], we match the Kolmogorov time ν/ε to the time scale of the flow gradients, (2⟨tr( SS T )⟩) −1/2 , and the Kolmogorov length (ν 3 /ε) 1/4 to ∼ 0.1 times the correlation length ℓ f of the random function u(x, t).
Also shown in Fig. 3(b) is the critical Reynolds number Re c where the discriminant ∆ in the harmonic-oscillator approximation changes sign (see appendix D in the main text).From the condition ∆ = 0, we obtain Re c = max(λ, 1) 12πA (g) C T |h(λ)|R . (S7) The function h(λ) is shown in Fig. S6, and R = ρ p /ρ f is the mass-density ratio.Eq.S7 is an implicit condition because C T depends on Re c .This condition was solved to obtain the value in Fig. 3(b), but the scale of Re c can be obtained immediately by using C T ∼ 1 for the particles considered in this paper.

II.2. Supplemental theory figure
Fig. S6 shows a plot of the function h(λ) featuring in the harmonic-oscillator equation, Eq. ( 5) in Appendix D in the main text. is the observed settling velocity of the sphere together with the standard error.vg(m) is its settling velocity from Clift & Gauvin (1971) empirical model [4].Error correspond to the relative error between vg(m) and vg(e).

FIG. 1 .
FIG. 1.Experimental setup.(a) Optical table with top cameras (TX and TY), and bottom cameras (BX and BY) named after the shown coordinate system (Z is the direction of gravity g), the settling chamber (SC), the synchronized pulsed LED unit (LED), controlled with a waveform generator (WG), and the oscilloscope (OSC).(b) Snapshots of a settling prolate spheroid recorded at 2932 frames per second.The tilt-angle -the angle between the particle symmetry axis and gravity -is shown in 5.1ms intervals in units of degrees.See Supplemental Material (SM)[45] for details.

FIG. 3 .
FIG. 3. (a) Bifurcation diagram.The dashed horizontal line distinguishes decay without oscillations (∆ > 0) from decay with oscillations (∆ < 0).Particles from group I in TableIare shown as red/•, group II as green/□, group III as blue/⋄.Particles from Ref.[32] as +, and fibers from Ref.[54] as a gray region.We approximated cylindrical fibers as slender prolate spheroids and estimated ∆ by setting CT = CF = 1, since the corresponding Rep are very small.(b) Standard deviation of tilt-angle fluctuations for spheroids settling in weakly turbulent air with dissipation rate ε, as a function of Rep. Shown are simulation results of the model described in Appendix E in symbols, for spheroids with λ = 0.2 but different volumes, Vp.Solid lines correspond to the overdamped approximation, neglecting particle inertia (Eq.(S6) in SM[45]).The harmonic-oscillator bifurcation (∆ = 0) is shown as a dashed vertical line.

FIG. 4 .
FIG. 4. Empirical coefficients CF and CT in Eq. (2) for the parameters in Fig. 2, obtained as described in Appendix A. (a) Force coefficient CF as a function of aspect ratio λ.The groups refer to Table I.(b) Torque coefficient CT .

1 .
Model for the effect of weak turbulent fluctuations

3 λ 2
+1 |F (λ)|J ⊥ /(m p A ∥ A ⊥ C ⊥ ), where A ∥ and A ⊥ are the components of A parallel and perpendicular to the symmetry axis of the spheroidal particle.The function F (λ) and the coefficients C F and C T define the torque model.They are described in Appendix A in the main text.

12π 2 125 ( 2 + 2 × 2
λ 2 + 2λ 4 ) A (p) C ⊥ m p A (g) C T F (λ)J ⊥ 1 + 3 16 C F Re p (3A (p) − A (g) ) 1 + 3 8 C F Re p A (g) if λ < 1 (disk-like) 1 if λ > 1 (rod-like) .(S6) FIG. S1.(a) Photo of the particle injector (PI) attached to the top cover of the settling chamber The cameras in the X-direction, TX and BX can be seen, as well as the mirror that reflects light through the SC in the Y -direction.(b) CAD drawing of the particle-injector (PI) components, that are installed on the top of the settling chamber (SC).(c) Schematic view of the setup showing the mirror (M) arrangements and illumination/imaging paths.

TABLE
[41,54]ded by the German Research Foundation (DFG) Walter Benjamin Position (project no.463393443).JG was supported by funding from the European Union Horizon 2020 Research and Innovation Programme under the Marie Sklodowska-Curie Actions, Grant Agreement No.675675.KG was supported by a grant from Vetenskapsrådet (no.2018-03974).BM was supported by Vetenskapsrådet (grant no.2021-4452), and acknowledges a Mary Shepard B. Upson Visiting Professorship with the Sibley School of Mechanical and Aerospace Engineering at Cornell.Statistical-model simulations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC).This research was also supported in part by the National Science Foundation under Grant No. NSF PHY-1748958.The Göttingen turret was manufactured with the support from the Mechanical Workshop and Electronic Workshop of the Max Planck Institute for Dynamics and Self-Organisation.We thank J.L. Pierson for pointing out Refs.[41,54]to us.

TABLE S1 .
Settling speed of spheres, based on total number of experiments Nexp for a sphere with volume Vp, diameter D, and Reynolds number Rep of the particle.vg(e)