Orientational dynamics of weakly inertial axisymmetric particles in steady viscous flows

The orientational dynamics of weakly inertial axisymmetric particles in a steady flow is investigated. We derive an asymptotic equation of motion for the unit axial vector along the particle symmetry axis, valid for small Stokes number St, and for any axisymmetric particle in any steady linear viscous flow. This reduced dynamics is analysed in two ways, both pertain to the case of a simple shear flow. In this case inertia induces a coupling between precession and nutation. This coupling affects the dynamics of the particle, breaks the degeneracy of the Jeffery orbits, and creates two limiting periodic orbits. We calculate the leading-order Floquet exponents of the limiting periodic orbits and show analytically that prolate objects tend to a tumbling orbit, while oblate objects tend to a log-rolling orbit, in agreement with previous analytical and numerical results. Second, we analyse the role of the limiting orbits when rotational noise is present. We formulate the Fokker-Planck equation describing the orientational distribution of an axisymmetric particle, valid for small St and general Peclet number Pe. Numerical solutions of the Fokker-Planck equation, obtained by means of expansion in spherical harmonics, show that stationary orientational distributions are close to the inertia-free case when Pe St<<1, whereas they are determined by inertial effects, though small, when Pe>>1/St>>1.


Introduction
Suspensions of rigid particles are abundant in both nature and technology, and are consequently studied in many disciplines of science. Examples are ash and ice crystals in the atmosphere [1], plankton in the ocean [2,3], or fibres in paper pulp [4]. In many situations the concentration of the suspended particles in the carrying fluid is small, so that particles do not affect the underlying fluid flow. In this case, considering the fluid flow field as given, a question of fundamental interest is to determine the resulting trajectory and orientation of a particle. The orientational dynamics of axisymmetric particles is often described by the Jeffery equation of motion [5], the rotational analogue to simple center-of-mass advection along streamlines. This equation of motion is valid in the limit of vanishing Stokes number St and particle Reynolds number Re (both numbers are defined below), representing particle and fluid inertia respectively. The Jeffery equation of motion is widely used to predict mechanical [6,4,7] as well as optical properties [8,9,10,11] of dilute suspensions of inertia-free rod-like or disk-like particles. It is also used to describe the orientational motion of single particles [3,12], and the alignment and tumbling of axisymmetric particles in turbulence [13,14].
In the case of simple shear flow an inertia-free particle rotates in one of infinitely many possible Jeffery orbits. In other words, the solution of the Jeffery equation is degenerate, forming a one-parameter family of closed orbits. Which orbit is followed depends upon the initial condition. By taking into account the inertia of the particle in numerical simulations, it was observed that attracting orbits exist [15,16], where particles tend to align their small axis towards the vorticity axis. In the special case of nearly spherical particles the existence of attracting orbits has been demonstrated analytically [17]. The significance of this inertial effect is that it breaks the degeneracy of the Jeffery orbits.
If, in addition, particles are subject to a random force, the resulting orientational distribution is expected to be determined by a balance between the dispersing random forces and the aligning inertial forces.
The goal of the present study is to investigate the orientational dynamics of axisymmetric particles with small inertia, to describe and quantify the attracting orbits, and to analyse their effect on the orientational distribution when noise is present.
To accomplish this goal we first derive an equation of motion for the orientation of an axisymmetric particle in steady flows in the absence of noise, valid for small values of the Stokes number St. We define the Stokes number by St = ms/bµ, where m is the particle mass, s is a typical flow-gradient rate, b is the minor particle length and µ the dynamic viscosity of the fluid (see Sec. 2 for details). We make use of the fact that St is small, so that the angular velocity of the particle is a small perturbation of the angular velocity in the absence of particle inertia. Our results, valid for any aspect ratio, generalise a result for nearly spherical spheroids reported in Ref. [17]. In Ref. [17] fluid inertia, as mea-sured by the particle Reynolds number, is also considered. The particle Reynolds number is defined by Re = (κρ f /ρ p )St, where ρ f and ρ p denote fluid and particle density. The 'shape factor' κ is κ = 1 for a disk-shaped particle, and κ = λ for a rod-like particle (λ is the aspect ratio of the particle). Finite-Re corrections are important for neutrally buoyant particles [ρ p /ρ f = O(1)], because the fluid inertia induces O(Re)-corrections in the hydrodynamic torque experienced by the particle [17,18]. The results obtained here are valid for heavy particles such that ρ p κρ f . The new equation of motion allows us to derive a Fokker-Planck equation determining the combined effects of particle inertia and noise upon the orientational distribution of the particles. This equation describes the evolution of an ensemble of weakly inertial particles in a steady flow, subject to random Brownian rotations. It is solved in the case of a simple shear flow by means of expansion in spherical harmonics. We find that the stationary orientational distributions differ significantly from the inertia-free case.
The rest of this paper is organised as follows. Asymptotic equations of motion valid for small St in the absence of noise are obtained in Sec. 2. Then, in Sec. 3, we investigate the case of a simple shear flow. We evaluate the convergence towards one of two possible periodic orbits by calculating the corresponding Floquet exponents. Finally, in Sec. 4, we formulate the Fokker-Planck equation describing the combined effect of noise and weak particle inertia. This equation is then solved numerically and the stationary distribution is discussed. Sec. 5 summarises our conclusions.

