Angular dynamics of a small particle in turbulence

We compute the angular dynamics of a neutrally buoyant nearly spherical particle immersed in an unsteady fluid. We assume that the particle is small, that its translational slip velocity is negligible, and that unsteady and convective inertia are small perturbations. We derive an approximation for the torque on the particle that determines the first inertial corrections to Jeffery's equation. These corrections arise as a consequence of local vortex stretching, and can be substantial in turbulence where local vortex stretching is strong and closely linked to the irreversibility of turbulence.

The angular dynamics of small non-spherical particles in flows is often discussed in terms of Jeffery's theory [1][2][3][4], neglecting the effects of particle and fluid inertia. It is assumed that the instantaneous torque on the particle vanishes at every instant in time, so that the angular velocity of a small spherical particle equals half of the fluid vorticity at the particle position. But this is no longer true when the particle is so large that inertial effects become important. It is straightforward to take into account particle inertia [3,[5][6][7]. This may be a good approximation for heavy particles, but for neutrally buoyant particles the acceleration of the surrounding fluid cannot be neglected. How to model the effect of fluid inertia on the angular motion of a suspended particle is an important question [4,7], yet difficult to answer. Past theoretical studies have therefore often concentrated on simple cases, for example on the angular dynamics of axisymmetric particles in shear flows, analysing the stability of Jeffery orbits under inertial perturbations [8][9][10][11].
In this Letter we estimate the first inertial contributions to the angular dynamics of a nearly spherical particle in a time-dependent and spatially varying flow. Our three most important assumptions are that the particle is small, almost neutrally buoyant, and nearly spherical (details are given below). We find that a spherical particle rotates at an angular velocity different from the local flow vorticity. The difference is caused by local vortex stretching. For nearly spherical particles we find additional contributions that depend in a more complex manner on the local fluid-velocity gradients. In turbulence these effects can alter the angular particle dynamics substantially. Recently the question was raised how the dynamics of tracer particles reflects the irreversibility of turbulence [12][13][14][15]. Our results show how the inertial angular dynamics of a small particles is linked to the breaking of time-reversal invariance in turbulence.
Formulation of the problem. Consider a small, nearly neutrally buoyant, spheroidal particle (with rotational symmetry axis n) in a fluid. We assume that the particle is nearly spherical, its aspect ratio λ is close to unity: λ ≡ 1 + where is a small parameter determining the eccentricity of the particle: > 0 for prolate particles, while < 0 for oblate particles. We assume that the particle is small: κ ≡ 2a/ 1. Here 2a is the length of the symmetry axis of the particle, and is the length over which the flow can be linearised near the particle.
We assume that inertial effects matter, but that they are weak. Convective inertia due to the fluid-velocity gradients is characterised by the shear Reynolds number Re s ≡ a 2 s/ν, where ν is the kinematic viscosity of the fluid, and s measures the magnitudes of the fluid-velocity gradients, the strain rate. The effect of unsteady inertia depends on the time scale τ c that describes how fast the boundary conditions change. The ratio of the magnitude of unsteady and convective inertia defines the Strouhal number Sl ≡ (sτ c ) −1 . The effect of particle inertia on the angular dynamics is determined by the Stokes number, the ratio of the rate of change of angular momentum and the torque: St ≡ (ρ p a 5 s)/(τ c µsa 3 ) = (ρ p /ρ f )Re s Sl. Here ρ f and ρ p are the mass densities of the fluid and the particle, and µ = ρ f ν is the dynamic viscosity.
We treat the effect of inertia perturbatively, this requires Re s and Re s Sl to be small (but not too small, see below). We disregard the effect of translational slip. This is justified for small particle Reynolds number, Re p ≡ av s /ν Re 1/2 s . Here v s is the slip velocity. To ensure that it is small enough we assume that the particle is approximately neutrally buoyant.
Equations of motion. The equations that govern the angular particle dynamics read: Here I is the moment-of-inertia tensor of the particle, with elements I ij = A I n i n j + B I (δ ij − n i n j ), A I and B I are the moments of inertia around and transverse to the axis n, and T is the hydrodynamic torque: The integral is over the particle surface S p , ds is the outward normal surface element, and σ(x, t) is the stress tensor of the fluid at position x, and r ≡ x−x p where x p is the particle position. The torque (2) is determined by the solution of Navier-Stokes equations. We decompose the stress as σ = σ ∞ + σ (1) , where σ ∞ is the stress tensor of the undisturbed Eulerian fluid velocity, denoted by U ∞ (x, t). The second term, σ (1) , is the contribution to the stress tensor from the disturbance flow. The torque is decomposed in a similar way, T = T ∞ + T (1) . We compute these two contributions to the torque separately.
Torque due to stress of undisturbed flow. The undisturbed Eulerian fluid velocity U ∞ (x, t) satisfies Navier-Stokes equations in the laboratory frame: The last equality defines the Lagrangian derivative along U ∞ , ∇ denotes the spatial derivative with respect to x, and the partial time derivative is evaluated at fixed x. The torque due to the undisturbed stress can be expressed as a volume integral using Eq. (3): To evaluate this expression further we expand DU ∞ (x, t)/Dt around the particle position x p : The components are relative to a fixed Cartesian basis in the laboratory frame. Substituting this expansion into Eq. (4) we find: The order is O(r 4 / 4 ) because terms odd in r vanish upon integration. The partial derivative in the integrand of Eq. (6) evaluates to ∇( where the elements of A ∞ p are the gradients of the undisturbed fluid velocity U ∞ at the particle centre: Finally, to perform the integral in Eq. (6) we use the definition of the moment-of-inertia tensor. We obtain: Here we have decomposed the matrix A ∞ p into its symmetric and antisymmetric parts, the strain-rate matrix (7) is valid for a spheroid with arbitrary aspect ratio (provided that κ 1).
FIG. 1. Schematic of vectors used to describe the particle motion, x is a position in the laboratory frame, O is its origin, and xp(t) is the particle position at time t. The perturbation theory leading to Eq. (13) assumes that the particle follows the Lagrangian fluid trajectory.
Using the vorticity equation DΩ ∞ Dt −S ∞ ·Ω ∞ = ν∇ 2 Ω ∞ evaluated at the particle position, we can express the difference in the first row of Eq. (7) in terms of the Laplacian of Ω ∞ which vanishes for a strictly linear flow. In general, however, the non-linearity of U ∞ (x, t) results in a non-zero value of ∇ 2 Ω ∞ at the particle position.
Disturbance torque. We calculate the torque due to the disturbed fluid in perturbation theory, assuming that inertial effects are small and neglecting translational slip. This dictates how the problem must be de-dimensionalised: t = τ c t , r = ar , U = saU , and σ = µsσ . In the remainder of this Letter we use these dimensionless variables. To simplify the notation we drop the primes, all equations below are written in dimensionless variables. The disturbance caused by the particle has the flow velocity w(r, t) ≡ U (r+x p , t)−U ∞ (r+x p , t). It is defined to be a function of r(t) ≡ x − x p (t) (Fig. 1). In dimensionless variables the disturbance problem reads: The partial time derivative is evaluated at fixed r, and we have linearised When is it justified to use this linear form in the disturbance problem (8)? The disturbance caused by the particle decays exponentially at distances larger than the Saffman length S ≡ 1/ √ Re s , so that we must require S < , where is the length scale over which the flow can be linearised. This condition is more restrictive than κ 1. In other words the shear Reynolds number Re s should not be too small, because convective inertia causes the disturbance to decay at S .
We use the reciprocal theorem [16] to find the hydrodynamic torque on the particle, given an 'auxiliary' Stokes solution in the same geometry [9,10,17,18]. In dimensionless variables the reciprocal theorem reads: The first term on the r.h.s. is the torque due to the undisturbed fluid stresses, Eq. (7), expressed in dimensionless variables. The second term is Jeffery's torque [1]: where K and H are Brenner's resistance tensors [19]. The third term in Eq. (9) is the torque due to the disturbance flow beyond the Stokes approximation. The tensor M(r) is determined by the known auxiliary Stokes solution (ũ = M ·ω), the Stokes flow around the particle rotating with angular velocitỹ ω in a quiescent fluid, see supplemental material [20]. The integration domain in (9) is the fluid volume outside the particle, and Re s f is defined as the r.h.s of the first Equation (8).
To simplify the calculation of the third term in Eq. (9) we assume that the particle is nearly spherical [18], λ = 1 + , and expand in the small parameter . To order 2 , for example, the moments of inertia around and transverse to the particle symmetry axis [18] read in dimensionless variables: ). For small values of the contribution of the volume integral in Eq. (9) can be evaluated. Details are given in the supplemental material [20].
For a spherical particle ( = 0) we find to order O(Re s ): We now use Eqs. (10) and (1) to determine the particle angular velocity ω to order O(Re s ). For a spherical particle I has the elements I ij = 8πδ ij /15, in dimensionless variables. In perturbation theory dω/dt = dΩ ∞ p /dt = DΩ ∞ /Dt| xp . We must also require that ρ p ≈ ρ f , so that Re p remains small. To order O(Re s ) we find: Eq. (11) is the main result of this Letter. We see that a very small particle in a viscous flow rotates with half the fluid vorticity, Ω ∞ p , as expected. The first inertial correction term, proportional to DΩ ∞ /Dt| xp , resembles the form of the slip velocity of a small particle subject to particle inertia and to the force due to the undisturbed pressure gradients. A small-St expansion gives for Sl = 1. Maxey used this approximation to conclude that heavy particles are centrifuged out of vortical regions in turbulence [22,23]. The expression for v s is of the same form as the second term on the r.h.s. of Eq. (11), since Re s = (ρ f /ρ p )St for Sl = 1. This term predicts that a particle that is slightly heavier than the fluid rotates a little bit more slowly, because it cannot keep up with the the fluid acceleration. A lighter particle, by contrast, rotates faster than Ω ∞ p . But these arguments disregard the S ∞ p ·Ω ∞ p -term in Eq. (11). This inertial term connects the angular particle dynamics to vortex stretching in the undisturbed flow. This term vanishes in steady linear  Eq. (13), along a Lagrangian path in turbulence (DNS at Re λ = 433 using JHU turbulence database [26,27]). b Vorticity along the same path. c Distribution of Y (DNS, lines). Also shown: experimental data (•) at Re λ = 50 read off from Fig. 2d in Ref. [28], rescaled by τK ≡ Tr A T · A −1/2 , see text. The experimental value (0.27 s −1 ) was taken from Ref. [29].
flows such as a simple shear [24,25], and for planar flows because vorticity is orthogonal to the flow plane.
Angular velocity and vortex stretching in turbulence. Fully developed turbulent flows exhibit large local vortex stretching rates, this is how the turbulent kinetic energy is dissipated at small scales. Irreversibility of turbulence [12][13][14][15] implies that the stretching rate Ω ∞ · (S ∞ · Ω ∞ ) does not average to zero [28,30]. This matrix element is positive on average because Ω ∞ tends to align with the middle eigenvector of S ∞ [12], and its eigenvalue is positive on average. Eq. (11) shows that the inertial correction to the angular velocity of a neutrally buoyant sphere in turbulence is determined by the stretching rate, and thus linked to the breaking of time-reversal invariance in turbulence.
But are the conditions of validity summarised above met for a neutrally buoyant particle in turbulence? The shear rate, defined by s = Tr S · S 1/2 , is of the same order as τ −1 K , the inverse of the Kolmogorov time τ K ≡ Tr A T · A −1/2 . This means that the shear Reynolds number is of order Re s ∼ a 2 /η 2 K , where η K = √ ντ K is the Kolmogorov length. We conclude that the particle must be smaller than the Kolmogorov length for the perturbation theory in Re s to be valid. For a nearly neutrally buoyant particle we can disregard gravity, and we assume that no other external forces act on the particle. Therefore the timescale τ c is that of the fluid and we take τ c = s −1 ∼ τ K , so that the Strouhal number is unity. The perturbation theory also requires that the Saffman length S ≡ a/ √ Re s ∼ η K is smaller than the length over which the flow can be linearised. In turbulence ∼ 10η K . So this condition is only marginally satisfied. Finally, translational slip is negligible when the Oseen length is much larger than the Saffman length. This is the case when Re 2 p Re s . For small neutrally buoyant particles in turbulence this condition is well satisfied.
We evaluate the inertial correction in Eq. (13) for a neutrally buoyant particle numerically, following a Lagrangian trajectory in fully developed turbulence. We use the JHU turbulence database [26,27] that contains a direct numerical simulation of forced, isotropic turbulence at Re λ = 433. The result is shown in Fig. 2a. Panel b shows the vorticity along the same path. We observe the stretching of a vortex tube at t ≈ 20 τ K . We see that the inertial correction to the particle angular velocity can be substantial during vortex stretching, a factor of order unity times Re s . Panel c shows that the distribution of the inertial correction has heavy, non-Gaussian tails that give rise to large values of the inertial correction. Comparison with experimental data at Re λ = 50 (•), from Fig. 2d in Ref. [28], shows that the tails are quite robust even at moderate values of Re s , a consequence of the universality of dissipative-range turbulent fluctuations [31].
Non-spherical particles. We have computed the angular velocity also for nearly spherical particles (| | 1). Neglecting inertial effects we obtain an -expansion of Jeffery's equation [1]: For Sl = 1 we find using the method described above (details in the supplemental material [20]): Eq. (14) shows that the inertial corrections to the angular velocity of non-spherical particles depend intricately on the relative alignment of the particle symmetry axis, of the vorticity, and of the eigensystem of the strain-rate matrix [2,3,32]. For a neutrally buoyant particle in incompressible isotropic homogeneous turbulence we can average the correction. Since δω contains a factor , we obtain the O( Re s )-result by averaging n uniformly over the unit sphere, n i n j = δ ij /3. This gives: Ω ∞ p · δω = (4/3) Re s Ω ∞ p · (S ∞ · Ω ∞ p ) . On average the form of the correction is similar to that in Eq. (11). We see that the inertial effect is weakened for slightly prolate neutrally buoyant particles, but increased for oblate particles.
Conclusions. We have computed the first inertial corrections to the angular velocity of a small, approximately neutrally buoyant particle in a space-and timedependent flow. Our main prediction, Eq. (13), expresses the inertial correction to the angular velocity of a small sphere in terms of a matrix element that determines vortex stretching. This shows that the inertial angular dynamics of a small neutrally buoyant sphere in turbulence picks up that time-reversal invariance is broken [12][13][14][15].
Our results demonstrate that convective and unsteady fluid inertia must be treated on an equal footing in turbulence. It is commonly argued that convective effects can be disregarded when the translational slip velocity v s is negligible. But we have shown here that substantial inertial corrections may arise from convective terms due to turbulent strains. Such terms are likely to be important in the translational problem too, so that the 'Maxey-Riley' equations [33,34] cannot be used to describe the translational dynamics of small spheres in turbulence.
Particle-tracking experiments [35][36][37] and particleresolving DNS required to test the predictions in this Letter have recently become possible [38,39], but they are still very challenging. In order to stimulate the substantial effort required for these measurements we now give a concrete suggestion for how Eq. (13) could be tested in either experiment or particle-resolving DNS. First, the particle should be neutrally buoyant, and its size a must be smaller than the Kolmogorov length η K . Small a ensures that Re 2 p Re s since Re p ∝ a 3 [Eq. (12)] while Re s ∝ a 2 . This may be a difficult technical requirement. It is easier to meet when the turbulence intensity is low, so that η K is larger, and Fig. 2c shows that the effect persists for lower turbulence intensities. Second, to determine the vorticity of the undisturbed flow, one must measure and interpolate the flow near the particle [35]. Then we suggest to consider the distribution of |ω| 2 /|Ω ∞ p | 2 − 1. The width of this distribution must approach zero for a perfect tracer particle. Fig. 2c shows that the first effect of inertia is to substantially widen the tails of this distribution, and to slightly shift its mean value. Finally, our results for nearly spherical particles indicate that disks may be more sensitive to inertial corrections than rods.
Can our results be generalised to cases where v s is not negligible? An important case is settling [40][41][42][43]. The settling of ice crystals, for instance, is important for rain initiation in cold turbulent clouds [44]. For a spatially constant flow the effect of non-zero v s was analysed by Lovalenti & Brady [16]. Can one use their methods to compute lift forces [45] on small particles in turbulence? This is a difficult problem because it requires singular perturbation theory [46]. Finally, larger particles pose other problems: they sense inertial-range turbulent fluctuations [47], and wakes may affect their dynamics [48]. sions and for making the preprint [49] available to us that considers a related question (the steady-state limit of the problem considered here). This work was supported by Vetenskapsrådet [grant number 2013-3992], Formas [grant number 2014-585], and by the grant 'Bottlenecks for particle growth in turbulent aerosols' from the Knut and Alice Wallenberg Foundation, Dnr. KAW 2014.0048. The numerical results in Fig. 2 use data from the JHU turbulence database [26,27].
All equations in this supplemental material are in dimensionless variables, t = τ c t , r = ar , U = saU , and σ = µsσ . For details see the Letter. For notational convenience we drop the primes in the following.

