Rotation of a spheroid in a simple shear at small Reynolds number

We derive an effective equation of motion for the orientational dynamics of a neutrally buoyant spheroid suspended in a simple shear flow, valid for arbitrary particle aspect ratios and to linear order in the shear Reynolds number. We show how inertial effects lift the degeneracy of the Jeffery orbits and determine the stabilities of the log-rolling and tumbling orbits at infinitesimal shear Reynolds numbers. For prolate spheroids we find stable tumbling in the shear plane, log-rolling is unstable. For oblate particles, by contrast, log-rolling is stable and tumbling is unstable provided that the aspect ratio is larger than a critical value. When the aspect ratio is smaller than this value tumbling turns stable, and an unstable limit cycle is born.


I. INTRODUCTION
In this article we describe the effect of weak inertia upon the orientational dynamics of a neutrally buoyant spheroid in a simple shear flow using perturbation theory. In the absence of inertial effects the rotation of a neutrally buoyant spheroid in a simple shear was determined by Jeffery who found that there are infinitely many degenerate periodic orbits 1 , the so-called 'Jeffery orbits'. In this limit the initial orientation determines in which way the particle rotates. Fluid and particle inertia lift this degeneracy, but little is known about how this comes about. A notable exception is the work by Subramanian and Koch who have solved the problem for rod-shaped particles in the slender-body approximation 2 . We discuss other theoretical results below in Section II.
The question is currently of great interest: several recent papers have reported results of direct numerical simulations (DNS) of the problem, using 'lattice Boltzmann' methods 3-6 . These studies reveal that fluid and particle inertia affect the orientational dynamics of a neutrally buoyant spheroid in a simple shear in intricate ways. The DNS are performed at moderate and large shear Reynolds numbers, defined as Re s = sa 2 /ν where a is the largest particle dimension, s is the shear strength and ν the kinematic viscosity of the suspending fluid. DNS at very small Reynolds numbers are difficult to perform. But this limit (Re s of order unity and smaller) is of particular interest. There is a long-standing question whether or not a nearly spherical prolate spheroid exhibits stable 'log-rolling' in this limit, so that its symmetry axis aligns with the vorticity axis. It was first suggested by Saffman that this is the case 7 , in an attempt to explain Jeffery's hypothesis 1 that spheroids rotate in orbits that minimise energy dissipation. But stable log-rolling of prolate spheroids has not been found in DNS, and it has been suggested that higher Re s -corrections may explain this discrepancy 6 . The small-Re s limit is of interest also because it provides stringent tests for DNS. These reasons motivated us to derive an equation of motion that takes into account the effect of weak fluid and particle inertia. Our main result is an approximate dynamical equation for the rotation of a neutrally buoyant spheroid suspended in a simple shear flow, valid for arbitrary aspect ratios and to first order in Re s (Eq. (42) in Section IV). In the slender-body limit this equation is of the same form as the one derived in Ref. 2. In the completely inertia-free case our results reduce to Jeffery's equation 1 . We find that corrections to this limit arise from both particle inertia (centrifugal and gyroscopic forces), as well as from fluid inertia (modifying the hydrodynamic torque on the particle). The particle-inertia corrections we report here are consistent with earlier numerical and analytical results 8,9 .
Fluid-inertia corrections are taken into account to first order in Re s using a reciprocal theorem 10 . Our approach is similar to the one adopted in Ref. 2 in the slender-body limit, but our equation of motion is valid for spheroids with arbitrary aspect ratios. By linear stability analysis we determine the stabilities of the periodic orbits of this equation at infinitesimal Re s as a function of the particle aspect ratio. The stability calculation details how the degeneracy of the Jeffery orbits for a neutrally buoyant spheroid in a simple shear is lifted by weak inertia.
We find that the log-rolling orbit is unstable for prolate particles. This explains why stable log-rolling is . Shows the Cartesian coordinate systemêj, j = 1, . . . , 3. Vorticity points in the negativeê3-direction. The flow-shear plane is spanned byê1 andê2. The unit vector n points along the symmetry axis of the spheroid. Its polar angle is denoted by θ, the azimuthal angle by ϕ.
The major axis length of the spheroid is denoted by a, the minor length by b. For prolate particles the aspect ratio is defined as λ = a/b, and for oblate particles as λ = b/a. not observed in DNS 3-6 at the smallest shear Reynolds numbers accessible in the simulations. Moreover we find that tumbling in the flow-shear plane is stable for prolate particles. As the aspect ratio tends to unity there is a bifurcation: for nearly spherical oblate particles log-rolling is stable and tumbling in the flow-shear plane is unstable. There is a second bifurcation for oblate particles. At a critical aspect ratio λ c ≈ 1/7.3 tumbling becomes stable and an unstable limit cycle is born. This means that the behaviour of a very flat disk depends on its initial orientation for λ < λ c . We discuss how the shape of the limit cycle changes as the aspect ratio tends to zero. The remainder of this article is organised as follows. In Section II we give an overview over the background of the problem. Section III summarises the method employed in this article, based on a reciprocal theorem 10 . We demonstrate how to calculate the effect of particle and fluid inertia to first order, and how we use the symmetries of the problem to make it tractable. Section IV summarises our results: the equation of motion and its stability analysis. We discuss the results in Section V and conclude with Section VI.
A brief account of the main results described in this article was given in Ref. 11. Here we describe the complete derivation. We also present additional results and discussion that could not be included in the shorter format: we quote precise asymptotic formulae for small and large aspect ratios, as well as for aspect ratios close to unity. We also characterise the limit cycle that arises for λ < λ c , and compute its linear stability.

II. BACKGROUND
The question of describing the rotation of a neutrally buoyant particle in a simple shear flow has a long history. Jeffery derived an expression for the torque on an ellipsoidal (tri-axial) particle neglecting inertial effects 1 . To obtain an equation of motion for small particles he assumed that the dynamics is overdamped, that the particle rotates so as to instantaneously achieve zero torque. This gives rise to Jeffery's equation that is commonly quoted for the special case of spheroidal (axisymmetric) particles. From this equation it follows that spheroids suspended in a simple shear tumble, they stay aligned with the flow direction for some time and then switch orientation by 180 degrees. The dynamics is degenerate in that there are infinitely many different periodic orbits, the so-called 'Jeffery orbits'. The initial orientation determines which particular orbit is selected. The goal of Jeffery's calculation was to compute the viscosity of a dilute suspension of spheroids, and Jeffery hypothesised that the particles select orbits that minimise energy dissipation.
Saffman 7 pointed out that inertial effects lift the degeneracy of the Jeffery orbits, and he described the orientational dynamics of a nearly spherical particle in a simple shear taking into account fluid inertia. For prolate particles he concluded that log-rolling is stable, that tumbling in the shear plane is unstable, and that the stabilities are reversed for oblate particles. These results are stated in terms of an effective drift for the particle orientation (towards the vorticity axis for prolate particles). This conclusion supports Jeffery's hypothesis. Saffman did not take into account particle inertia. His method of calculation rests on a joint expansion in small eccentricity and Re s .
Harper & Chang 12 addressed the problem in a different way, modeling the dynamics of a rod in a simple shear in terms a dumb-bell, that is two spheres connected by an invisible rigid rod. The spheres are subject to Stokes drag and hydrodynamic lift forces 13 . This approximation neglects hydrodynamic interactions between the two spheres, as well as the unsteady term in the Navier-Stokes equations. Harper & Chang arrive at the opposite conclusion, namely that log-rolling is unstable. Since their result pertains to the slender-body limit the question is how the stability of the log-rolling orbit depends on the particle aspect ratio.
It was subsequently shown by Hinch & Leal 14 how weak rotational diffusion breaks the degeneracy of the Jeffery orbits, and their results form the basis for a large part of the work during the last decades on the rheology of dilute suspensions, see Refs. 15 and 16 for reviews.
Recently there has been a surge of interest in determining the effect of weak inertia upon a spheroid tumbling in a simple shear flow in the absence of rotational diffusion. Subramanian & Koch 2 derived an effective equation of motion for a neutrally buoyant rod in the slender-body limit to first order in fluid and particle inertia. Their calculation uses a reciprocal theorem 10 and takes into account the unsteady term in the Navier Stokes equation as well as particle inertia. The authors arrive at qualitatively the same conclusion as Harper & Chang, namely that the orientation of the rod eventually drifts towards the flow-shear plane.
In a second paper Subramanian & Koch 17 repeated Saffman's calculation for a neutrally buoyant nearly spherical particle. They used a different method, similar to the one used in Ref. 2, and come to the same conclusion as Saffman, that log-rolling is stable for nearly-spherical prolate particles.
Recent DNS 3-6 have explored the stability of log-rolling and tumbling orbits, mostly at moderate and large Reynolds numbers, and only for a small number of aspect ratios. The simulations show unstable log-rolling for prolate particles at the smallest Reynolds numbers studied. We note that You, Phan-Thien & Tanner 18 misquote Saffman when they describe their numerical results on the rotation of a spheroid in a Couette flow at Reynolds numbers of the order of 10 and larger. In the introduction of Ref. 18 it is implied that Saffman's theory 7 predicts that nearly spherical prolate particles tend to the flow-shear plane.

III. METHOD
In this section we give a brief but complete summary of our calculation. The most technical details and tabulations are deferred to appendices. We start with notation and the relevant dimensionless parameters determining the physics. Then we give the governing equations and explain how to express the hydrodynamic torque through a reciprocal theorem 2,10,19,20 . Finally we explain the perturbation scheme and list the symmetries that severely constrain the form of the solution.

A. Notation
The calculations described in this paper involve vectors and tensors in three spatial dimensions. We employ index notation with the implicit summation convention for repated indices, and we use the Kronecker (δ ij ) and Levi-Civita (ε ijk ) tensors.

B. Units and dimensionless numbers
The physics of the problem is governed by three dimensionless numbers: the shear Reynolds number Re s (measuring fluid inertia), the Stokes number St (measuring particle inertia) and the particle aspect ratio λ.
We work with dimensionless variables. The length scale is given by the particle major axis a. The velocity scale is taken to be sa where s is the shear rate. The explicit time dependence of the flow (the time scale for the unsteady fluid inertia) scales as ∼ 1/s since it is determined by the particle angular velocity because to lowest order the unsteadiness arises from the particle motion. The corresponding scale for pressure is µs, and force and torque are measured in units of µsa 2 and µsa 3 , respectively.
From these scales the dimensionless parameters are formed. As mentioned in the Introduction, the shear Reynolds number is defined as (1) Here ρ f is the density and µ = ρ f ν is the dynamic viscosity of the surrounding fluid (ν is the kinematic viscosity). The Stokes number, measuring the particle inertia, is given by the ratio of the typical rate of change of angular momentum and the typical torque: Here ρ p is the particle density. For a neutrally buoyant particle, ρ p = ρ f , we have St = Re s . We define the particle aspect ratio λ as the ratio between the length along the symmetry axis and the length transverse to the symmetry axis ( Fig. 1). That is λ > 1 denotes prolate particles while λ < 1 denotes oblate particles. Because we measure length in units of the major particle axis a, the aspect ratio λ of a prolate particle is a/b, while the aspect ratio of an oblate particle is b/a, where b denotes the length of the minor axis of the particle (see Fig. 1).

C. Equations of motion
Let n i denote the components of the unit vector n pointing in the direction of the particle symmetry axis (Fig. 1), and ω i the components of the angular velocity of the particle. Newton's second law for the orientational degrees of freedom for an axisymmetric particle reads: Dots denote time derivatives, and I ij are the elements of the moment-of-inertia tensor of the particle, and T i is the torque exerted on the particle. The moment-of-inertia tensor of an axisymmetric particle with axis of symmetry n is on the form where A I and B I correspond to the moments-of-inertia around and transverse to the symmetry axis. Using the dimensionless variables introduced in Section III B we have for a prolate spheroid (λ > 1) and for an oblate spheroid (λ < 1) We rewrite the equation of motion (3) aṡ In the final step we used the definition (4) of I ij and the equation of motion (3) forṅ i . In this paper, the T i are the elements of the hydrodynamic torque exerted on the particle by the fluid. In Section III D we formulate the hydrodynamic torque to O(Re s ) via the reciprocal theorem. In Section III E we perturbatively compute the resulting angular velocity to order O(Re s ) and O(St).

D. Calculation of the hydrodynamic torque to order Res
The straightforward approach to determine the torque on a particle in a fluid is to solve Navier-Stokes equations for the velocity and pressure fields, then compute the stress tensor and finally integrate the stress tensor over the surface of the particle. The reciprocal theorem 2,10,19 offers an alternative, and often more convenient, route to the hydrodynamic forces. In particular, we may avoid solving for the complete flow field. In this Section we specify the Navier-Stokes problem we need to solve, and explain how we use the reciprocal theorem to simplify the calculations.
Navier-Stokes problem for the disturbance flow. We consider a particle with boundary S immersed in an linear ambient flow (u ∞ , p ∞ ). Throughout this paper we express the ambient flow as Here S ∞ and O ∞ are the symmetric and antisymmetric parts of the flow gradient, given by In dimensionless variables (Section III B) the Navier-Stokes equations read Note that the unsteady and convective inertia terms come with the same prefactor in this problem. This happens because the timescale of the particle motion is the same as the timescale of the flow. The boundary condition is no-slip on the surface of the particle We introduce the disturbance field (u i , p ) from the particle If we assume that (u ∞ , p ∞ ) satisfies the Navier-Stokes equations we have the disturbance problem and the boundary conditions is expressed in the slip angular velocity Ω i = Ω ∞ i − ω i as Finally, when applying the reciprocal theorem we shall use that, by definition, the divergence of the stress tensor satisfies the following equalities: The Stokes solution. This paper concerns a spheroidal particle suspended in a linear flow. We thus need explicit solutions to Eq. (14) at Re s = 0 in this geometry. We use a finite multipole expansion 10,21 (see Appendix A). In our notation they read n B jklm = n j δ kl n m + n k δ jl n m + n j δ km n l + n k δ jm n l − 4n j n k n l n m , n C jklm = −δ jk δ lm + δ jl δ km + δ kl δ jm + δ jk n l n m + δ lm n j n k − n j δ kl n m − n k δ jl n m − n j δ km n l − n k δ jm n l + n j n k n l n m .
Here A R , B R , C R , A S , B S , C S and α are known constants that depend on the particle aspect ratio λ. The exact definition of the spheroidal multipoles Q and the values of all constants are given in Appendix A, see in particular Table III.
The reciprocal theorem. This theorem 2,10,19,20 relates integrals of the velocity and stress fields of two incompressible and Newtonian fluids. The idea is the following. Let one set of fields represent the actual problem of interest, the primary problem. Then choose the second set of fields to be an auxiliary problem with known solution, such that an integral in the theorem relates to hydrodynamic torque of the primary problem. Provided that all integrals in the theorem converge and can be evaluated, we can solve the resulting equations for the hydrodynamic torque.
The reciprocal theorem for the two sets (ũ i ,σ ij ) and (u i , σ ij ) can be stated as Here dF i = σ ij ξ j dS is the differential force from the fluid on the surface element with normal vector ξ j dS. The volume integrals are to be taken over the entire fluid volume outside the particle, and the surface integrals over all surfaces bounding the fluid volume, with surface normals pointing out of the fluid volume.
In the following we apply the reciprocal theorem to the calculation of the hydrodynamic torque on a particle.
Calculation of the torque. We choose the auxiliary problem (ũ i ,σ ij ) to be the Stokes flow around an identical particle rotating with an angular velocityω i in an otherwise quiescent fluid. Its solution is given by Eq. (17) with u ∞ = 0. The primary problem is the disturbance problem defined in Eq. (14). Inserting the boundary conditions into the reciprocal theorem yields We also used that ∂ jσij = 0. This equality holds becauseũ i is a Stokes flow. Both primary and auxiliary velocity fields vanish as |r| → ∞, therefore both integration surfaces are only the particle surface S. Note that the surface integrals are to be taken with surface normals out of the fluid domain, so that dF i is the differential force exerted on the particle by the fluid. In the integrals we identify the hydrodynamic torque on the particle, it is given by It follows: The auxiliary torqueT j together with the surface integral add up to the Jeffery torque 1 T (0) j : The constant c ξ is given in Table III in Appendix A. The contribution evaluates to zero for any linear flow u ∞ i . It follows that Eq. (21) becomes Sinceũ i is linear inω j , this variable can be eliminated. We finally obtain: Thus far we have made no approximation, and Eq. (25) is exact, the difficulty lies in evaluating the Navier-Stokes disturbance flow u . This is a complicated non-linear problem since T (0) j ,Ũ ij and f i (u ) all depend on the direction n and upon the angular velocity ω of the particle. The flow equations thus couple non-linearly to the rigid body equations of motion for the particle. In the following we solve this system of equations in perturbation theory valid to first order in St and Re s .

E. Perturbative calculation of the particle angular velocity
In this section we determine the angular velocity ω of the particle to lowest order in St and Re s , assuming that both St and Re s are small, so that Re s St is negligible. We recall the equation of motion (7) for the particle orientation, and insert the expression for the hydrodynamic torque obtained in Section III D: Now we expand the angular velocity as Next, we insert these expansions into the equation of motion (7) and collect terms of equal order in St and Re s : In the last term it is understood that the volume integral need only be evaluated to O(1), so that we may use the Stokes flow solutions for u . The first equation gives the Jeffery angular velocity ω The dynamics of n is to lowest order given bẏ From Table III in Appendix A we infer that for both prolate and oblate spheroids This shows that Eq. (31) is Jeffery's equation 1 for the orientational dynamics of a spheroid in a simple shear. The two remaining equations in (29) may be inverted to Eq. (28) together with Eqs. (30) and (33) yield the effective angular velocity under the effect of weak particle and fluid inertia. From the equation of motion (7) we define the effective vector fielḋ This vector field describes the time evolution of n. The first term is the Jeffery vector field (31). The two new terms represent the effects of particle inertia and fluid inertia. The terms due to particle inertia are straightforward to evaluate directly, but the volume integral in Eq. (33) is very tedious to evaluate. To make the calculation feasible we exploit the symmetries of the problem.

F. Symmetries of the effective equation of motion
Both correction terms in Eq. (34) are quadratic in the ambient flow gradient tensor A ∞ ij . In other words, they are on the formṅ where the tensorial coefficients C (i) ijklm are composed of the remaining available tensor quantities: n p and δ pq (ε ijk is already used in O ∞ ij ). We make an exhaustive enumeration of all possible combinations, and then use the symmetries listed in Table I to remove or combine items in the list. For example we start by letting C (1) ijklm = P η (P) 1 n p1 δ p2p3 δ p4p5 + η (P) 2 n p1 n p2 n p3 δ p4p5 + η (P) 3 n p1 n p2 n p3 n p4 n p5 where the sum is over all 5! = 120 permutations P of (i, j, k, l, m), and η (P) i are unique coefficients for each term. We include only odd powers of n i as any even terms would break the particle inversion symmetry. We then insert this enumeration into the first term of Eq. (35), and contract and apply the first three symmetries in Table I until we reach a list of unique candidate terms. In this case the only two unique terms turn out to be O ij O jk n k and n i n j O jk O kl n l . Finally, we use the fact that the equation of motion may not change the magnitude of the unit vector n. This constraint forces the coefficients of the two unique terms to be the same magnitude but opposite sign. Upon renaming the coefficients ±α 1 we get the first term in Eq. (37). The other terms are derived similarly by inserting (36) into the other terms in Eq. (35). The result contains only six independent terms:ṅ Here the scalar functions α 1 , . . . , α 6 are linear in St and Re s , and depend on the aspect ratio λ in a nonlinear (and unknown) way. These coefficients are determined by evaluating the vector field in Eq. (34) for six independent directions of n and solving the resulting system of linear equations for α 1 , . . . , α 6 .
Thus, for the case of the simple shear flow it suffices to evaluate the effective vector field in Eq. (34), in particular the volume integral in Eq. (33), with four independent values of n in order to solve for the unknown scalar coefficients β 1 , . . . , β 4 .

G. Evaluation of the volume integral in Eq. (33)
The volume integral in Eq. (33) contains four distinct terms: ∂ t u i , represents unsteady fluid inertia, and the three terms u ∞ j ∂ j u i + u j ∂ j u ∞ i + u j ∂ j u i represent convective fluid inertia. We compute these four terms using the explicit Stokes-flow solutions (17). While the Stokes flow has no explicit time dependence, both particle direction n and angular velocity ω do. Thus each occurence of n k and ω k has to be differentiated to compute the contribution due to unsteady fluid inertia. The differentiation and tensor contractions are implemented by a custom set of pattern matching rules in Mathematica R . The calculation is both long and error prone. We have therefore automated every possible step, including solving the Stokes-flow equations.
We demonstrate the remainder of the procedure by a small example. Consider the contribution in thê e 3 -direction of Eq. (33) due to unsteady fluid inertia: We first perform the time derivatives on (17) in the manner explained above. Then we insert the components of n, and the explicit form of the shear flow (38). At this point we can explicitly perform the sum over all repeated indices. The result in this example consists of 858 terms, after collecting terms with same spatial dependence. The terms have a prefactor that stem from the Stokes-flow coefficients (see Appendix A), and a spatial dependence coming from r i and the spheroidal integrals J n m and K n m (see Appendix B). For n = [1/2, √ 3/2, 0] a typical term looks like this: We note that the only spatial dependence on the azimuthal angle around the symmetry axis of the body comes from factors r i . We introduce a rotated coordinate system in which r i = R ji r j , such that r 1 is along the particle symmetry axis (see Appendix B). This change of basis enables integration of one spatial coordinate. After this operation 260 terms still remain which we program Mathematica R to express in spheroidal coordinates (Appendix C) and integrate over the remaining two spatial coordinates.
As a consistency check we have also evaluated the volume integral numerically over all three spatial dimensions by converting to spheroidal coordinates and choosing a specific value of λ. For extreme values of λ the numerics are difficult, nevertheless they serve as a check for a wide range of aspect ratios (see markers in Fig. 2).
We compute the contributions to β α from three sources: particle inertia, unsteady fluid inertia and convective fluid inertia. Although the result is only valid for neutrally buoyant particles (Re s = St), it is interesting to consider the contributions separately: The contribution from particle inertia is straightforward to compute and can be expressed in closed form as The coefficients on the r.h.s. of these equations are tabulated for both prolate and oblate spheroids in Table III in Appendix A. The coefficients in Eq. (44) are shown as dotted lines in Fig. 2. The expressions for the contributions from fluid inertia are very lengthy and not particularly instructive. We therefore present the full result graphically as function of aspect ratio λ in Fig. 2. In addition we give the asymptotic behavior of all contributions to β α in Table II in three limiting cases: thin oblate particles  The asymptotic results for λ → 0, λ → ∞, and for | | → 0 are shown as red dashed lines in Fig. 2.

B. Linear stability analysis at infinitesimal Res
The effective equations of motion (42) have two special polar angles θ across which no orbit may pass, regardless of the values of β α . These angles are θ = 0 (the vorticity direction) and at θ = π/2 (the flow-shear plane). In the Jeffery dynamics (Re s = St = 0) the two orbits are called 'log-rolling' and 'tumbling', and they are both marginally stable, just like all other Jeffery orbits. When the β α are non-zero but infinitesimal, the log-rolling and tumbling Jeffery orbits still exist for any finite aspect ratio, but their stabilities change.
We quantify how particle and fluid inertia lift the degeneracy of the Jeffery orbits by computing the stability exponents γ for the log-rolling (γ LR ) and tumbling (γ T ) orbits. The stability exponent is the exponential growth rate over one period of the orbit: where T p = 4π/ √ 1−Λ 2 is the Jeffery period. As Re s → 0 we find For Re s = St these two exponents are shown as function of particle aspect ratio in Fig. 3. Also shown are their limiting behaviours in the thin oblate limit (λ → 0) in the nearly spherical limit ( → 0) and in the thin prolate limit (λ → ∞) Fig. 3 shows that prolate spheroids of all aspect ratios are unstable at the log-rolling position, and stable at the tumbling orbit. For nearly spherical particles there is a bifurcation: log-rolling and tumbling switch stabilities. For oblate spheroids the log-rolling position is stable for any aspect ratio. For oblate particles there is a second bifurcation at λ c ≈ 1/7.3 where the tumbling orbit becomes stable. Clearly, this behavior is caused by the convective inertia of the fluid (see the dash-dotted line in Fig. 3). For sufficiently oblate particles both log-rolling and tumbling orbits are stable, and the long-time dynamics depend on the initial orientation of the particle. Between the two now stable orbits a new unstable limit cycle is born, separating the two basins of attraction. Fig. 4 shows how the shape of this limit cycle depends upon the particle aspect ratio. Close to the bifurcation the limit cycle lies in the neighbourhood of the tumbling orbit. But as λ → 0 the limit cycle approaches the log-rolling orbit. We have computed the stability exponent of the limit cycle at infinitesimal Re s by numerically integrating Eqs. (42). The result is shown in Fig. 5. We see that γ LC > 0, and its magnitude is of the same order as that of γ T .

V. DISCUSSION
Effective equation of motion. Eq. (42) is an effective equation of motion for the orientational dynamics of a neutrally buoyant spheroid in a simple shear flow. How the dynamics depends upon the particle aspect ratio is determined by four coefficients β 1 , . . . , β 4 . Fig. 2 shows the four functions β α (λ). Limiting behaviours of the β α are tabulated in Table II. We see that the β-coefficients tend to zero as λ → ∞, but they approach constants as λ → 0. In both limits the contribution from particle inertia must tend to zero because the volume of the particle does. The effects of fluid inertia vanish as λ → ∞ because the particle effectively disappears in the slender-body limit, the perturbation caused by the particle decreases as ∼ 1/ log λ as the asymptotic form in Table II shows. We remark that the leading-order term in this asymptotic form makes a substantial correction to the slender-body theory for aspect ratios of order 30.
An oblate particle, on the other hand, always presents no-slip boundaries to the fluid, with an area of the order of ∼ a 2 as λ → 0. Therefore the contribution of fluid inertia approaches a constant. We note that the asymptotic forms of the coefficients β α listed in Table II yield accurate values for λ < 1/30 and λ > 30, as  We see in Fig. 2 that the particle-inertia contribution to the coefficients β α is always much smaller than the fluid-inertia contributions. In general both unsteady and convective fluid inertia contribute, and it would be qualitatively wrong to neglect one of these terms. This is due to the fact that the timescale of the particle motion is the same as the timescale of the flow, and it raises the question under which circumstances both effects may matter for the tumbling of small particles in unsteady flows, and in particular in turbulence.
Linear stability analysis at infinitesimal Re s . The stability exponents of tumbling and log-rolling orbits are shown in Fig. 3. We find that the log-rolling orbit is unstable for prolate spheroids of any aspect ratio, tumbling is stable for prolate spheroids, and no other orbit exist at infinitesimal Re s . For moderately oblate particles with aspect ratios λ > λ c ≈ 1/7.3 the stabilities are reversed: log-rolling is stable, tumbling is unstable, and no other periodic orbits exist for infinitesimal Re s . At λ = λ c there is a bifurcation where an unstable periodic orbit is born close to the tumbling orbit, which in turn becomes stable. As λ becomes even smaller, the unstable orbit moves closer to the log-rolling orbit (Fig. 4). We remark that the asymptotic forms (47) and (49) of the stability exponents yield very accurate approximations for the log-rolling exponent, save for aspect ratios close to unity. For the tumbling exponent the asymptotes do not work equally well.
Our results are in agreement with results of recent DNS studies 3-6,18 determining the orientational dynamics of a neutrally buoyant spheroid in a simple shear flow. These studies are conducted for a number of different aspect ratios with shear Reynolds numbers ranging from moderate to large. At the smallest values of Re s accessible in the DNS no stable log-rolling is found for prolate spheroids of any aspect ratio. For oblate particles with aspect ratio λ = 1/5 DNS show stable log-rolling and unstable tumbling at the smallest Re s that were simulated 6 , also in agreement with our results. There are no simulations for particles for λ < λ c at small Re s . Saffman 7 predicted that log-rolling is stable for nearly spherical prolate particles, at variance with the behaviour described above. We do not know why the original calculation fails to give the correct stability of log-rolling. Since no details of the calculation are given it is difficult to figure out the precise origin of this discrepancy. Subramanian & Koch 17 also computed the stability of the log-rolling orbit for nearly spherical particles and came to the same conclusion as Saffman, different from ours. We have compared the smalllimit of our calculation to the results of Ref. 17 and find that the particle-inertia correction to the equation of motion agrees, Eqs. (3.15) and (3.16) in Ref. 17. But the fluid-inertia correction does not satisfy the symmetries of the problem. We believe that this explains the discrepancy.
We have independently calculated the stability of log-rolling for nearly spherical particles by expanding the particle-angular velocity jointly in and Re s , using spherical harmonics as a basis set 22 . The results of this calculation agree to order with the results presented above. Further we have checked that the particle-inertia correction in Eq. (42) is consistent with the results obtained in Ref. 9. We also compared the slender-body limit of our results to the prediction of Subramanian & Koch for the dynamics of slender fibres 2 and found that the fluid-inertia corrections agree (up to a factor of 8π).
These observations indicate that the results presented in this paper are correct, explain the results of DNS and resolve the puzzle concerning the stability of log-rolling of spheroids in a simple shear at small Re s .
A new benchmark for DNS at small Re s . Recently a number of groups have developed DNS codes based on the lattice Boltzmann method to simulate the dynamics of particles in flows 3-6 . Much effort is spent on validating the model, studying for instance the effects changing grid size, time step, size of the simulation box, and so forth. The benchmark adopted is often the question whether Jeffery orbits are seen for a neutrally buoyant spheroid in a simple shear at small Reynolds numbers. But the limit Re s = 0 can never be strictly reached in the simulations. DNS at small values of Re s (specifically: in the linear regime), by contrast, allow precise comparisons with the results obtained in this paper. One could for instance compare trajectories, stability exponents, and period times. We thus expect that our results can serve as benchmarks for present and future DNS codes.

VI. CONCLUSIONS
In this paper we have derived an effective equation of motion for the orientational dynamics of a neutrally buoyant spheroid suspended in a simple shear flow. The equation is valid for arbitrary aspect ratios and to linear order in Re s , at small but finite shear Reynolds numbers. The effective equation of motion allows us to determine how the degeneracy of the Jeffery orbits is lifted by weak inertial effects. We have determined the bifurcations that occur at infinitesimal Re s as the particle aspect ratio changes. For prolate spheroids log-rolling is unstable, for oblate spheroids it is stable. Tumbling in the shear plane is stable for prolate particles and unstable for nearly spherical oblate particles. For thin disks with aspect ratios λ < 1/7.3, both log-rolling and tumbling are stable. An unstable limit cycle separates the basins of attraction of the periodic orbits.
Our results imply that tumbling and log-rolling orbits survive a finite perturbation whose magnitude depends on the aspect ratio λ. It would be of interest to derive a bifurcation diagram in the λ-Re s -plane for small Re s . We plan to determine how the small-Re s region of this diagram connects to the intricate bifurcation patterns that were found by Rosén, Lundell & Aidun 5 at larger shear Reynolds numbers. We expect that the results summarised here can guide numerical computations with the lattice Boltzmann method that become difficult at small Re s and large aspect ratios. both for prolate and oblate spheroids. In order to extract the six independent equations for the six remaining coefficients we exploit that the boundary condition must be satisfied for any choice of n j , Ω j and S ∞ jk . First, with S ∞ jk = 0 we contract Eq. (A4) with n i , ε ijk n j Ω k and ε ipq n q ε pjk n j Ω k . Secondly, with Ω j = 0, we contract Eq. (A4) with n i , S ∞ ij n j and finally ε ijk n j S ∞ kl n l . These six equations together have only one solution. We tabulate the resulting expressions for both oblate and prolate spheroids in Table III. Computing the torque on a body due to this flow is straightforward, because by construction 21 the torque on a body due to the rotlet flow u i = G R ij,k ε jkl A l is T R l = −16πA l , where A l is the rotlet strength. The minus sign is due to the fact that the torque is exerted on body by the flow. To compute the torque from the spheroidal rotlet (A10) we linearly superpose the contributions from all the contained rotlets. The torque from the flow u i = Q R ij,k ε jkl B l is therefore The factor c ξ depends only on the aspect ratio of the particle (see Table III).  10 . We tabulate them here for convenience, and because our conventions differ slightly from those adopted in Ref. 10. We remark that some of the coefficients tabulated here assume imaginary values. All physical quantities come out to be real-valued.