Theory for the effect of fluid inertia on the orientation of a small particle settling in turbulence

Ice crystals settling through a turbulent cloud are rotated by turbulent velocity gradients. In the same way, turbulence affects the orientation of aggregates of organic matter settling in the ocean. In fact most solid particles encountered in Nature are not spherical, and their orientation affects their settling speed, as well as collision rates between particles. Therefore it is important to understand the distribution of orientations of non-spherical particles settling in turbulence. Here we study the angular dynamics of small prolate spheroids settling in homogeneous isotropic turbulence. We consider a limit of the problem where the fluid torque due to convective inertia dominates, so that rods settle essentially horizontally. Turbulence causes the orientation of the settling particles to fluctuate, and we calculate their orientation distribution for prolate spheroids with arbitrary aspect ratios for large settling number Sv (a dimensionless measure of the settling speed), assuming small Stokes number St (a dimensionless measure of particle inertia). This overdamped theory predicts that the orientation distribution is very narrow at large Sv, with a variance proportional to ${\rm Sv}^{-4}$. By considering the role of particle inertia, we analyse the limitations of the overdamped theory, and determine its range of applicability. Our predictions are in excellent agreement with numerical simulations of simplified models of turbulent flows. Finally we contrast our results with those of an alternative theory predicting that the orientation variance scales as ${\rm Sv}^{-2}$ at large Sv.


Introduction
The settling of particles in turbulence is important in a wide range of scientific problems. An example is the settling of small ice crystals in turbulence, a process that is considered in the context of rain formation from cold cumulus clouds [1,2,3]. Further examples are the settling of small aggregates of organic matter ('marine snow') [4], and the dynamics of swimming microorganisms [5,6,7] in the turbulent ocean.
The settling of spherical particles in turbulence has been intensively studied. Maxey and collaborators [8,9,10] found that turbulence increases the settling speed of small spherical particles. This pioneering work has led to many experimental and numerical studies, using direct numerical simulation (DNS) of turbulence, and it is a question of substantial current interest [11,12]. An important question is how frequently particles collide as they settle in turbulence [13,14]. The collision rate is influenced by spatial inhomogeneities in the particle-number density due to the effect of particle inertia. There is substantial recent progress in understanding this two-particle problem [15,16,17,18,19]. The conclusion is that settling may increase or decrease spatial clustering of spherical particles, and that it tends to decrease the relative velocities of nearby particles because settling reduces the frequencies of 'caustics', singularities in the inertial-particle dynamics [15].
Most solid particles encountered in Nature and in Engineering are not spherical, yet less is known about the settling of non-spherical particles in turbulence, and their settling depends in an essential way on their orientation. In a fluid at rest the orientation of a slowly settling non-spherical particle is determined by weak torques induced by the convective inertia of the fluid -set in motion by the moving particle. For a single, isolated particle in a quiescent fluid this effect is well understood [20,21,22,23]: convective fluid inertia due to slip between the particle and the fluid velocity causes non-spherical particles to settle with their broad side first. For axisymmetric rods, for example, symmetry dictates that the angular dynamics has two fixed orientations: either the rod is aligned with gravity (tip first) or perpendicular to gravity. At weak inertia, only the latter orientation is stable, so that the rod settles with its long edge first. But when there is turbulence, then turbulent vorticity and strain exert additional torques that cause fluctuations in the orientations of the settling crystals [1,24].
To understand the angular motion of a non-spherical particle settling in turbulence is in general a very complex problem, because there are many dimensionless parameters to consider. There is particle shape (shape parameter Λ), and the effect of particle inertia is measured by the Stokes number St. The importance of settling is determined by Sv, a dimensionless measure of the settling speed. The significance of fluid inertia is quantified by two Reynolds numbers, the particle Reynolds number Re p (convective inertia due to slip between particle and fluid velocity), and the shear Reynolds number Re s (convective inertia due to fluid-velocity gradients). The nature of the turbulent velocity fluctuations is determined by the Taylor-scale Reynolds number Re λ .
If the particles are so small that they just follow the flow and that any inertial corrections to the fluid torque are negligible (Re p = Re s = 0), then the angular dynamics of small crystals in turbulence is well understood [25,26,27,28,29,30,31,32,33,7,34]. The particle orientation responds to local vorticity and strain through Jeffery's equation [25]. The effect of particle inertia is straightforward to take into account [35], but the role of fluid inertia is more difficult to describe, even in the absence of settling. In certain limiting cases fluid-inertial effects are well understood. The most important example is that of a small neutrally buoyant (Re s = St) spheroid moving in a timeindependent linear shear flow, so that the centre-of-mass of the particle follows the flow (Re p = 0). Neglecting inertial effects (Re s = 0) and angular diffusion, the angular dynamics degenerates into a one-parameter family of marginally stable orbits, the socalled Jeffery orbits [25]. Fluid inertia breaks this degeneracy and gives rise to certain stable orbits [36,37,38,39]. Much less is known when Re p is not zero. Candelier, Mehlig & Magnaudet [40] recently showed how to compute the effect of a small slip upon the force and torque on a non-spherical particle in a general linear time-independent flow, by generalising Saffman's result [41,42] on the lift upon a small sphere in a shear flow, valid in the limit where Re p √ Re s 1. The results summarised in the previous paragraph pertain to time-independent flows. Time-dependent spatially inhomogeneous flows present new challenges, and very little is known about the effect of fluid inertia for such flows, in particular for turbulence. In some studies, therefore, effects of fluid inertia were simply neglected [43,44,45,46,47]. These models predict that the breaking of isotropy due to gravity causes a bias in the orientation distribution of the settling particles, so that rods tend to settle tip first, parallel to gravity. For small particles it is safe to neglect Re s [48]. But experiments and numerical simulations of slender particles settling in a vortex flow [49] and in turbulence [50] show that convective inertial torques due to settling can make a qualitative difference to the orientation distribution.
In this paper we therefore consider the effect of the convective inertial torques on the orientation of small spheroids settling in turbulence. Following Ref. [49], our model assumes that the hydrodynamic torque is approximately given by the sum of Jeffery's torque and the convective inertial torque in a homogeneous, time-independent flow. For nearly spherical particles this convective torque was calculated by Cox [20], and for slender bodies by Khayat & Cox [21]. Their results were generalised to spheroids with arbitrary aspect ratios in Ref. [22].
Our goal is to analyse how the turbulent-velocity fluctuations affect the orientation distribution of a prolate spheroid settling through turbulence. We assume that the particles are small enough so that convective-inertia effects due to the fluid-velocity gradients are negligible, that inertial effects on the centre-of-mass motion are small (small St and Re p ), but that the settling number Sv is large enough so that the fluidinertia torque dominates the angular dynamics.
We find an approximate theory for the angular distribution of settling spheroids using a statistical model [51,45] for the turbulent fluctuations. The theory is valid for large Sv and small St, in the overdamped limit, and its predictions are in excellent agreement with results of numerical simulations of the statistical model, and with simulations using a kinematic-simulation (KS) model [52,53] for the turbulent flow. We find that the variance of the orientation scales as Sv −4 in the limit of large settling number Sv, for small enough Stokes number St, and the theory determines how the prefactor depends on the shape of the spheroid. In the slender-body limit, the Sv −4 -scaling of the variance was also found in Ref. [54] using an approach equivalent to ours.
We contrast our results with a theory for the orientation variance derived by Klett [24] for nearly spherical particles. This theory predicts that the variance is proportional to Sv −2 . At first sight this may appear to be at variance with the overdamped theory, but we show that the overdamped approximation breaks down into several different regimes when particle inertia begins to matter. At very large values of Sv, when the time scale at which the fluid-velocity gradients decorrelate is the smallest time scale of the inertial dynamics, our numerical simulations show a Sv −2 -scaling, as suggested by Klett's theory. But the theory is difficult to justify because it neglects particle inertia in the centre-of-mass dynamics. We show that translational particle inertia has a significant effect upon the angular dynamics, so that it must be taken into account as soon as the overdamped approximation breaks down.
The remainder of this paper is organised as follows. In Section 2 we describe our model: the approximate equations of motion and the statistical model for the turbulentvelocity fluctuations. In Section 3 we show results of numerical simulations of our model. We describe how and why the results differ from those in Refs. [43,44,45,46,47], and explain the intuition behind our theory for small St and large Sv. The overdamped theory is described in Section 4. Section 5 discusses the effect of particle inertia, and Section 6 contains our conclusions as well as an outlook.