Equation of motion
The aim of this Section is to derive an equation of motion for the orientation of a small axisymmetric particle, valid to first order in St. The orientation is described by a unit vector n pointing in the direction of the symmetry axis of the particle. In the following we derive an equation of the formṅ = F 0 (n) + StF 1 (n) + O(St 2 ), where the dot denotes the time derivative, and the functions F i depend on the flow gradient A = ∇u. We decompose A into its symmetric and antisymmetric parts The matrix O is related to the vorticity ∇ ∧ u of the flow: Ox = Ω ∧ x, where Ω = ∇ ∧ u/2.
We start with Newton's equation of motion for the orientational degrees of freedom for a rigid body. Let ω denote the angular velocity of the particle, theṅ where M is the moment of inertia, and T the torque applied to the body. The hydrodynamic torque on a particle depends on the flow configuration at the particle position. Therefore, when studying inhomogenous flows such as turbulent channel flows, Eq. (1) must be completed with equations for the particle center-of-mass motion [19]. Here we consider homogeneous flows where the center-of-mass motion may be neglected. With the hydrodynamic torque calculated by Jeffery [5], Eq. (1) describes the rotation of both axisymmetric [15] and asymmetric [16] particles. In the following we consider the case of axisymmetric particles. An axisymmetric object has two principal moments of inertia, around and transverse to the axis of symmetry. The moment of inertia tensor is therefore on the form where I is the identity matrix, and X I , Y I are constants depending on the shape of the particle. Here and in the following our notation is similar to that used in the book by Kim and Karrila [20]. The difference, and the reason we keep the tildes, is that here the quantities have dimensions. For instance a spheroid with lengths a and b and mass m has X I = 2mb 2 /5 and Y I = m(a 2 + b 2 )/5. Further, the hydrodynamic torque depends linearly on the flow gradient via resistance tensors which are determined solely by the shape of the particle [6,5,20] Here X C , Y C and Y H are shape constants. The relation to the dimensionless shape functions given in [20] is, for instance, X C = 8πµa 3 X C . The shape functions for spheroids are given on p. 64 (prolate case) and p. 68 (oblate case) of Ref. [20]. Inserting the explicit forms Eq. (2) and Eq.
An axisymmetric body has two different lengths a > b, and we choose to base the Stokes number on the lesser of the two, such that St = ms/bµ. Here m is the mass of the particle, s is a typical shear rate and µ is the dynamic viscosity of the surrounding fluid. We introduce dimensionless variables t = t /s, A = sA , ω = sω , and drop the primes. We now assume that for small St the angular velocity ω is almost the same as in the unperturbed case: Inserting into Eq. (4) and comparing order by order in St we find Here the dimensionless coefficients The lowest order, ω 0 is found from Eq. (5) to be which we recognise as the Jeffery angular velocity [20]. In particular, for a spheroid of aspect ratio λ = a/b the material constant α S /α Y = (λ 2 − 1)/(λ 2 + 1) ≡ Λ. To recover the well-known Jeffery equation forṅ, note thaṫ Solving for the next order ω 1 in Eq. (6) we find This expression for the correction to the angular velocity is valid for any axisymmetric body, specified by the parameters α i . For the case of a simple shear flow and nearly spherical spheroids Eq. (7) agrees with earlier results (Eq. (3.12) in Ref. [17] is equal to Eq. (7) to first order in Λ). In the case of spheroids, both prolate and oblate, we insert the appropriate shape functions from [20] and observe that We find our final result for the equation of motion for n to order O(St) aṡ Eq. (8) is valid as an approximation to Eq. (1) when particles are weakly inertial. The vector n is advected in an effective vector field F 0 (n) + StF 1 (n) that differs from Jeffery's equation by a term linear in St. This approximation does not capture the formation of caustics that occur when particles with different orientational histories arrive at the same orientation but with different angular velocities [21]. From Eq. (8) we see that the geometrical factor multiplying the inertial correction is 1/α Y . Therefore, in the following, we use the effective Stokes number

Simple shear flow: stability analysis
The simple shear flow, u(r, t) = syx, is interesting because it is the local linearisation of any flow with parallel streamlines, like pipe flows. It is also the flow in an ideal Couette rheometer used to measure the viscosity of fluids. The solutions of Eq.   Eq. (8) is a drift towards a stable periodic orbit. In this section we will characterise this orbit drift for general particle shapes λ by computing the Floquet exponents of the limiting orbits.
The modified vector field on the right hand side of Eq. (8) allows for two limiting orbits. First the equatorial orbit, where the n-vector moves in the plane perpendicular to vorticity (usually referred to as tumbling), and second the polar orbit, where the n-vector is fixed parallel to the vorticity (referred to as log rolling). We introduce a spherical coordinate system where θ is the polar angle, measured from the vorticity vector, and ϕ is the azimuthal angle, measured from the flow direction, see Fig. 2. Then the two limiting orbits correspond to θ = 0 and θ = π/2. In spherical coordinates, Eq. (8) takes the following form for a simple shear flow: sin θ cos θ cos 2 (ϕ) 4Λ sin 2 θ sin 2 ϕ − Λ + 1 .
The Floquet exponent is obtained by time-averaging Eq. (9) along the unperturbed periodic orbit (ϕ 0 , θ 0 ): The last equality follows from the fact that in one period, ϕ travels monotonically through ϕ = 0 to −2π. Upon evaluating Eq. (10) we obtain for the polar, log rolling orbit and for the equatorial, tumbling orbit A positive Floquet exponent γ indicates an unstable orbit, and a negative an attracting, stable orbit. The exponents (divided by St) are shown as a function of particle shape λ in Fig. 3. We observe that the stable limiting orbit for flat, diskshaped particles is the polar, log rolling orbit, while rod-shaped particles are attracted to the equatorial, tumbling orbit. This is consistent with the intuition that the most stable rotation is that with the mass distributed far from the axis of rotation. Our result explains the orbit drift observed in [15] at small values of St. In Ref. [17] this drift (at small Stokes number and vanishing particle Reynolds number) was explained for nearly spherical particles. Here we generalise the result to arbitrary axisymmetric particles.

Competition between noise and inertia
The indeterminacy of the Jeffery orbits poses a problem for predictions involving the orientational distribution of particles. Jeffery [5] hypothesised that orbits corresponding to minimal energy dissipation should be preferred. Another approach, pursued by a number of authors [23,24,6], shows that thermal noise, after a long time, removes the dependence on initial conditions and produces a stationary distribution of the n vector. In particular, in [24] it is shown that the orientational distribution converges to a function independent of noise strength, as the noise strength is reduced. On the contrary, in the preceding sections we have shown that in the absence of noise a weakly inertial particle will rotate into a specific stable orbit. In this section we consider what effect weak particle inertia have on the orientational distribution.
The time scale for inertia to have an effect is equal to the inverse Floquet exponent, s/γ ∼ 1/ St. The corresponding time scale for diffusion is sτ D ∼ Pe, where the Péclet number Pe = s/D measures the noise strength, or diffusion constant D = k B T/ Y C , relative to the shear strength [20]. Thus, when τ D γ ∼ Pe St 1 we expect noise to determine the stationary distribution, as described by previous authors [6]. On the other hand, when Pe St 1 we expect a distribution dominated by the inertial effects, with a peak around the stable orbit described in the previous section. In order to quantify these expectations we need to compute moments of the stationary distribution for various Pe, St and λ.
The Fokker-Planck equation describing the evolution of the orientational prob-ability density P(n, t) is Here ∂ n is the gradient on the sphere defined as ∂ n ≡ (I − nn T )∇. The boundary condition of Eq. (11) is that the probability density P must integrate to unity over the sphere. The Fokker-Planck equation can be solved numerically by expansion in spherical harmonics, the details of this procedure are described in Appendix A. For a discussion of the solutions of Eq. (11) we compute the stationary average sin 2 θ , where θ is the polar angle from the vorticity vector of the simple shear flow (see Fig. 2). Recall that weakly inertial rod-shaped particles exhibit a stable orbit at θ = π/2, and disk-shaped particles at θ = 0. Moreover, values of sin 2 θ have been tabulated (Table 6c in Ref. [6]) for different values of Pe and λ. This allows us to validate our method at St = 0.
In Fig. 4 we show the angular average sin 2 θ as a function of Pe for different values of St, and for two particle shapes: prolate λ = 5 and oblate λ = 1/5. We note the following features. At strong noise, that is small Pe, the moment is independent of St and approaches the uniformly random result sin 2 θ = 2/3. Next we see that the curve for St = 0 (red curve) levels out and converges at a finite value of Pe, as predicted in Ref. [24]. Finally we see that for finite St, the curves coincide with the St = 0 curve up to PeSt ≈ 1, but then deviates and instead converges towards sin 2 θ = 1 (prolate particle), or sin 2 θ = 0 (oblate particle) when PeSt 1.
The limit of St → 0 is singular: the stationary distribution at low noise and weak inertia differs substantially from the one predicted by a St = 0 calculation, although the time to arrive at the steady state becomes longer for weaker inertia. The most important observation is that the stationary average at low noise strength but without inertia is the same for prolate and oblate particles, whereas the inertial correction has opposite directions for the two different shapes. We emphasise this observation in Fig. 5

Conclusions
The orientational dynamics of weakly inertial axisymmetric particles in a steady flow has been investigated. We have derived an approximate equation of motion for the unit axial vector along the particle symmetry axis, assuming that the Stokes number is small. The fact that we have obtained this equation of motion in explicit form (valid for any axisymmetric and weakly inertial particle in any steady linear viscous flow) has made it possible to address two important questions, both pertaining to the case of a simple shear flow.
First, we have shown that axisymmetric particles drift towards either the logrolling (disk-shaped particles) or the tumbling orbit (rod-like particles). We have quantified this result for any particle aspect ratio by calculating the Floquet exponents of the periodic orbits to leading order in St. Our result generalises a result for nearly spherical particles in Ref. [17] to arbitrary axisymmetric particles.
Second, the approximate equation of motion derived in Sec. 2 has allowed us to formulate the Fokker-Planck equation describing the orientational distribution of axisymmetric, weakly inertial particles in a shear flow when rotational noise is present. Numerical solutions of the Fokker-Planck equation show that stationary orientational distributions are close to the inertia-free case when Pe St 1, whereas they are determined by inertial effects, though small, when Pe 1/ St

1.
An important open question is to describe the effect of fluid inertia (at small but finite particle Reynolds number) upon the orientational motion for particles with arbitrary aspect ratios, in competition with the effects of particle inertia and rotational noise.
The calculations summarised in this article are restricted to axisymmetric particles. But our approach is readily generalised to asymmetric (triaxial) particles.

Inserting the expansion yields
In order to solve Eq. (A.7) it remains to compute the matrix elements of the opera-torĴ. This is achieved by expressingĴ as a combination of the angular momentum operators. How the angular momentum operators act on the spherical harmonics is well known, see for example Ref. [25]. The angular momentum operatorL is given bŷ In terms ofL we have The states |l, m are the eigenfunctions ofL 2 andL 3 with eigenvalues l(l + 1) and m. Now, to evaluate the drift term ∂ nṅ we need to know howL 1 ,L 2 , andn act on |l, m . The first two are known through the use of ladder operators, defined bŷ and their effect isL Next, the drift term we consider is a polynomial inn 1 ,n 2 andn 3 (the components ofn), we therefore need to evaluate the action of a monomialn α The action ofŶ q p is computed with the Clebsch-Gordan coefficients (see for example Ref. [26], p.216) × l 1 l 2 ; 00|l 1 l 2 ; l0 l 1 l 2 ; m 1 m 2 |l 1 l 2 ; lm .
The Clebsch-Gordan coefficients are denoted by l 1 l 2 ; m 1 m 2 |l 1 l 2 ; lm (see e.g Eq. (3.7.44) of Ref. [26]), and they are available in Mathematica by the function ClebschGordan[l1,m1,l2,m2,l,m]. Finally, we order the operators, so that as many terms as possible cancel before we actually begin evaluating the operators. In particular we want to reordern i againstL i , and we make use of the commutator where ε i jk is the Levi-Civita tensor. Let us consider an example. In the special case ofĴ = Pe −1 ∂ 2 n we see that J = −Pe −1L 2 . It follows: J|l, m = −Pe −1L 2 |l, m = −l(l + 1)Pe −1 |l, m .
This means that the diffusion operator exponentially suppresses all modes with p > 0, and the lowest mode p = 0 is determined by the normalisation condition Eq. (A.8). This is the solution of the diffusion equation on the sphere, with the uniform distribution as steady state. Now take the for the drift term the standard Jeffery equationṅ = On + Λ(Sn − nn T Sn) = F 0 (n), and call this operatorĴ = −∂ n F 0 (n) + Pe −1 ∂ 2 n . We find The fact that the operatorĴ is real implies that the matrix elements must have the symmetry p, q|Ĵ|l, m = (−1) q p, −q|Ĵ|l, m .
We can understand the reason for this because the system of differential equations in Eq. (A.7) needs to preserve the condition that P is real, as stated in Eq. (A.5). It implies we have only to compute half of the matrix elements. Upon including the correction Eq. (8) due to weak particle inertia, the expression for the operatorĴ becomes too lengthy to include here. However, we have implemented the identities and rules described in this appendix in a Mathematica notebook [27]. Given an operatorĴ expressed in terms ofn andL, it converts the expression into sums of angular momentum operators as exemplified above. We also provide an additional Mathematica notebook that, given the matrix elements, assembles a sparse matrix J, up to a desired order l max . This matrix can then be used to solve the truncated version of Eq. (A.7): where c is a vector with elements c m l . Alternatively, the matrix can be solved for the stationary solution of c. All results shown in this paper were computed using l max = 400, which leads to very good convergence in all cases shown. Only when solutions approach delta peaks, as for example the limiting stable orbits in the large PeSt case, the expansion procedure converges very slowly.