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 and log-rolling is unstable. For oblate spheroids, by contrast, log-rolling is stable and tumbling is unstable provided that the particle is not too disk-like (moderate asphericity). For very flat oblate spheroids, both log-rolling and tumbling are stable, separated by an unstable limit cycle.

## 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 inertias 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 Sec. II.

The question is currently of great interest: several recent papers have reported results of direct numerical simulations (DNSs) of the problem, using “lattice Boltzmann” methods.^{3–6} These studies reveal that fluid and particle inertias affect the orientational dynamics of a neutrally buoyant spheroid in a simple shear in intricate ways. The DNSs 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. DNSs 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 Sec. 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) and 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 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 Sec. 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 inertias 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 Sec. V and conclude with Sec. 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 and 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°. 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 the 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 the particle inertia. His method of calculation rests on a joint expansion in small eccentricity and Re_{s}.

Harper and 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 subjected 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 and 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 and 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 and 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 inertias. 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 and Chang, namely, that the orientation of the rod eventually drifts towards the flow-shear plane.

In a second paper, Subramanian and 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 came to the same conclusion as Saffman that log-rolling is stable for nearly spherical prolate particles.

Recent DNSs^{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 Yu, Phan-Thien and 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 repeated 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 the 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

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

**pointing in the direction of the particle symmetry axis (Fig. 1) and**

*n**ω*the components of the angular velocity of the particle. Newton’s second law for the orientational degrees of freedom for an axisymmetric particle reads

_{i} Dots denote time derivatives, and *I _{ij}* are the elements of the moment-of-inertia tensor of the particle and

*T*are the components of the torque exerted on the particle. The elements of the moment-of-inertia tensor of an axisymmetric particle with axis of symmetry

_{i}**take the form**

*n* where *A ^{I}* and

*B*correspond to the moments-of-inertia around and transverse to the symmetry axis. Using the dimensionless variables introduced in Sec. III B, we have for a prolate spheroid (

^{I}*λ*> 1),

and for an oblate spheroid (*λ* < 1),

We rewrite equation of motion (3) as

In the final step, we used definition (4) of *I _{ij}* and equation of motion (3) for $ n \u0307 i $. In this paper,

*T*are the components of the hydrodynamic torque exerted on the particle by the fluid. In Sec. III D, we formulate the hydrodynamic torque to

_{i}*O*(Re

_{s}) via the reciprocal theorem. In Sec. III E, we perturbatively compute the resulting angular velocity to orders

*O*(Re

_{s}) and

*O*(St).

### D. Calculation of the hydrodynamic torque to order Re_{s}

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, to compute the stress tensor, and finally to 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 that we need to solve and explain how we use the reciprocal theorem to simplify the calculations.

#### 1. Navier-Stokes problem for the disturbance flow

We consider a particle with boundary $S$ immersed in a linear ambient flow (*u*^{∞}, *p*^{∞}). Throughout this paper, we express the components of the ambient flow as

or equivalently with $ \epsilon i k j \Omega k \u221e = O ij \u221e $,

Here, *S*^{∞} and *O*^{∞} are the symmetric and antisymmetric parts of the flow gradient with elements

In dimensionless variables (Sec. 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 time scale of the particle motion is the same as the time scale of the flow. The boundary condition is no-slip on the surface of the particle

We introduce the disturbance field (** u**′,

*p*′) from the particle, defined by

If we assume that (*u*^{∞}, *p*^{∞}) satisfies the Navier-Stokes equations, we have the disturbance problem

and the boundary conditions are expressed in the slip angular velocity $ \Omega i = \Omega i \u221e \u2212 \omega i $ as

Finally, when applying the reciprocal theorem, we shall use that, by definition, the divergence of the stress tensor satisfies the following equalities:

#### 2. 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, the solutions read

where

Here, *A ^{R}*,

*B*,

^{R}*C*,

^{R}*A*,

^{S}*B*,

^{S}*C*, and

^{S}*α*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.

Expressions common to both prolate and oblate spheroids | ||

$\alpha = 1 8 \lambda 2 \u2212 1 $ | $ A R = \lambda 2 \u2212 1 4 C \u2212 \lambda 3 + \lambda $ | $ B R = \lambda 2 \u2212 1 \lambda 2 + 1 4 \u2212 2 C \lambda 2 + C + \lambda 3 \u2212 \lambda $ |

$ C R = \lambda 2 \u2212 1 3 / 2 4 \u2212 2 C \lambda 2 + C + \lambda 3 \u2212 \lambda $ | $ A S = \lambda 2 \u2212 1 3 / 2 4 2 C \lambda 2 + C \u2212 3 \lambda 3 + 3 \lambda $ | $ C S = \lambda 2 \u2212 1 3 / 2 2 3 C + 2 \lambda 5 \u2212 7 \lambda 3 + 5 \lambda $ |

$ B S =\u2212 \lambda 2 \u2212 1 3 / 2 C \lambda + \lambda 4 \u2212 3 \lambda 2 + 2 8 \u2212 2 C \lambda 2 + C + \lambda 3 \u2212 \lambda \u2212 3 C \lambda + \lambda 4 + \lambda 2 \u2212 2 $ | ||

Expressions particular to prolate and oblate spheroids | ||

Oblate (λ < 1) | Prolate (λ > 1) | |

C | $\u2212 1 \u2212 \lambda 2 \u2009 cot \u2212 1 \lambda 1 \u2212 \lambda 2 $ | $ \lambda 2 \u2212 1 \u2009 coth \u2212 1 \lambda \lambda 2 \u2212 1 $ |

d | $2 1 \u2212 \lambda 2 $ | $ 2 \lambda 2 \u2212 1 \lambda $ |

c | $ i d 2 $ | $ d 2 $ |

c_{ξ} | $ 64 3 i\pi 1 \u2212 \lambda 2 3 / 2 $ | $\u2212 64 \pi \lambda 2 \u2212 1 3 / 2 3 \lambda 3 $ |

A ^{I} | $ 8 \pi \lambda 15 $ | $ 8 \pi 15 \lambda 4 $ |

B ^{I} | $ 4 \pi 15 \lambda \lambda 2 + 1 $ | $ 4 \pi \lambda 2 + 1 15 \lambda 4 $ |

Expressions common to both prolate and oblate spheroids | ||

$\alpha = 1 8 \lambda 2 \u2212 1 $ | $ A R = \lambda 2 \u2212 1 4 C \u2212 \lambda 3 + \lambda $ | $ B R = \lambda 2 \u2212 1 \lambda 2 + 1 4 \u2212 2 C \lambda 2 + C + \lambda 3 \u2212 \lambda $ |

$ C R = \lambda 2 \u2212 1 3 / 2 4 \u2212 2 C \lambda 2 + C + \lambda 3 \u2212 \lambda $ | $ A S = \lambda 2 \u2212 1 3 / 2 4 2 C \lambda 2 + C \u2212 3 \lambda 3 + 3 \lambda $ | $ C S = \lambda 2 \u2212 1 3 / 2 2 3 C + 2 \lambda 5 \u2212 7 \lambda 3 + 5 \lambda $ |

$ B S =\u2212 \lambda 2 \u2212 1 3 / 2 C \lambda + \lambda 4 \u2212 3 \lambda 2 + 2 8 \u2212 2 C \lambda 2 + C + \lambda 3 \u2212 \lambda \u2212 3 C \lambda + \lambda 4 + \lambda 2 \u2212 2 $ | ||

Expressions particular to prolate and oblate spheroids | ||

Oblate (λ < 1) | Prolate (λ > 1) | |

C | $\u2212 1 \u2212 \lambda 2 \u2009 cot \u2212 1 \lambda 1 \u2212 \lambda 2 $ | $ \lambda 2 \u2212 1 \u2009 coth \u2212 1 \lambda \lambda 2 \u2212 1 $ |

d | $2 1 \u2212 \lambda 2 $ | $ 2 \lambda 2 \u2212 1 \lambda $ |

c | $ i d 2 $ | $ d 2 $ |

c_{ξ} | $ 64 3 i\pi 1 \u2212 \lambda 2 3 / 2 $ | $\u2212 64 \pi \lambda 2 \u2212 1 3 / 2 3 \lambda 3 $ |

A ^{I} | $ 8 \pi \lambda 15 $ | $ 8 \pi 15 \lambda 4 $ |

B ^{I} | $ 4 \pi 15 \lambda \lambda 2 + 1 $ | $ 4 \pi \lambda 2 + 1 15 \lambda 4 $ |

#### 3. 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 $ ( u \u0303 i , \sigma \u0303 ij ) $ and $ ( u i \u2032 , \sigma ij \u2032 ) $ can be stated as

Here, d*F _{i}* =

*σ*d

_{ij}ξ_{j}*S*is the differential force from the fluid on the surface element with normal vector

*ξ*d

_{j}*S*. 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.

#### 4. Calculation of the torque

We choose the auxiliary problem $ ( u \u0303 i , \sigma \u0303 ij ) $ to be the Stokes flow around an identical particle rotating with an angular velocity $ \omega \u0303 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

Here, *f _{i}*(

**′) is defined in Eq. (16). We also used that $ \u2202 j \sigma \u0303 ij =0$. This equality holds because $ u \u0303 i $ is a Stokes flow. Both primary and auxiliary velocity fields vanish as |**

*u***| → ∞; 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 d**

*r**F*is the differential force exerted on the particle by the fluid. In the integrals, we identify the hydrodynamic torque on the particle. Its components are given by

_{i}It follows

The auxiliary torque $ T \u0303 j $ together with the surface integral adds up to the Jeffery torque^{1} $ T j ( 0 ) $,

evaluates to zero for any time-independent linear flow $ u i \u221e $. It follows that Eq. (21) becomes

Since $ u \u0303 i $ is linear in $ \omega \u0303 j $, this variable can be eliminated. We finally obtain

where

Thus, so 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 j ( 0 ) $, $ U \u0303 ij $, and

*f*(

_{i}**′) all depend on the direction**

*u***and upon the angular velocity**

*n***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 the 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 equation of motion (7) for the particle orientation and insert the expression for the hydrodynamic torque obtained in Sec. III D,

Now, we expand the angular velocity as

Next, we insert these expansions into 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 $ \omega i ( 0 ) $,

The dynamics of ** n** is to the lowest order given by

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

Equation (28) together with Eqs. (30) and (33) yields the effective angular velocity under the effect of weak particle and fluid inertia. From equation of motion (7), we define the effective vector field

This vector field describes the time evolution of ** n**. The first term is 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 \u221e $. In other words, they are on the form

where the tensorial coefficients $ C i j k l m ( i ) $ are composed of the remaining available tensor quantities: *n _{p}* and

*δ*(ε

_{pq}_{ijk}is already used in $ O ij \u221e $). 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

where the sum is over all 5! permutations P of (*i*, *j*, *k*, *l*, *m*) and $ \eta i ( P ) $ 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*and

_{ij}O_{jk}n_{k}*n*. Finally, we use the fact that the equation of motion may not change the magnitude of the unit vector

_{i}n_{j}O_{jk}O_{kl}n_{l}**. This constraint forces the coefficients of the two unique terms to be the same magnitude but opposite signs. Upon renaming the coefficients ±**

*n**α*

_{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 non-linear (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}.

$ S i i \u221e =0$ | Incompressible flow |

$ S ij \u221e = S ji \u221e $ | S^{∞} symmetric |

$ O ij \u221e =\u2212 O ji \u221e $ | O^{∞} anti-symmetric |

$ n i n \u0307 i =0$ | Dynamics preserves magnitude |

$ n i \u2192\u2212 n i \u27f9 n \u0307 i \u2192\u2212 n \u0307 i $ | Particle inversion symmetry |

$ S i i \u221e =0$ | Incompressible flow |

$ S ij \u221e = S ji \u221e $ | S^{∞} symmetric |

$ O ij \u221e =\u2212 O ji \u221e $ | O^{∞} anti-symmetric |

$ n i n \u0307 i =0$ | Dynamics preserves magnitude |

$ n i \u2192\u2212 n i \u27f9 n \u0307 i \u2192\u2212 n \u0307 i $ | Particle inversion symmetry |

In the particular case of a simple shear flow, we have explicitly (see Fig. 1 for the geometry)

We observe that for the simple shear, $ O ij \u221e O jk \u221e =\u2212 S ij \u221e S jk \u221e $ and $ S ij \u221e O jk \u221e =\u2212 O ij \u221e S jk \u221e $. Then, the form of the equation of motion simplifies to

with

### G. Evaluation of the volume integral in Eq. (33)

The volume integral in Eq. (33) contains four distinct terms: $ \u2202 t u i \u2032 $ represents unsteady fluid inertia and the three terms $ u j \u221e \u2202 j u i \u2032 + u j \u2032 \u2202 j u i \u221e + u j \u2032 \u2202 j u i \u2032 $ represent convective fluid inertia. We compute these four terms using explicit Stokes-flow solutions (17). While the Stokes flow has no explicit time dependence, both particle direction ** n** and angular velocity

**do. Thus, each occurrence of**

*ω**n*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

_{k}^{®}. 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 the $ e \u02c6 3 $-direction of Eq. (33) due to unsteady fluid inertia

We first perform the time derivatives in (17) in the manner explained above. Then, we insert the components of ** n** and the explicit form of 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 stems from the Stokes-flow coefficients (see Appendix A) and a spatial dependence coming from

*r*and the spheroidal integrals $ J m n $ and $ K m n $ (see Appendix B). For $n= [ 1 / 2 , 3 / 2 , 0 ] $, a typical term looks like the following:

_{i} We note that the only spatial dependence on the azimuthal angle around the symmetry axis of the body comes from factors of *r _{i}*. We introduce a rotated coordinate system in which $ r i = R ji r j \u2032 $, such that $ r 1 \u2032 $ 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 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).

## IV. RESULTS

### A. Effective equation of motion

We parametrize the vector ** n** in a spherical coordinate system (

*θ*,

*φ*) with

*θ*the polar angle and

*φ*the azimuthal angle (Fig. 1),

In these coordinates, Eq. (39) is expressed as

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 unsteady and convective fluid inertias ($ \beta \alpha ( U ) $ and $ \beta \alpha ( C ) $) 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 behaviors of all contributions to *β*_{α} in three limiting cases: thin oblate particles (*λ* → 0), thin prolate particles (*λ* → ∞), and nearly spherical particles. For nearly spherical particles, we define a small parameter *ϵ* as follows:

The asymptotic results for *β*_{α} in the limits *λ* → 0, *λ* → ∞, and |*ϵ*| → 0 are summarised in Table II and are shown as red lines (dashed-dotted) in Fig. 2.

Thin oblate particles (λ → 0)
. | ||||
---|---|---|---|---|

. | Total . | Unsteady . | Convective . | Particle . |

β_{1} | $ 11 30 $ | $ 1 5 $ | $ 1 6 $ | 0 |

β_{2} | $ 1 10 $ | $\u2212 1 20 $ | $ 3 20 $ | 0 |

β_{3} | $\u2212 1 5 $ | $\u2212 3 20 $ | $\u2212 1 20 $ | 0 |

β_{4} | $\u2212 1 3 $ | $\u2212 3 20 $ | $\u2212 11 60 $ | 0 |

Nearly spherical particles (|ϵ| ≪ 1) | ||||

Total | Unsteady | Convective | Particle | |

β_{1} | $ 137 \u03f5 2 294 $ | 0 | $ 163 \u03f5 2 490 $ | $ 2 \u03f5 2 15 $ |

β_{2} | $ 2 \u03f5 21 + 81 \u03f5 2 245 $ | $ 62 \u03f5 2 525 $ | $ \u03f5 35 + 37 \u03f5 2 294 $ | $ \u03f5 15 + 13 \u03f5 2 150 $ |

β_{3} | $\u2212 2 \u03f5 7 \u2212 229 \u03f5 2 735 $ | $\u2212 58 \u03f5 2 525 $ | $\u2212 37 \u03f5 105 \u2212 227 \u03f5 2 1470 $ | $ \u03f5 15 \u2212 7 \u03f5 2 150 $ |

β_{4} | $ 8 \u03f5 21 \u2212 103 \u03f5 2 735 $ | 0 | $ 11 \u03f5 35 \u2212 229 \u03f5 2 2450 $ | $ \u03f5 15 \u2212 7 \u03f5 2 150 $ |

Thin prolate particles (λ → ∞) | ||||

Total | Unsteady | Convective | Particle | |

β_{1} | $ 7 30 log \u2009 2 \lambda \u2212 45 $ | $ 1 8 log \u2009 2 \lambda \u2212 12 $ | $ 13 120 log \u2009 2 \lambda \u2212 180 $ | 0 |

β_{2} | $ 1 10 log \u2009 2 \lambda \u2212 15 $ | $ 1 8 log \u2009 2 \lambda \u2212 12 $ | $ 1 20 log \u2009 2 \lambda \u2212 60 $ | 0 |

β_{3} | 0 | 0 | 0 | 0 |

β_{4} | 0 | 0 | 0 | 0 |

Thin oblate particles (λ → 0)
. | ||||
---|---|---|---|---|

. | Total . | Unsteady . | Convective . | Particle . |

β_{1} | $ 11 30 $ | $ 1 5 $ | $ 1 6 $ | 0 |

β_{2} | $ 1 10 $ | $\u2212 1 20 $ | $ 3 20 $ | 0 |

β_{3} | $\u2212 1 5 $ | $\u2212 3 20 $ | $\u2212 1 20 $ | 0 |

β_{4} | $\u2212 1 3 $ | $\u2212 3 20 $ | $\u2212 11 60 $ | 0 |

Nearly spherical particles (|ϵ| ≪ 1) | ||||

Total | Unsteady | Convective | Particle | |

β_{1} | $ 137 \u03f5 2 294 $ | 0 | $ 163 \u03f5 2 490 $ | $ 2 \u03f5 2 15 $ |

β_{2} | $ 2 \u03f5 21 + 81 \u03f5 2 245 $ | $ 62 \u03f5 2 525 $ | $ \u03f5 35 + 37 \u03f5 2 294 $ | $ \u03f5 15 + 13 \u03f5 2 150 $ |

β_{3} | $\u2212 2 \u03f5 7 \u2212 229 \u03f5 2 735 $ | $\u2212 58 \u03f5 2 525 $ | $\u2212 37 \u03f5 105 \u2212 227 \u03f5 2 1470 $ | $ \u03f5 15 \u2212 7 \u03f5 2 150 $ |

β_{4} | $ 8 \u03f5 21 \u2212 103 \u03f5 2 735 $ | 0 | $ 11 \u03f5 35 \u2212 229 \u03f5 2 2450 $ | $ \u03f5 15 \u2212 7 \u03f5 2 150 $ |

Thin prolate particles (λ → ∞) | ||||

Total | Unsteady | Convective | Particle | |

β_{1} | $ 7 30 log \u2009 2 \lambda \u2212 45 $ | $ 1 8 log \u2009 2 \lambda \u2212 12 $ | $ 13 120 log \u2009 2 \lambda \u2212 180 $ | 0 |

β_{2} | $ 1 10 log \u2009 2 \lambda \u2212 15 $ | $ 1 8 log \u2009 2 \lambda \u2212 12 $ | $ 1 20 log \u2009 2 \lambda \u2212 60 $ | 0 |

β_{3} | 0 | 0 | 0 | 0 |

β_{4} | 0 | 0 | 0 | 0 |

### B. Linear stability analysis at infinitesimal Re_{s}

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 inertias 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\pi / 1 \u2212 \Lambda 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 dashed-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

### A. Effective equation of motion

Equation (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 and 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 shown in Fig. 2.

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 inertias contribute, and it would be qualitatively wrong to neglect one of these terms. This is due to the fact that the time scale of the particle motion is the same as the time scale 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.

### B. 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 exists 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, except 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, DNSs 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 and 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 small-*ϵ* limit of our calculation to the results of Ref. 17 and found 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 and 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}.

### C. 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 and 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. DNSs 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 and 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 the 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 and 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.

## Acknowledgments

This work was supported by grants from Vetenskapsrådet and the Göran Gustafsson Foundation for Research in Natural Sciences and Medicine. Support from the MPNS COST Action MP1305 “Flowing matter” is gratefully acknowledged.

### APPENDIX A: SOLUTIONS TO STOKES’ EQUATION

In this appendix, we solve the steady Stokes’ equation for an arbitrarily aligned spheroid in a general linear flow *u*^{∞} = *A*^{∞}** r**. The calculation is a special case of the calculation by Jeffery.

^{1}However, instead of the ellipsoidal harmonics that Jeffery used, we employ a finite multipole expansion, following Chwang and Wu.

^{21}The purpose of this appendix is to derive an explicit closed form expression for the Stokes flow field, suitable for evaluation in the reciprocal theorem. For a more general description of the method, we refer to the book by Kim and Karrila.

^{10}

#### 1. Formulation of the problem

Stokes’ equation reads:

with no-slip boundary conditions on the surface *S* of the particle

Here, *ω _{j}* is the angular velocity of the particle. Furthermore, it is assumed that the flow remains unperturbed at infinitely far away from the particle

We solve for the disturbance flow $ u i \u2032 = u i \u2212 u i \u221e $ that satisfies Stokes’ equation (A1) with boundary conditions,

We decompose the linear background flow $ u i \u221e $ into its symmetric and antisymmetric parts, defining the vector $ \Omega i \u221e $ and strain $ S ij \u221e $ by

Finally, in terms of the “slip angular velocity” $ \Omega i = \Omega i \u221e \u2212 \omega i $, the problem to be solved reads

#### 2. Multipoles

We solve Eq. (A6) by a finite multipole expansion.^{10,21} The multipoles are the Green’s function for the Stokes’ equation and its derivatives. In this appendix, we use the shorthand notation $ G i j , k \u2261 \u2202 k G ij $. The multipoles needed to solve for the fluid velocity field around particles in a linear flow are

The following two higher-order multipoles are required in the reciprocal theorem. We include them for reference,

Note that we use the “Oseen tensor” notation. The Green’s function for the Stokes’ equation is in fact $ G ij = G ij /8\pi $. It is convenient to split the dipole contribution $ G i j , k $ into its antisymmetric (“rotlet”) and symmetric (“stresslet”) parts. They are

#### 3. Spheroidal multipoles

Whereas the flow around a spherical particle may be represented by multipoles anchored at a single point, representing the flow around a spheroidal particle requires a weighted line distribution of multipoles.^{10,21} We therefore define the “spheroidal multipoles” as the following distributions, note especially the different weights for higher-order multipoles:

The constant *c* is related to the spheroidal geometry. Prolate and oblate coordinates are obtained by rotating an ellipse around its major or minor axis. We call the distance between the foci of the underlying ellipse *d*, then *c* = *d*/2 for prolate coordinates, and *c* = *id*/2 for oblate coordinates (see definition of coordinate systems in Appendix C.)

In order to write down explicit tensor expressions for the spheroidal multipoles, we introduce the integrals $ I m n $, $ J m n $, and $ K m n $ by

The spatial variation of the functions $ I m n $ depends upon |** r**|

^{2}and

**only. Further properties and evaluation of the integrals are discussed in Appendix B. With $ J m n $ and $ K m n $, we express the spheroidal multipoles explicitly, for example, the spheroidal rotlet,**

*r*⋅*n* The integrals $ I m n $ play the same part in spheroidal geometry as do 1/*r ^{m}* in spherical geometry. The spheroidal stresslet and quadrupole are given by

#### 4. Solution by a finite multipole expansion

The spheroidal multipoles are functions that satisfy Stokes’ equation, and a suitable linear combination of them also satisfies the no-slip boundary condition on the surface of a spheroid with symmetry axis ** n**. The remaining problem is to determine the coefficients for this linear combination.

Following Kim and Karrila,^{10} we use the following ansatz for the disturbance flow field:

where

Given the ambient strain $ S ij \u221e $ and angular slip velocity $ \Omega i = \Omega i f \u2212 \omega i p $, we must determine seven unknown scalars, which may depend upon the particle shape: *A ^{R}*,

*B*,

^{R}*C*,

^{R}*A*,

^{S}*B*,

^{S}*C*, and

^{S}*α*. When the coefficients are known, Eq. (A14) is the sought Stokes solution.

In order to match linear boundary condition Eq. (A4), we need the combinations of $ J m n $ and $ K m n $ in the ansatz to be constant on the particle surface, much like the scalar function 1/*r ^{m}* is in spherical geometry.

Upon examination, the functions $ J 3 0 $ and $ K 5 0 $ are constant on the spheroidal surface. Further, the functions $ J 3 1 $ and $ K 5 1 $ can be written as $ J 3 1 = n j r j J 3 \u2032 1 $ and $ K 5 1 = n j r j K 5 \u2032 1 $, where $ J 3 \u2032 1 $ and $ K 5 \u2032 1 $ are constant on the spheroidal surface. The remaining spheroidal functions $ J 5 n $ and $ K 7 n $ which appear in ansatz (A14) are more complicated. However, it turns out that they appear only in the combinations $ J 5 n \u221210\alpha K 7 n $. We therefore choose

With this choice of *α*, it holds that, on the surface of the spheroid,

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 \u221e $. First, with $ S jk \u221e =0$, we contract Eq. (A4) with

*n*, ε

_{i}_{ijk}

*n*Ω

_{j}_{k}, and ε

_{ipq}

*n*ε

_{q}_{pjk}

*n*Ω

_{j}_{k}. Second, with Ω

_{j}= 0, we contract Eq. (A4) with

*n*, $ S ij \u221e n j $, and finally $ \epsilon ijk n j S kl \u221e n l $. These six equations together have only one solution. We tabulate the resulting expressions for both oblate and prolate spheroids in Table III.

_{i}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 ij , k R \epsilon jkl A l $ is $ T l R =\u221216\pi 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 spheroidal rotlet (A10), we linearly superpose the contributions from all the contained rotlets. The torque from the flow $ u i = Q ij , k R \epsilon jkl B l $ is therefore

The factor *c*_{ξ} depends only on the aspect ratio of the particle (see Table III).

### APPENDIX B: SPHEROIDAL INTEGRALS

In order to solve Stokes’ equation and evaluating the volume integrals in the reciprocal theorem, we need to evaluate integrals on the form

First, when matching boundary conditions, we must evaluate the integrals with ** r** on the surface of the spheroidal particle. Second, when evaluating the reciprocal theorem, we need to integrate products of two or three $ I m n $ multiplied with the components of the spatial coordinate

**over the entire fluid volume outside the particle. Therefore, we express the functions $ I m n $ in a spheroidal coordinate system with symmetry axis $ x \u02c6 \u2032 $ along**

*r***. This is accomplished by a rotational change of variables**

*n***′ =**

*r***, $ x \u02c6 \u2032 =Rn$, where the latter equality defines a rotation**

*Rr***. The absolute value (distance) between**

*R***and**

*r**ξ*

**is preserved by a rotation and the integral is transformed into**

*n* This form is equivalent to the integrals *B*_{m,n} in Chwang and Wu.^{21} Geometrically, Eq. (B1) represents a line source along the direction ** n**. The rotation

**places the line source along the**

*R**x*′-axis in an auxiliary coordinate system. The result is a function of |

**|**

*r*^{2}and $ x \u2032 = x \u02c6 \u2032 \u22c5 r \u2032 =n\u22c5r$.

Explicit expressions for $ I m n $ may be found by direct integration or by a recursion formula.^{21} Since we require only a finite number of integrals, we simply perform the direct integration once and for all and save the result in a table.

Finally, when evaluating the term corresponding to unsteady fluid inertia in the volume integral of the reciprocal theorem, we need to compute the derivatives of $ I m n $ with respect to the moving vector ** n**. By differentiating Eq. (B1), we derive the following formula:

### APPENDIX C: SPHEROIDAL COORDINATES

Both oblate and prolate spheroidal coordinates are extensions of a two-dimensional elliptic coordinate system (*ξ*_{1}, *ξ*_{2}). The *ξ*_{1}-coordinate represents concentric ellipses, while *ξ*_{2} represents the corresponding hyperbolas. Their intersections give unique coordinates in the *x*-*y*-plane. An azimuthal angle of revolution *ϕ* denotes the extension into three dimensions.

#### 1. Oblate spheroidal coordinates

Start with the *x*-*y*-plane and place an ellipse of focal distance *d* with its minor axis along the *x*-axis. Now, revolve the ellipse by 2*π* around the *x*-axis to produce an oblate spheroid. Then, *ξ*_{1} represents concentric oblate spheroidal surfaces, *ξ*_{2} represents the corresponding hyperbolic surfaces, and we call *ϕ* the angle of revolution. The coordinate equations are

The coordinate ranges are 0 ≤ *ξ*_{1} < ∞, −1 ≤ *ξ*_{2} ≤ 1, and 0 ≤ *ϕ* ≤ 2*π*, and the volume element $dV= 1 8 d 3 \xi 1 2 + \xi 2 2 d \xi 1 d \xi 2 d\varphi $. In this paper, we treat oblate spheroids with dimensionless major axis length unity and minor axis length *λ*. These lengths determine the focal distance *d* as

and the particle surface is parameterised by

#### 2. Prolate spheroidal coordinates

Start with the *x*-*y*-plane and place an ellipse of focal distance *d* with its major axis along the *x*-axis. Now, revolve the ellipse by 2*π* around the *x*-axis to produce a prolate spheroid. Then, *ξ*_{1} represents concentric prolate spheroidal surfaces, *ξ*_{2} represents the corresponding hyperbolic surfaces, and we call *ϕ* the angle of revolution. The coordinate equations are

The coordinate ranges are 1 ≤ *ξ*_{1} < ∞, −1 ≤ *ξ*_{2} ≤ 1, and 0 ≤ *ϕ* ≤ 2*π*, and the volume element $dV= 1 8 d 3 \xi 1 2 \u2212 \xi 2 2 d \xi 1 d \xi 2 d\varphi $. In this paper, we treat prolate spheroids with dimensionless major axis length unity and minor axis length 1/*λ*. These lengths determine the focal distance *d* as

and the particle surface is parameterised by