Particle equation of motion
Newton's equations of motion for a single non-spherical particle read: Here g is the gravitational acceleration with directionĝ = g/|g|, x p is the position of the particle, v p its centre-of-mass velocity, m p the particle mass, and the dots denote time derivatives. We assume that the particle is axisymmetric, so that its orientation is characterised by the unit vector n along the symmetry axis of the particle. The angular velocity of the particle is denoted by ω p , and I p (n) is its rotational inertia tensor per unit-mass in the lab frame. For a spheroid, the elements of I p (n) are given by [55] (I p ) ij (n) = I ⊥ (δ ij − n i n j ) + I n i n j , I ⊥ = 1 + λ 2 5 a 2 ⊥ , I = 2 5 a 2 ⊥ , (2) where λ ≡ a /a ⊥ is the aspect ratio of the spheroid, 2a is the length of the symmetry axis, and 2a ⊥ is the diameter of the spheroid. Prolate spheroids correspond to λ > 1, whereas oblate spheroids have λ < 1.
The difficulty lies in computing the hydrodynamic force f and torque τ on the particle. In the Stokes approximation, unsteady and convective inertial effects are neglected. In this creeping-flow limit [55], the force and torque upon the spheroid are linearly related to the slip velocity W ≡ v p − u, to the angular slip velocity ω p − Ω, and to the fluid strain S: Here µ is the dynamic viscosity of the fluid, u ≡ u(x p , t) is the undisturbed fluid velocity at the particle position x p , Ω ≡ 1 2 ∇ ∧ u is half the vorticity of the undisturbed fluidvelocity field at the particle position, and S is the strain-rate matrix, the symmetric part of the matrix of the undisturbed fluid-velocity gradients (its antisymmetric part is denoted by O). The tensors A, C, and H are translational and rotational resistance tensors. Their forms are determined by the shape of the particle. Eq. (3) shows that the tensor A relates the hydrodynamic force f (0) to the slip velocity W . For an axisymmetric particle with fore-aft symmetry the tensor takes the form The resistance coefficients A ⊥ and A depend on the shape of the particle. For a spheroid, they are given by [55]: For a sphere one has A ⊥ = A = 1, so that f (0) simplifies to the usual expression for Stokes force on a sphere moving with velocity v p through a fluid with velocity u. In the creeping-flow limit, the steady slip velocity W of a spheroid subject to a gravitational force m p g is obtained by setting the accelerationv p to zero in Eq. (1a): Here 1 is the unit matrix, and τ p ≡ (2a a ⊥ ρ p )/(9νρ f ) is the particle response time in Stokes' approximation with kinematic viscosity ν = µ/ρ f , fluid-mass density ρ f , and particle-mass density ρ p . The slip velocity depends on the orientation n of the particle.
For an axisymmetric particle with fore-aft symmetry, the rotational resistance tensors take the form: C ij ≡ C ⊥ (δ ij − n i n j ) + C n i n j and H ijk = H 0 ijl n k n l .
Here ijl is the anti-symmetric Levi-Civita tensor, and we use the Einstein summation convention: repeated indices are summed from 1 to 3. For a spheroid, the rotational resistance coefficients read [55]: Expressions (3) to (8) determine the hydrodynamic force and torque in the creeping-flow limit. Fluid-inertia effects are neglected in f (0) and τ (0) . There are two distinct corrections when fluid-inertia effects are weak but not negligible, due to the undisturbed fluid-velocity gradients, S and O, and due to the slip velocity W . The former are parameterised by the shear Reynolds number Re s , the latter by the particle Reynolds number Re p : Here a = max{a ⊥ , a } is the largest dimension of the particle, and W ⊥ is an estimate of the slip velocity: the magnitude of the velocity of a small slender spheroidal particle settling under gravity in a quiescent fluid with its symmetry axis perpendicular to gravity. From Eq. (6) we see that W (0) In the definition of Re s , the parameter s is a characteristic shear rate. In turbulence it is on average of the order s ∼ τ −1 K where τ K is the Kolmogorov time Here the average · · · is over Lagrangian fluid trajectories, and E is the turbulent dissipation rate per unit mass. This yields the estimate [48] Re s ∼ (a/η K ) 2 , where is the Kolmogorov length [56]. Thus the shear Reynolds number is small for small particles. Now consider the effect of convective inertia. Following Ref. [49] we assume that the torque on the particle is given by the sum of Jeffery's torque and the instantaneous convective-inertia torque in a homogeneous flow. This approximation can be strictly justified for a steady linear flow in the limit √ Re s Re p 1. In this limit the singular perturbation problem that determines the fluid-inertia torque simplifies: the Saffman length (∝ Re −1/2 s ) is much larger than the Oseen length (∝ Re p −1 ). This implies that the leading convective-inertial corrections to the torque are those corresponding to a quiescent fluid, and a similar argument can be made for the convective-inertia contribution to the force. While there is no general theory explaining how the convectiveinertia contributions to the force and the torque are affected by spatial inhomogeneities in time-dependent flows, the results of Ref. [49] show that the simple model used here can successfully explain important features of the orientation distribution of rods settling in a vortex flow.  25), as a function of the particle aspect ratio λ.
The leading-order inertial force correction reads for a spheroid moving in a quiescent fluid [57,21]: with W = |W | andŴ = W /W . For a spheroid, the leading-order inertial contribution to the torque was calculated in Ref. [22]: The shape factor F (λ) is given in Ref. [22]. It is also shown in Fig. 1(a). Combining Eqs. (1), (2), (3) with Eqs. (12,13) yields the equations of motion for our model. We use the Kolmogorov time τ K and the Kolmogorov length η K to dedimensionalise the equations of motion, This gives (after dropping the primes): Eqs. (14) have four independent dimensionless parameters: Here Λ is the shape parameter that appears in Jeffery's equation, and Sv is the settling number [58], a dimensionless measure of the settling speed. It is proportional to the particle size squared, a 2 , just as the Stokes number.
The shape-dependent prefactors in Eq. (14) are [in addition to those given by Eqs. (4) and (5)] as well as The Reynolds number Re p does not appear explicitly in Eqs. (14) because we dedimensionalised the equations of motion with the Kolmogorov scales τ K and η K . If we use an estimate of the slip velocity instead (such as W (0) ⊥ ), then Re p features in the dimensionless equations of motion. The latter convention is used in Refs. [21,22], and more generally in perturbative calculations of weak inertial effects on the motion of particles in simple flows [59,41,42,40]. These two different choices must lead to equivalent equations of motion, but our scheme has the advantage that it emphasises the different roles played by f (1) and τ (1) for small particles in turbulence. Eq. (14b) shows that the fluid-inertia contribution to the force, f (1) , is multiplied by the dimensionless prefactor a/η K . This means that f (1) makes only a small contribution for small enough particles, which we do not expect to qualitatively change the results derived below. In the following we therefore neglect this contribution (although it could be taken into account in simulations and theory).
More importantly, the fluid-inertia contribution to the torque in Eq. (14d) has no such factor. The fluid-inertia torque is of the same order as the Jeffery torque. This implies that the fluid-inertia contribution to the torque is potentially much more significant than the fluid-inertia contribution to the force. At large Sv in particular the particle settles rapidly, so that W is large. In this limit one therefore expects the fluid-inertia torque τ (1) to dominate over Jeffery's torque τ (0) , so that the inertial torque cannot be neglected (as was done in Refs. [43,44,45,46,47]). It is argued in Ref. [60] that the orientation bias predicted in Refs. [43,44,45] can possibly be observed in small-Re λ flow, but not at high Re λ .
In the following we neglect the contribution from f (1) . At the same time we assume that the settling speed is so large that the fluid-inertia torque τ (1) dominates the angular dynamics. If there was no flow, the particles would settle with their broad side first in this limit. The question is how turbulent fluctuations modify the orientation distribution of the settling particles.

Statistical model
In our theory we use a statistical model [51] to represent the turbulent fluctuations. We model the incompressible homogeneous and isotropic turbulent fluid-velocity field u(x, t) as a Gaussian random function with correlation length , correlation time τ , and rms magnitude u 0 (here and in Section 2.3 we write the equations in dimensional form because we want to make explicit how these scales are related to the Kolmogorov scales). Following Ref. [51] we express the fluid-velocity field u(x, t) in three spatial dimensions (3D) as The components A j of the vector field A are Gaussian random functions with mean zero, A j (x, t) = 0, and with correlation functions We choose the normalisation N 3 = 1/ √ 6 so that u 0 = |u| 2 . Below we also quote results for a two-dimensional (2D) version of this model. In this case we take with N 2 = 1/ √ 2, and where ∂ j represents the derivative with respect to the spatial coordinate x j . As the equation of motion for the 2D model we use Eq. (14) with n and the translational dynamics constrained to the flow plane.
The statistical model has an additional dimensionless parameter, the Kubo number [51] Ku = u 0 τ / . In the limit of large Ku the Gaussian random function u(x, t) models small-scale fluid-velocity fluctuations in the dissipative range of homogeneous isotropic turbulence. Evaluating Eq. (10) in the statistical model gives (Section 5.1 in Ref. [51]): where d is the spatial dimension. The spatial correlation length satisfies 2 = u 2 1 / (∂ 1 u 1 ) 2 , which defines the Taylor length scale [56] in turbulence. The length scale is related to the Kolmogorov length by [56,61] where C is a constant of order unity. The ratio /η K (or alternatively the Taylorscale Reynolds number Re λ ) constitutes a sixth dimensionless parameter of the model, in addition to the Kubo number and the four parameters listed in Eq. (15). In all statistical-model simulations described in this paper we set Ku = 10 and /η K = 10, and we determine the parameters τ and of the statistical model from Eqs. (21,22). The statistical model is constructed to approximate the dissipative-range fluctuations of 3D turbulence [51]. We note that the predictions of the 2D and 3D statistical models are essentially similar, but the two-dimensional model is easier to analyse, and it can be simulated more accurately. Two-dimensional and three-dimensional turbulence, by contrast, exhibit significantly different fluid-velocity fluctuations.

Kinematic-simulation model
To demonstrate the robustness of our theory we also compare its predictions to results of numerical simulations using a different model for the turbulent flow, namely the Kinematic-Simulation (KS) model [52]. The KS model has been shown to reproduce qualitatively many features of turbulent transport, and it provides a convenient way to represent a flow with a wide range of spatial scales, such as turbulence, albeit in a simplified manner. In short, we discretise Fourier space in geometrically spaced shells, up to a largest wavenumber. The largest and smallest length scales of the flow are L and η, respectively. The total number of shells is denoted by N k . We choose the characteristic wave vector in shell n as: k n = k 1 (L/η) (n−1)/(N k −1) . In each cell, we pick one wave vector, k n . The flow is then simply constructed as a sum of Fourier modes: The Fourier coefficients are chosen so that k n · a n = k n · b n = 0 (incompressibility), and with magnitude a 2 represents the Kolmogorov spectrum [56]. The frequency ω n in Eq. (23) is taken to be ω n = 1 2 . Further details about the implementation of this model for u(x, t) can be found in Ref. [53]. Figure 2 shows orientation distributions obtained by numerical simulations of Eqs. (14) for the three-dimensional statistical model described in Section 2.2. Shown are distributions of n g = n ·ĝ for a range of different Stokes numbers. We see that the particles settle with their broadside approximately aligned with gravity. For rods this means that n ⊥ĝ, so that n g = 0, and for disks n ĝ, so that n g = 1. These are the stable orientations for prolate and oblate particles settling in a quiescent fluid [21,22].

Orientation distributions
Compare the distributions in Fig.2 to those shown in Fig. 1 of Ref. [45]. There, by contrast, the rods tend to settle tip first, and disks tend to settle edge first. The reason for the difference is that the effect of the fluid-inertia torque was neglected in Ref. [45], whereas in the present work we choose parameters where this torque dominates the angular dynamics.
When the Stokes number is small we expect that the vector n spends most of its time close to a stable fixed point of the angular dynamics. But we expect that the turbulent velocity gradients modify the fixed point, so that it is no longer simply n g = 0 (rods) or n g = 1 (disks). Since the turbulent velocity gradients change as functions of time, the fixed-point orientation becomes time dependent too. In the overdamped limit (small Stokes numbers) we expect that the particle orientation follows the fixedpoint orientation quite closely. This allows us to derive a theory for the orientation distribution in this limit, described in the following Section.

Overdamped limit
The model (14) is very difficult to analyse in general. Therefore, to simplify the analysis, we consider a limit of the problem where the relaxation time of n is much faster than the time scale on which the gradients change as the particle moves through the flow. This corresponds to the overdamped limit of the problem, St → 0 in Eqs. (14). It was shown by experiments and numerical simulations in Ref. [49] that this limit quantitatively describes the orientation distribution of rods settling in a two-dimensional vortex flow, and in the slender-body limit this approach was also used in Refs. [54,62].
We also assume that Sv is large enough so that the fluid-inertia torque dominates the angular dynamics. This allows us to take into account turbulent fluctuations perturbatively. It also means that we can approximate the instantaneous slip velocity by W (0) (n), Eq. (6). In this limit we find: with n g = n ·ĝ, as defined in Section 3. The overdamped equation for the dynamics of the vector n corresponding to Eq. (24b) readṡ n = On + Λ[Sn − (n · Sn)n] + A Sv 2 n g (ĝ − n g n) . (24c) To simplify the notation we introduced the parameter Fig. 1(b) shows how A depends on the particle-aspect ratio λ.

Two-dimensional dynamics in the overdamped limit
We consider the 2D model first because it is much easier to analyse than the threedimensional model. We assume that the gravitational acceleration points into theê 1direction, and define φ to be the angle (0 ≤ φ < π) between n and this axis, so that n g = n ·ĝ = cos φ. For prolate particles (λ > 1 or equivalently Λ > 0) the overdamped angular dynamics (24c) becomes in two spatial dimensions: This two-dimensional overdamped equation of motion for the angular dynamics is essentially equivalent to model M2 in Ref. [49], used there for simulations of the angular dynamics of rods settling in a two-dimensional vortex flow. Apart from the fact that Ref. [49] considers a different flow, it describes small cylindrical particles with slightly different resistance tensors, and it approximates the n-dependence of the settling velocity. Equation (26) shows that the fluid-inertia torque has the same angular dependence as the S 11 -component of the strain, but in general the sign may differ. When S 11 > 0, the strain tends to align the rod withê 1 , the direction of gravity. The fluid-inertia torque acts against alignment with this direction. To quantify this statement, consider the fixed points of the angular dynamics (26). In the limit |A |Sv 2 → ∞ the inertial torque dominates the angular dynamics, so that the fluid-velocity gradients do not matter. In this limit the fixed points are φ * 1 = 0 and φ * 2 = π/2. For a prolate particle (λ > 1) φ * 1 = 0 is unstable while φ * 2 = π/2 is stable. This is the limit considered in Ref. [21], a slender rod falling in a quiescent fluid: since φ * 2 is stable the rod settles with its broad side first. For an oblate particle the stabilities are reversed [22].
What is the effect of the turbulent flow? In general this question is difficult to answer. But if the angle φ relaxes much more quickly than the fluid-velocity gradients change along the particle path, then the problem becomes tractable. Assuming that the gradients are constant, we can find exact expressions for the two fixed points of Eq. (26), for arbitrary aspect ratios and fluid-velocity gradients. We take λ > 1 and expand the stable fixed point around π/2 assuming that |A |Sv 2 is large: Here B ij are the elements of the matrix B = O + ΛS. Eq. (27) shows how the fixedpoint orientation changes as the turbulent velocity gradients evolve. We expect that the orientation of a settling rod follows these fixed-point orientations closely in the overdamped limit, provided that its angular relaxation time is smaller than the time scale on which the flow (and thus φ * 2 ) changes. We now analyse the angular dynamics of the settling particles in this 'persistent limit' [63]. Fig. 3 shows examples of how the fixed point φ * 2 (t) of the angular dynamics fluctuates as the particle settles through the turbulent flow and encounters different fluid-velocity gradients. The data are obtained by numerical simulation of the 2D model described in Section 2, for small Stokes numbers. Also shown is the instantaneous angle φ(t) obtained in these simulations. We see that the orientation dynamics follows the fixed point φ * 2 quite closely when St is small. In the overdamped limit the relaxation time τ φ of the angular dynamics (in units of τ K ) is given by the inverse of the stability exponent σ of the fixed point φ * 2 . From (26) we find to first order in (|A |Sv 2 ) −1 that σ ∼ −|A |Sv 2 . This gives When Sv is large, the fluid-velocity gradients seen by the settling particle change at the settling time scale τ s , the time it takes a particle settling with an angle φ = π/2 at a settling velocity given by Eq. (6) to fall one correlation length We therefore conclude that the persistent limit requires: This indicates that the persistent approximation works in the overdamped limit when |A |Sv is large enough. In the opposite limit, for small values of Sv, the settling time scale τ s is larger than the Lagrangian time scale, so that the fluid-velocity gradients change at the Lagrangian time scale, of order unity in units of τ K . Hence we must demand τ φ 1 to ensure that the persistent approximation works. This corresponds to the condition |A |Sv 2 1 .
In the persistent limit, the overdamped angular dynamics (26) responds so rapidly that the orientation of the particle follows the instantaneous fixed point of the dynamical system (26) quite closely. In this case the orientation distribution of the settling particle is determined by the distribution of φ * 2 , and thus by the distribution of fluid-velocity gradients encountered by the particle, through Eq. (27). This distribution may differ from the distribution of fluid-velocity gradients at a fixed spatial position (preferential sampling [51]). But in the overdamped limit preferential sampling of the fluid-velocity gradients is expected to be weak. We have checked that it is negligible for data shown in this paper. If we consider only the leading correction in Eq. (27), then the orientation distribution is determined by the distribution P B (B 12 ) of B 12 : In the two-dimensional statistical model the distribution P B (B 12 ) is Gaussian with variance σ 2 B = 1 8 (2 + Λ 2 ). This means that the distribution of φ is Gaussian too: with variance Eq. (32) shows that the distribution of φ simply reflects that of the fluid-velocity gradients, in the overdamped and persistent limit. The corresponding distribution of n g = n ·ĝ is:

Three-dimensional dynamics in the overdamped limit
In this Section we show how to obtain the distribution of n g = n ·ĝ for the threedimensional statistical model, in the same overdamped and persistent limit considered above. The calculation is analogous to the one described in Section 4.1. Let p = n−n gĝ . Using p 2 = 1 − n 2 g we express the equation of motion (24c) of n g aṡ n g =ĝ ·ṅ =ĝ · On + Λ[ĝ · Sn − (n · Sn)n g ] + A Sv 2 n g (1 − n 2 g ) = O gp +Λ[(1−2n 2 g )S gp + n g (1−n 2 g )S gg −n g S pp ] + A Sv 2 n g (1 − n 2 g ) . (36) Here the subscripts g and p denote contractions withĝ and p. In the limit of |A |Sv 2 → ∞, n * g = 0 is the stable fixed point for prolate particle (λ > 1). To determine how the fixed point changes due to fluid-velocity fluctuations we seek an expansion in (|A |Sv 2 ) −1 as in Section 4.1, of the form n * g ∝ 1/(|A |Sv 2 ) + . . .. We obtain to leading order: Assuming that the orientation of p is uncorrelated from the fluid-velocity gradients, we obtain for the variance of the distribution of n g : where σ 2 B is the variance of the distribution of B 12 (the gravitational acceleration points in theê 1 -direction). We also used that p 2 = 1 − n 2 g ≈ 1. This is a good approximation because in the limit we consider n g is small for prolate particles. Assuming that p and the fluid-velocity gradients are uncorrelated implies that the distribution of n g is Gaussian in the statistical model: and the variance evaluates to (40) Figure 5 shows results for the distribution of n g from simulations of the three-dimensional statistical model. Panel (a) shows results for small Stokes numbers, the parameters are the same as in Fig. 4(a). Also shown are the results of the theory, Eqs. (39) and (40). In this case St is small enough and Sv large enough so that the theory works very well. Panel (b) shows the orientation distribution for different Stokes numbers, to demonstrate how the theory fails when the Stokes number becomes larger. The behaviour is similar to that described in Section 4.1: the distribution widens as St increases. Eq. (38) says that the variance of the distribution of n g is inversely proportional to the fourth power of Sv, σ 2 ng ∝ Sv −4 , for large values of the settling number provided that the Stokes number is small enough. In Fig. 6(a) this prediction is compared with results of simulations of the three-dimensional statistical model. Shown is the variance of n g as a function of Sv, for two Stokes numbers. When the Stokes number is small we see that the prediction (40) works well for large Sv, as expected. Fig. 6(b) shows the kurtosis β 2 = n 4 g / n 2 g 2 , measuring the flatness of the distribution P (n g ). As predicted by the theory, the kurtosis approaches the Gaussian limit (β 2 = 3) for large settling numbers, at small enough Stokes numbers.
When Sv → 0 the variance tends to 1 3 and β 2 → 9 5 , indicating that the persistent approximation fails because Eq. (31) is no longer satisfied. In this limit the distribution of n g becomes uniform and independent of the Stokes number, because the angular dynamics is isotropic when gravitational settling is weak. Fig. 6(c) shows results for the variance from numerical simulations using the KS model (Section 2.3), for three different values of the Stokes number. The results are very similar to those obtained using the statistical model [ Fig. 6(a)]. There is good agreement with the overdamped theory, Eq. (38), at large Sv for small enough St. We determined σ 2 B from the KS simulations, so there are no fitting parameters in Fig. 6(c). The good agreement shows that the overdamped theory is robust, insensitive to the details of the spectrum of the velocity fluctuations. Fig. 6 also shows numerical data for larger values of St. For small Sv this makes little difference, the distribution is uniform. For larger Sv the numerical results first follow Eq. (38) or (40). But as Sv increases further, the overdamped theory starts to fail, the earlier the larger the Stokes number. This indicates that particle inertia begins to become important.

Effect of particle inertia
We saw in the previous Section that the overdamped theory breaks down at large Sv. To understand when and why the overdamped theory fails one must check the full inertial dynamics. We analyse the 2D statistical model first.

Two-dimensional model
Consider the angular dynamics in the absence of flow, to estimate the time scales that are important for the angular dynamics. When u = 0, the dynamics of the phase-space coordinate z ≡ (v px , v py , φ, ω) has the stable fixed point z * = (Sv/A ⊥ , 0, π/2, 0), gravity in the direction ofê 1 . The stability matrix follows from Eq. (14): where A was defined in Eq. (17). The relaxation time following from Eq. (41) is given by τ φ = max(−1/ σ i ), the maximal stability time of J. Here σ i denotes the real part of the i-th eigenvalue of J. One eigenvalue of this matrix is σ = −A ⊥ /St. We have computed the other eigenvalues numerically and analytically in limiting cases. We find that the time scale τ φ interpolates between Eq. (28) for small St and ∼St/A ⊥ for large St, for a fixed value of Sv. If we fix St, by contrast, then we find that the time scale τ φ interpolates between Eq. (28) for small Sv and ∼St/A ⊥ for large Sv. We expect that the overdamped approximation fails when the inertial estimate for the relaxation time of the angular dynamics, τ φ ∼ St/A ⊥ , becomes larger than the overdamped estimate Eq. (28). This means that the overdamped approximation requires Conversely, when Eq. (42) is not satisfied then particle inertia matters, so that the overdamped approximation must fail [ Fig. 6(a)]. For a quantitative comparison, Fig. 7(a) shows numerical results for the variance of the orientation distribution obtained from simulations of the two-dimensional model. We see that the overdamped approximation breaks down for values of Sv larger than ∼ A ⊥ /(|A |St), as predicted by Eq. (42). We observe that the variance decreases more slowly as Sv increases further. Fig. 7(a) also reveals that there is yet another, asymptotic regime at very large values of Sv -so large that it is difficult to achieve small Re p at the same time (Section 6). It is nevertheless of interest to analyse this regime, because it reveals the ingredients that a theory describing effects of particle inertia must contain. Fig. 7(a) suggests that for very large values of Sv. Our simulations indicate that the prefactor c 1 depends upon /η K , St, and upon λ (not shown). We surmise that this regime describes particles settling so rapidly that the settling time scale τ s is the smallest time scale in the system. This cannot hold unless τ φ ∼ St/A ⊥ is much larger than τ s , and this crossover occurs at

Sv St
We expect Eq. (43) to be accurate for values of Sv much larger than those given by Eq. (44). This condition is also shown in Fig. 7, and we see that the large-Sv regime starts at values of Sv approximately satisfying (44). Since condition (42) is violated in this regime, particle inertia must be taken into account. A difficulty is that particle inertia changes the translational as well as the angular dynamics. Thus it is no longer guaranteed that W = W (0) (n) (assumed in the overdamped theory of Section 4). This means that particle inertia is expected to modify the angular dynamics in at least two ways. Firstly, it introduces the time derivative d 2 dt 2 δφ into the angular dynamics. Secondly, the fluctuations of the torque change because W = W (0) (n) when particle inertia matters. This is discussed in Section 5.2. Fig. 7(b) shows how the variance of δφ depends on particle shape, for fixed Sv and St. There are four regimes. First, in the limit λ → ∞ the distribution is uniform and independent of the Stokes number. In this regime the dynamics is overdamped [condition (42)], but the persistent approximation fails because Eq. (31) is not satisfied. Second, for intermediate aspect ratios, both conditions are satisfied, so that the theory [Eqs. (33) and (34)] is accurate. Third, at λ becomes smaller, the overdamped approximation breaks down. In this regime particle inertia must be taken into account. Fourth, as λ → 1 the orientation distribution must become uniform. This cross-over happens very rapidly: for spheres (λ = 1) the orientation distribution is uniform, but already for λ ∼ 1.05 there is strong alignment.

Klett's small-angle expansion
Klett [24] proposed a theory for the orientation variance of nearly spherical particles settling in turbulence, including particle inertia in the angular dynamics. He uses that the orientation variance is very small for large values of Sv. This suggests to expand the equations of motion in small deviations of the angle φ = acos(n ·ĝ) from its equilibrium value: φ = φ * +δφ where φ * = π 2 for prolate particles. Klett assumes that W = W (0) (n) [Eq. (6)] and expands the angular dynamics for nearly spherical particles in δφ.
We can derive an equation of motion consistent with his by expanding Eqs. (14) to leading order in δφ, assuming that W = W (0) (n), and retaining only the leading terms in (|A |Sv 2 ) −1 . In this way we obtain for a prolate particle of arbitrary aspect ratio in three spatial dimensions: When we expand the geometrical coefficients in Eq. (45) for small Λ we find that the prefactors of the terms on the l.h.s. of this equation are almost identical, in this limit, to those in Eq. (17) of Ref. [24]. Slight discrepancies arise in the δφ-term because we use the expression for the inertial torque from Ref. [22], while Klett uses the form obtained by Cox [20] (the relative error of the prefactors is of the order of 10 −3 [22]). At any rate, Eq. (45) is simply a damped driven harmonic oscillator, with implicit solution Here Ω 0 = [C ⊥ /(2I ⊥ St)] 4|A |Sv 2 I ⊥ St/C ⊥ − 1. Note that we discarded terms related to the initial angle, because they cannot be important for the steady-state variance of δφ in the limit of large Sv. Squaring Eq. (46) and averaging over realisations of the turbulent fluctuations in the statistical model we obtain for large Sv where c 0 is a function of /η K , St, and of the aspect ratio λ. We neglected a Sv −3 contribution to δφ 2 because it is exponentially suppressed. Eq. (47) fails to describe the large-Sv behaviour (43), shown as the thick black dashed line in Fig. 7. This means that Eq. (45) cannot be used to estimate the large-Sv width of the orientation distribution, or to compute deviations from the overdamped theory. Which approximation causes Eq. (45) to fail? Since the variance is small for large Sv, δφ remains small at all times. Therefore we see no reason to doubt that the smallangle expansion is valid. This leads us to conclude that the assumption W = W (0) (n) breaks down, in agreement with our conclusions in the previous Section. To check this, we artificially imposed the constraint W = W (0) (n) in simulations of the twodimensional statistical model. The resulting large-Sv variance follows Eq. (47), and thus fails to give the correct scaling, Eq. (43). This demonstrates that it is important to allow W to deviate from W (0) (n) when particle inertia matters.
Klett's theory is difficult to justify from first principles because it assumes that W = W (0) (n). However, he obtains that δφ 2 ∝ Sv −2 , assuming that the fluid-velocity gradients on the r.h.s. of Eq. (45) are just white noise in time. In view of Eq. (44) it is possible that a first-principles theory may yield just that. But fluctuations of W − W (0) (n t ) yield additional time-dependent terms in the angular equation of motion that are expected to change the properties of the noise driving the angular dynamics, resulting in a different prediction for the orientation variance. More importantly, Fig. 7 demonstrates that δφ 2 ∝ Sv −2 applies only in the unphysical limit of very large Sv, and that particle inertia causes a complex parameter dependence of the orientation variance at smaller values of Sv, with a number of different regimes to consider.

Conclusions
Convective fluid inertia affects the orientation of a small axisymmetric particle settling in a turbulent flow. In Refs. [43,44,45,46,47] this effect was neglected. Here we considered a limit of the problem where it is dominant, but where turbulent fluctuations still matter. Our goal was to compute the distribution of orientations of a spheroid in turbulence, to work out how the torques due to convective fluid inertia and due to the turbulent velocity gradients affect the orientation distribution. In general the angular dynamics of the settling particle is very complicated. Here we looked at a limit in which the problem becomes tractable: we assumed small Stokes number (a dimensionless measure of particle inertia) and large settling number (dimensionless settling speed). For small Stokes numbers the dynamics is overdamped. For large values of the settling number, the angular dynamics becomes persistent: it relaxes much more rapidly than the fluid-velocity gradients change. In this limit the angular dynamics follows the fixed points determined by the instantaneous fluid-velocity gradients, and our theory for the orientation distribution relates the shape of the distribution to that of the instantaneous fluid-velocity gradients encountered by the settling particle. Our predictions are in excellent agreement with numerical statistical-model simulations, and with simulations using KS turbulence at large Sv and small enough St.
At large Sv the orientation distribution is very narrowly centered around the orientation the settling particle would assume in a quiescent fluid, in the absence of flow. The overdamped theory predicts that the variance of the distribution is proportional to Sv −4 for large Sv, and it determines how the prefactor depends on aspect ratio λ of the particle. In the limit λ → ∞ the variance was computed in Ref. [54].
We demonstrated that the overdamped theory breaks down at finite Stokes numbers, when the settling number exceeds a threshold determined by St. In this regime particle inertia matters. Klett [24] proposed a theory for the orientation variance for nearly spherical particles, taking into account particle inertia in the angular dynamics. His theory assumes that this dynamics is driven by the fluid-velocity gradients experienced by the settling particle, and that these gradients are uncorrelated in time so that diffusion approximations can be applied. Klett's theory predicts that the variance is proportional to Sv −2 , and we do observe this scaling for very large Sv, so large that the settling time is the smallest time scale of the inertial dynamics. But to derive a theory from first principles it is necessary to take into account particle inertia not only in the angular dynamics but also in the centre-of-mass motion, resulting in additional fluctuating terms in the angular equation of motion that are expected to change the orientation variance. More importantly, our simulations also show that particle inertia gives rise to a complex dependence of the orientation variance on particle shape, on the Stokes number, and upon the settling number. When the variance is small, it may be possible to derive a theory for the variance using small-angle approximations. But this remains a question for the future.
Here we applied our theory only to prolate particles. It is of interest to consider oblate particles too, because flat disks and slender rods have qualitatively different shape factors (Fig. 1). We therefore expect that the effect of particle inertia on the angular dynamics of flat disks can be quite different from that on slender rods. Also, we considered only the leading order in the inverse settling number, but the overdamped theory allows us to take into account higher-order corrections in this parameter. Such corrections change the relation between the fixed points of the angular dynamics and the fluid-velocity gradients experienced by the particle. This modifies the form of the distribution of n g , and it may explain the overshooting seen in Fig. 6(b) at moderate values of Sv, but the details remain to be worked out.
Here we analysed a limit of the problem where the fluid-inertia torque dominates the angular motion. In Refs. [43,44,45,46,47], by contrast, this torque was neglected. The question is thus whether one can find regions where inertial torque does not dominate. This is considered in Ref. [60]. The simulations described there show that the fluidinertia torque can be smaller than Jeffery's torque only when Re λ is small. In a very turbulent flow, when Re λ is large, the torque induced by fluid inertia is always dominant. More precisely, when the ratio of the correlation length over the Kolmogorov length is large, /η K ∝ Re λ 1/2 1, then the only possible orientation bias corresponds to nonspherical particles settling with their broad sides down, the limit considered here.
The experiments measuring the orientations of rods settling in a vortex flow described in Ref. [49] are performed in the overdamped limit. In the future we intend to apply the theory outlined in Section 4 to spheroids settling in a two-dimensional vortex flow, using the fact that the fixed points of the angular dynamics can be found explicitly as functions of the fluid-velocity gradients in two spatial dimensions. We will analyse the effect of particle shape by considering the angular dynamics of flat disks settling in such flows. Figure 1 indicates that the behaviour could be quite different from that of rods, because the shape factors are so different. This two-dimensional system is well suited to study the effects of finite Stokes numbers in more detail, because the two-dimensional dynamics is much simpler than the three-dimensional turbulent dynamics.
The overdamped theory [Eq. (38)] assumes that Sv is large, and that St is small enough. Since Sv = St gτ 2 K /η K = St/Fr, this requires some discussion. Here Fr = η K /(gτ 2 K ) is the Froude number [58]. We conclude that the Froude number must be small for the overdamped theory to work quantitatively. In turbulence Fr ∼ E 3/4 /(gν 1/4 ) where E is the dissipation rate per unit mass. Using ν ∼ 10 −5 m 2 s −1 and g = 10 ms −2 we find that Fr ranges from 0.002 at E = 1 cm 2 s −3 to 0.3 at E = 1000 cm 2 s −3 . So we require modest values of the dissipation rate per unit mass, E , for the theory to work quantitatively. This is the limit where gravity dominates over the turbulent fluctuations, the limit we intended to describe.
In the future it is necessary to address possible shortcomings of our model which approximates the inertial contributions to force and torque by those for a homogeneous steady flow. Even in the steady case it remains an open question how to model the torque when Re p and √ Re s are of the same order, even if both dimensionless numbers are small. Furthermore, turbulent flow is unsteady. While it is common practice to use steady approximations for the instantaneous force and torque (as we do here) it is not known how to compute contributions to the torque due to unsteadiness for general inhomogeneous flows. We expect that the methods presented in Ref. [40] can be generalised to treat at least spatially linear, unsteady flows. Finally, to justify our model for the inertial torque it is necessary that Re p is small. At the same time we assumed that Sv is large. From the definitions (9) and (15) of these dimensionless numbers we see that Re p = (a/η K )(Sv/A ⊥ ). To satisfy both requirements we must therefore assume the particles to be much smaller than the Kolmogorov length. Since η K ∼ (ν 3 /E ) 1/4 this condition is more easily met when E is small. In the slender-body limit, Khayat & Cox [21] obtained an improved approximation for the inertial torque, valid for larger Re p , which was tested in Ref. [49] and was found to agree better with the experiments at larger Re p . But corresponding corrections for other particle shapes are not yet known.