I. EVALUATION OF JEFFERY'S TORQUE
Eq. (10) in the Letter gives an expression for Jeffery's torque [S1]: written in terms of Brenner's resistance tensors K and H [S2].
Here ω s ≡ ω − Ω ∞ p is the angular slip velocity, Ω ∞ p = 1 2 ∇ ∧ U ∞ | xp , and S ∞ p is the symmetric part of the matrix A ∞ p of the undisturbed fluid-velocity gradients evaluated at x p . The anti-symmetric part is denoted by O ∞ p . This matrix is related to Ω ∞ p by O ∞ p · r = Ω ∞ p ∧ r. Expanding the resistance tensors in the nearly-spherical limit one finds [S3]: Here is a small parameter related to the particle eccentricity (see Ref. [S7]). For a sphere we have K ij = 8πδ ij and H = 0.

II. COMPUTATION OF DISTURBANCE TORQUE USING RECIPROCAL THEOREM
The hydrodynamic torque is given by The integral is over the particle surface S p , and the vector ds points in the direction of the outward surface normal. The stress tensor σ is determined by the solution of the Navier-Stokes equations. In dimensionless variables it has elements σ ij = −δ ij p + 2S ij where p is the pressure and S ij are the elements of S ∞ . We decompose σ as σ = σ ∞ + σ (1) , where σ ∞ is the undisturbed stress tensor and σ (1) is the contribution to the stress tensor from the disturbance flow. The hydrodynamic torque is split up in the same way: We use the reciprocal theorem [S4] to compute T (1) , following Refs. [S5-S7]. The reciprocal theorem relates integrals of two different Newtonian and incompressible flows: Here w is the disturbance flow, it obeys Eq. (9) in the Letter. The fields u and σ are the fluid velocity and the stress tensor of an 'auxiliary' Stokes problem in the same geometry, namely the Stokes flow produced by particle rotating with angular velocity ω in a quiescent fluid: u = ω ∧ r for r ∈ S p and u → 0 as r → ∞ .
The hydrodynamic (auxiliary) torque on the particle is denoted by T : The surface integrals in Eq. (S5) are over the particle surface, the vector ds points in the direction of the outward surface normal. The volume integrals are over the volume outside the particle. The surface integrals in Eq. (S5) are evaluated by substituting the boundary conditions for the disturbance flow and the auxiliary flow on the particle surface. The volume integral on the l.h.s. of Eq. (S5) vanishes by virtue of Eq. (S6a).
For the divergence of the disturbance stress tensor in the volume integral on the r.h.s., we substitute Eq. (9) in the Letter. The vector f is defined by We find: The surface integrals can be expressed in terms of the torques T (1) and T : This is an exact relation between the sought disturbance torque and the known auxiliary torque. The disturbance velocity w, however, is not known. To lowest order in Re s we can use the Stokes approximation w (0) in evaluating the volume integral. The function w (0) is the solution of the Stokes problem: w (0) = ω − Ω ∞ p ∧ r − S ∞ p · r for r ∈ S p and w (0) → 0 as r → ∞ . (S11b) The first two terms on the r.h.s. of (S10) evaluate to ω · T (0) . Adding T ∞ [Eq. (S12)] we obtain Eq. (10) in the Letter, with u = M · ω.

III. EVALUATION OF TORQUE AND ANGULAR VELOCITY
A. Sphere in a steady flow We first consider a simpler problem than that addressed in the Letter. What is the torque on a sphere in the steady state where the sphere rotates with the steady angular velocity ω = Ω ∞ p + O(Re s )? The undisturbed torque T ∞ follows from Eq. (8) in the Letter The disturbance torque is evaluated using the reciprocal theorem (Section II). The solution of the auxiliary problem (S6) for a spherical particle reads: The last equality defines the elements of the matrix M used in Eq. (8) in the Letter. The corresponding torque T evaluates to and the surface integral in Eq. (S10) vanishes: Sp S ∞ p · r · σ · ds = 0 . (S15) It remains to evaluate the volume integral in Eq. (S10). The Stokes solution (the solution of (S11b) for ω = Ω ∞ p ) reads:  In deriving this expression we neglected terms of the order of ω s and ω 2 s . This is justified because we know that the slip angular velocity is of order in the creeping-flow limit, see below. Setting = 0 in Eq. (S36) we find again Eq. (S17) derived in the previous Section using a different method (coordinate system fixed in the laboratory frame). We note that ω s = 0 to lowest order in Re s for = 0.
Returning to the non-spherical case we compute the angular velocity using the ansatz: To order O( ) 2 we obtain: