---
abstract: 'The orientational dynamics of weakly inertial axisymmetric particles in a steady flow is investigated. We derive an asymptotic equation of motion for the unit axial vector along the particle symmetry axis, valid for small Stokes number ${\ensuremath{\mbox{St}}}$, and for any axisymmetric particle in any steady linear viscous flow. This reduced dynamics is analysed in two ways, both pertain to the case of a simple shear flow. In this case inertia induces a coupling between precession and nutation. This coupling affects the dynamics of the particle, breaks the degeneracy of the Jeffery orbits, and creates two limiting periodic orbits. We calculate the leading-order Floquet exponents of the limiting periodic orbits and show analytically that prolate objects tend to a tumbling orbit, while oblate objects tend to a log-rolling orbit, in agreement with previous analytical and numerical results. Second, we analyse the role of the limiting orbits when rotational noise is present. We formulate the Fokker-Planck equation describing the orientational distribution of an axisymmetric particle, valid for small ${\ensuremath{\mbox{St}}}$ and general Péclet number ${\ensuremath{\mbox{Pe}}}$. Numerical solutions of the Fokker-Planck equation, obtained by means of expansion in spherical harmonics, show that stationary orientational distributions are close to the inertia-free case when ${\ensuremath{\mbox{Pe}}}{\ensuremath{\mbox{St}}}\ll1$, whereas they are determined by inertial effects, though small, when ${\ensuremath{\mbox{Pe}}}\gg 1/{\ensuremath{\mbox{St}}}\gg1$.'
address:
- 'Department of Physics, University of Gothenburg, SE-41296 Gothenburg, Sweden'
- 'Department of Mathematics and Mechanics, LUSAC-ESIX, University of Caen, France'
author:
- 'J. Einarsson'
- 'J. R. Angilella'
- 'B. Mehlig'
title: Orientational dynamics of weakly inertial axisymmetric particles in steady viscous flows
---

shear flow, Jeffery orbits, particle inertia, Floquet analysis, Fokker-Planck equation

Introduction
============

Suspensions of rigid particles are abundant in both nature and technology, and are consequently studied in many disciplines of science. Examples are ash and ice crystals in the atmosphere [@Rolf2012], plankton in the ocean [@Stommel1949; @Gua12], or fibres in paper pulp [@Lun11]. In many situations the concentration of the suspended particles in the carrying fluid is small, so that particles do not affect the underlying fluid flow. In this case, considering the fluid flow field as given, a question of fundamental interest is to determine the resulting trajectory and orientation of a particle. The orientational dynamics of axisymmetric particles is often described by the Jeffery equation of motion [@Jeff22], the rotational analogue to simple center-of-mass advection along streamlines. This equation of motion is valid in the limit of vanishing Stokes number ${\ensuremath{\mbox{St}}}$ and particle Reynolds number ${\ensuremath{\mbox{Re}}}$ (both numbers are defined below), representing particle and fluid inertia respectively. The Jeffery equation of motion is widely used to predict mechanical [@Bre74; @Lun11; @pet99] as well as optical properties [@Gauthier1998; @wil09; @wil11; @Bez10] of dilute suspensions of inertia-free rod-like or disk-like particles. It is also used to describe the orientational motion of single particles [@Gua12; @ein12], and the alignment and tumbling of axisymmetric particles in turbulence [@Pum11; @Par12].

In the case of simple shear flow an inertia-free particle rotates in one of infinitely many possible Jeffery orbits. In other words, the solution of the Jeffery equation is degenerate, forming a one-parameter family of closed orbits. Which orbit is followed depends upon the initial condition. By taking into account the inertia of the particle in numerical simulations, it was observed that attracting orbits exist [@Lundell2010; @Lundell2011], where particles tend to align their small axis towards the vorticity axis. In the special case of nearly spherical particles the existence of attracting orbits has been demonstrated analytically [@Sub06]. The significance of this inertial effect is that it breaks the degeneracy of the Jeffery orbits.

If, in addition, particles are subject to a random force, the resulting orientational distribution is expected to be determined by a balance between the dispersing random forces and the aligning inertial forces.

The goal of the present study is to investigate the orientational dynamics of axisymmetric particles with small inertia, to describe and quantify the attracting orbits, and to analyse their effect on the orientational distribution when noise is present.

To accomplish this goal we first derive an equation of motion for the orientation of an axisymmetric particle in steady flows in the absence of noise, valid for small values of the Stokes number ${\ensuremath{\mbox{St}}}$. We define the Stokes number by ${\ensuremath{\mbox{St}}}= ms/b\mu$, where $m$ is the particle mass, $s$ is a typical flow-gradient rate, $b$ is the minor particle length and $\mu$ the dynamic viscosity of the fluid (see Sec. 2 for details). We make use of the fact that ${\ensuremath{\mbox{St}}}$ is small, so that the angular velocity of the particle is a small perturbation of the angular velocity in the absence of particle inertia. Our results, valid for any aspect ratio, generalise a result for nearly spherical spheroids reported in Ref. [@Sub06]. In Ref. [@Sub06] fluid inertia, as measured by the particle Reynolds number, is also considered. The particle Reynolds number is defined by ${\ensuremath{\mbox{Re}}}= (\kappa\rho_f/\rho_p){\ensuremath{\mbox{St}}}$, where $\rho_f$ and $\rho_p$ denote fluid and particle density. The shape factor $\kappa$ is $\kappa=1$ for a disk-shaped particle, and $\kappa=\lambda$ for a rod-like particle ($\lambda$ is the aspect ratio of the particle). Finite-${\ensuremath{\mbox{Re}}}$ corrections are important for neutrally buoyant particles \[$\rho_p/\rho_f=\mathcal O(1)$\], because the fluid inertia induces $\mathcal O({\ensuremath{\mbox{Re}}})$-corrections in the hydrodynamic torque experienced by the particle [@Sub06; @lov93]. The results obtained here are valid for heavy particles such that $\rho_p \gg \kappa\rho_f$.

The new equation of motion allows us to derive a Fokker-Planck equation determining the combined effects of particle inertia and noise upon the orientational distribution of the particles. This equation describes the evolution of an ensemble of weakly inertial particles in a steady flow, subject to random Brownian rotations. It is solved in the case of a simple shear flow by means of expansion in spherical harmonics. We find that the stationary orientational distributions differ significantly from the inertia-free case.

The rest of this paper is organised as follows. Asymptotic equations of motion valid for small ${\ensuremath{\mbox{St}}}$ in the absence of noise are obtained in Sec. 2. Then, in Sec. 3, we investigate the case of a simple shear flow. We evaluate the convergence towards one of two possible periodic orbits by calculating the corresponding Floquet exponents. Finally, in Sec. 4, we formulate the Fokker-Planck equation describing the combined effect of noise and weak particle inertia. This equation is then solved numerically and the stationary distribution is discussed. Sec. \[sec:conc\] summarises our conclusions.

Equation of motion
==================

The aim of this Section is to derive an equation of motion for the orientation of a small axisymmetric particle, valid to first order in ${\ensuremath{\mbox{St}}}$. The orientation is described by a unit vector ${\ensuremath{\mbox{\boldmath$n$}}}$ pointing in the direction of the symmetry axis of the particle. In the following we derive an equation of the form $\dot{{\ensuremath{\mbox{\boldmath$n$}}}}={\ensuremath{\mbox{\boldmath$F$}}}_0({\ensuremath{\mbox{\boldmath$n$}}}) + {\ensuremath{\mbox{St}}}{\ensuremath{\mbox{\boldmath$F$}}}_1({\ensuremath{\mbox{\boldmath$n$}}}) + \mathcal O ({\ensuremath{\mbox{St}}}^2)$, where the dot denotes the time derivative, and the functions $F_i$ depend on the flow gradient ${\ensuremath{\mathbb{A}}} = {\ensuremath{\mbox{\boldmath$\nabla$}}} {\ensuremath{\mbox{\boldmath$u$}}}$. We decompose ${\ensuremath{\mathbb{A}}}$ into its symmetric and antisymmetric parts $$\begin{aligned}
    &{\ensuremath{\mathbb{O}}} = \frac{1}{2}({\ensuremath{\mathbb{A}}} - {\ensuremath{\mathbb{A}}}{^{\mathrm T}}),\quad
    {\ensuremath{\mathbb{S}}} = \frac{1}{2}({\ensuremath{\mathbb{A}}} + {\ensuremath{\mathbb{A}}}{^{\mathrm T}}),\quad
    {\ensuremath{\mathbb{A}}} = {\ensuremath{\mathbb{O}}} + {\ensuremath{\mathbb{S}}}.\end{aligned}$$ The matrix ${\ensuremath{\mathbb{O}}}$ is related to the vorticity $\nabla {\wedge}{\ensuremath{\mbox{\boldmath$u$}}}$ of the flow: ${\ensuremath{\mathbb{O}}} {\ensuremath{\mbox{\boldmath$x$}}} = {\ensuremath{\mbox{\boldmath$\Omega$}}} {\wedge}{\ensuremath{\mbox{\boldmath$x$}}}$, where ${\ensuremath{\mbox{\boldmath$\Omega$}}} = \nabla {\wedge}{\ensuremath{\mbox{\boldmath$u$}}}/2$.

We start with Newton’s equation of motion for the orientational degrees of freedom for a rigid body. Let ${\ensuremath{\mbox{\boldmath$\omega$}}}$ denote the angular velocity of the particle, then $$\begin{aligned}
    \dot{{\ensuremath{\mbox{\boldmath$\omega$}}}} &= {\ensuremath{\mathbb{M}}}^{-1}\left(-\dot{{\ensuremath{\mathbb{M}}}}{\ensuremath{\mbox{\boldmath$\omega$}}} + {\ensuremath{\mbox{\boldmath$T$}}}\right) {\nonumber}\\
    \dot{{\ensuremath{\mbox{\boldmath$n$}}}} &= {\ensuremath{\mbox{\boldmath$\omega$}}} {\wedge}{\ensuremath{\mbox{\boldmath$n$}}}{\label{eqn:fulleqns}}\end{aligned}$$ where ${\ensuremath{\mathbb{M}}}$ is the moment of inertia, and ${\ensuremath{\mbox{\boldmath$T$}}}$ the torque applied to the body. The hydrodynamic torque on a particle depends on the flow configuration at the particle position. Therefore, when studying inhomogenous flows such as turbulent channel flows, [Eq. (\[eqn:fulleqns\])]{} must be completed with equations for the particle center-of-mass motion [@Mar10]. Here we consider homogeneous flows where the center-of-mass motion may be neglected. With the hydrodynamic torque calculated by Jeffery [@Jeff22], [Eq. (\[eqn:fulleqns\])]{} describes the rotation of both axisymmetric [@Lundell2010] and asymmetric [@Lundell2011] particles. In the following we consider the case of axisymmetric particles.

An axisymmetric object has two principal moments of inertia, around and transverse to the axis of symmetry. The moment of inertia tensor is therefore on the form $$\begin{aligned}
    {\ensuremath{\mathbb{M}}} &= \widetilde{X}^I{{\ensuremath{\mbox{\boldmath$n$}}} {\ensuremath{\mbox{\boldmath$n$}}}{^{\mathrm T}}}+ \widetilde{Y}^I ({{\ensuremath{\mathbb{I}}}}- {{\ensuremath{\mbox{\boldmath$n$}}} {\ensuremath{\mbox{\boldmath$n$}}}{^{\mathrm T}}}),{\label{eqn:momentofinertia}}\end{aligned}$$ where ${{\ensuremath{\mathbb{I}}}}$ is the identity matrix, and $\widetilde{X}^I, \widetilde{Y}^I$ are constants depending on the shape of the particle. Here and in the following our notation is similar to that used in the book by Kim and Karrila [@Kim91]. The difference, and the reason we keep the tildes, is that here the quantities have dimensions. For instance a spheroid with lengths $a$ and $b$ and mass $m$ has $\widetilde{X}^I=2mb^2/5$ and $\widetilde{Y}^I=m(a^2+b^2)/5$. Further, the hydrodynamic torque depends linearly on the flow gradient via resistance tensors which are determined solely by the shape of the particle [@Bre74; @Jeff22; @Kim91] $$\begin{aligned}
    {\ensuremath{\mbox{\boldmath$T$}}} &= \left[\widetilde{X}^C {{\ensuremath{\mbox{\boldmath$n$}}} {\ensuremath{\mbox{\boldmath$n$}}}{^{\mathrm T}}}+ \widetilde{Y}^C({{\ensuremath{\mathbb{I}}}}- {{\ensuremath{\mbox{\boldmath$n$}}} {\ensuremath{\mbox{\boldmath$n$}}}{^{\mathrm T}}})\right]({\ensuremath{\mbox{\boldmath$\Omega$}}}-{\ensuremath{\mbox{\boldmath$\omega$}}})+\widetilde{Y}^H {\ensuremath{\mbox{\boldmath$n$}}} {\wedge}{\ensuremath{\mathbb{S}}} {\ensuremath{\mbox{\boldmath$n$}}}.{\label{eqn:torque}}\end{aligned}$$ Here $\widetilde{X}^C, \widetilde{Y}^C$ and $\widetilde{Y}^H$ are shape constants. The relation to the dimensionless shape functions given in [@Kim91] is, for instance, $\widetilde{X}^C=8\pi\mu a^3 X^C$. The shape functions for spheroids are given on p.64 (prolate case) and p.68 (oblate case) of Ref. [@Kim91]. Inserting the explicit forms [Eq. (\[eqn:momentofinertia\])]{} and [Eq. (\[eqn:torque\])]{} into [Eq. (\[eqn:fulleqns\])]{} yields $$\begin{aligned}
    \dot{{\ensuremath{\mbox{\boldmath$\omega$}}}} &= \left[\frac{\widetilde{X}^C}{\widetilde{X}^I}{{\ensuremath{\mbox{\boldmath$n$}}} {\ensuremath{\mbox{\boldmath$n$}}}{^{\mathrm T}}}+ \frac{\widetilde{Y}^C}{\widetilde{Y}^I}({{\ensuremath{\mathbb{I}}}}-{{\ensuremath{\mbox{\boldmath$n$}}} {\ensuremath{\mbox{\boldmath$n$}}}{^{\mathrm T}}})\right]({\ensuremath{\mbox{\boldmath$\Omega$}}} - {\ensuremath{\mbox{\boldmath$\omega$}}}) {\nonumber}\\
    &\quad+ \frac{\widetilde{Y}^H}{\widetilde{Y}^I}{\ensuremath{\mbox{\boldmath$n$}}}{\wedge}{\ensuremath{\mathbb{S}}} {\ensuremath{\mbox{\boldmath$n$}}} -\frac{\widetilde{X}^I-\widetilde{Y}^I}{\widetilde{Y}^I}({\ensuremath{\mbox{\boldmath$\omega$}}} {\wedge}{\ensuremath{\mbox{\boldmath$n$}}}){\ensuremath{\mbox{\boldmath$n$}}}{^{\mathrm T}}{\ensuremath{\mbox{\boldmath$\omega$}}}.{\label{eqn:dimomega}}\end{aligned}$$ An axisymmetric body has two different lengths $a > b$, and we choose to base the Stokes number on the lesser of the two, such that ${\ensuremath{\mbox{St}}}=ms/b\mu$. Here $m$ is the mass of the particle, $s$ is a typical shear rate and $\mu$ is the dynamic viscosity of the surrounding fluid. We introduce dimensionless variables $t=t'/s$, ${\ensuremath{\mathbb{A}}} = s{\ensuremath{\mathbb{A}}}'$, ${\ensuremath{\mbox{\boldmath$\omega$}}} = s {\ensuremath{\mbox{\boldmath$\omega$}}}'$, and drop the primes.

We now assume that for small ${\ensuremath{\mbox{St}}}$ the angular velocity ${\ensuremath{\mbox{\boldmath$\omega$}}}$ is almost the same as in the unperturbed case: $$\begin{aligned}
    {\ensuremath{\mbox{\boldmath$\omega$}}} = {\ensuremath{\mbox{\boldmath$\omega$}}}_0 + {\ensuremath{\mbox{St}}}{\ensuremath{\mbox{\boldmath$\omega$}}}_1 + \mathcal{O}({\ensuremath{\mbox{St}}}^2).\end{aligned}$$ Inserting into [Eq. (\[eqn:dimomega\])]{} and comparing order by order in ${\ensuremath{\mbox{St}}}$ we find $$\begin{aligned}
    0 &= \left[\alpha^X {{\ensuremath{\mbox{\boldmath$n$}}} {\ensuremath{\mbox{\boldmath$n$}}}{^{\mathrm T}}}+ \alpha^Y({{\ensuremath{\mathbb{I}}}}-{{\ensuremath{\mbox{\boldmath$n$}}} {\ensuremath{\mbox{\boldmath$n$}}}{^{\mathrm T}}}) \right]({\ensuremath{\mbox{\boldmath$\Omega$}}} - {\ensuremath{\mbox{\boldmath$\omega$}}}_0) + \alpha^S {\ensuremath{\mbox{\boldmath$n$}}} {\wedge}{\ensuremath{\mathbb{S}}} {\ensuremath{\mbox{\boldmath$n$}}}, {\label{eqn:order0}}\\
    \dot{{\ensuremath{\mbox{\boldmath$\omega$}}}}_0 &= -\left[\alpha^X {{\ensuremath{\mbox{\boldmath$n$}}} {\ensuremath{\mbox{\boldmath$n$}}}{^{\mathrm T}}}+ \alpha^Y({{\ensuremath{\mathbb{I}}}}-{{\ensuremath{\mbox{\boldmath$n$}}} {\ensuremath{\mbox{\boldmath$n$}}}{^{\mathrm T}}}) \right]{\ensuremath{\mbox{\boldmath$\omega$}}}_1 + \alpha^I ({\ensuremath{\mbox{\boldmath$\omega$}}}_0{\wedge}{\ensuremath{\mbox{\boldmath$n$}}}){\ensuremath{\mbox{\boldmath$n$}}}{^{\mathrm T}}{\ensuremath{\mbox{\boldmath$\omega$}}}_0 {\label{eqn:order1}}.\end{aligned}$$ Here the dimensionless coefficients $\alpha^i$, $i=X$, $Y$, $S$, $I$ are defined by $$\begin{aligned}
    &\frac{s}{{\ensuremath{\mbox{St}}}}\alpha^X \equiv \frac{\widetilde{X}^C}{\widetilde{X}^I},
    \quad \frac{s}{{\ensuremath{\mbox{St}}}}\alpha^Y \equiv \frac{\widetilde{Y}^C}{\widetilde{Y}^I}, 
    \quad \frac{s}{{\ensuremath{\mbox{St}}}}\alpha^S \equiv \frac{\widetilde{Y}^H}{\widetilde{Y}^I},
    \quad \alpha^I \equiv \frac{\widetilde{Y}^I-\widetilde{X}^I}{\widetilde{Y}^I}.\end{aligned}$$ The lowest order, ${\ensuremath{\mbox{\boldmath$\omega$}}}_0$ is found from [Eq. (\[eqn:order0\])]{} to be $$\begin{aligned}
    {\ensuremath{\mbox{\boldmath$\omega$}}}_0 &= {\ensuremath{\mbox{\boldmath$\Omega$}}} + \frac{\alpha^S}{\alpha^Y} {\ensuremath{\mbox{\boldmath$n$}}} {\wedge}{\ensuremath{\mathbb{S}}} {\ensuremath{\mbox{\boldmath$n$}}}\end{aligned}$$ which we recognise as the Jeffery angular velocity [@Kim91]. In particular, for a spheroid of aspect ratio $\lambda=a/b$ the material constant $\alpha^S/\alpha^Y=(\lambda^2-1)/(\lambda^2+1)\equiv \Lambda$. To recover the well-known Jeffery equation for $\dot{{\ensuremath{\mbox{\boldmath$n$}}}}$, note that $$\begin{aligned}
    \dot{{\ensuremath{\mbox{\boldmath$n$}}}} &= {\ensuremath{\mbox{\boldmath$\omega$}}}_0 {\wedge}{\ensuremath{\mbox{\boldmath$n$}}} + \mathcal O ({\ensuremath{\mbox{St}}})
    = {\ensuremath{\mathbb{O}}} {\ensuremath{\mbox{\boldmath$n$}}} + \frac{\alpha^S}{\alpha^Y}({\ensuremath{\mathbb{S}}} {\ensuremath{\mbox{\boldmath$n$}}} - {\ensuremath{\mbox{\boldmath$n$}}} ({\ensuremath{\mbox{\boldmath$n$}}}{^{\mathrm T}}{\ensuremath{\mathbb{S}}} {\ensuremath{\mbox{\boldmath$n$}}}))  + \mathcal O ({\ensuremath{\mbox{St}}}).\end{aligned}$$ Solving for the next order ${\ensuremath{\mbox{\boldmath$\omega$}}}_1$ in [Eq. (\[eqn:order1\])]{} we find $$\begin{aligned}
    {\ensuremath{\mbox{\boldmath$\omega$}}}_1 &= \frac{\alpha^I}{\alpha^Y}({\ensuremath{\mbox{\boldmath$\omega$}}}_0 {\wedge}{\ensuremath{\mbox{\boldmath$n$}}}){\ensuremath{\mbox{\boldmath$n$}}}{^{\mathrm T}}{\ensuremath{\mbox{\boldmath$\omega$}}}_0 {\nonumber}\\ 
    &\quad- \left[\frac{\alpha^S}{\alpha^Y \alpha^X}{{\ensuremath{\mbox{\boldmath$n$}}} {\ensuremath{\mbox{\boldmath$n$}}}{^{\mathrm T}}}+ \frac{\alpha^S}{\alpha^Y \alpha^Y}({{\ensuremath{\mathbb{I}}}}-{{\ensuremath{\mbox{\boldmath$n$}}} {\ensuremath{\mbox{\boldmath$n$}}}{^{\mathrm T}}})\right]\left(({\ensuremath{\mbox{\boldmath$\omega$}}}_0 {\wedge}{\ensuremath{\mbox{\boldmath$n$}}}){\wedge}{\ensuremath{\mathbb{S}}} {\ensuremath{\mbox{\boldmath$n$}}} \right) {\nonumber}\\
    &\quad- \frac{\alpha^S}{\alpha^Y \alpha^Y}{\ensuremath{\mbox{\boldmath$n$}}} {\wedge}{\ensuremath{\mathbb{S}}}({\ensuremath{\mbox{\boldmath$\omega$}}}_0 {\wedge}{\ensuremath{\mbox{\boldmath$n$}}}).{\label{eqn:w1general}}\end{aligned}$$ This expression for the correction to the angular velocity is valid for any axisymmetric body, specified by the parameters $\alpha^i$. For the case of a simple shear flow and nearly spherical spheroids [Eq. (\[eqn:w1general\])]{} agrees with earlier results (Eq. (3.12) in Ref. [@Sub06] is equal to [Eq. (\[eqn:w1general\])]{} to first order in $\Lambda$). In the case of spheroids, both prolate and oblate, we insert the appropriate shape functions from [@Kim91] and observe that $$\begin{aligned}
    &\alpha^I = \frac{\alpha^S}{\alpha^Y} = \frac{\lambda^2-1}{\lambda^2+1}\equiv \Lambda.\end{aligned}$$ We find our final result for the equation of motion for ${\ensuremath{\mbox{\boldmath$n$}}}$ to order $\mathcal O({\ensuremath{\mbox{St}}})$ as $$\begin{aligned}
    \dot{{\ensuremath{\mbox{\boldmath$n$}}}} &= {\ensuremath{\mathbb{O}}} {\ensuremath{\mbox{\boldmath$n$}}} + \Lambda({{\ensuremath{\mathbb{S}}}{\ensuremath{\mbox{\boldmath$n$}}} - {\ensuremath{\mbox{\boldmath$n$}}}({\ensuremath{\mbox{\boldmath$n$}}}{^{\mathrm T}}{\ensuremath{\mathbb{S}}}{\ensuremath{\mbox{\boldmath$n$}}})}) {\nonumber}\\
    &+{\nonumber}\frac{{\ensuremath{\mbox{St}}}\Lambda}{\alpha^{Y}}\bigg[-\left( {({\ensuremath{\mathbb{O}}}{\ensuremath{\mathbb{O}}}+{\ensuremath{\mathbb{S}}}{\ensuremath{\mathbb{O}}}+\Lambda{\ensuremath{\mathbb{S}}}{\ensuremath{\mathbb{S}}}){\ensuremath{\mbox{\boldmath$n$}}} - {\ensuremath{\mbox{\boldmath$n$}}}({\ensuremath{\mbox{\boldmath$n$}}}{^{\mathrm T}}({\ensuremath{\mathbb{O}}}{\ensuremath{\mathbb{O}}}+{\ensuremath{\mathbb{S}}}{\ensuremath{\mathbb{O}}}+\Lambda{\ensuremath{\mathbb{S}}}{\ensuremath{\mathbb{S}}}){\ensuremath{\mbox{\boldmath$n$}}})} \right)\\{\nonumber}&\qquad\qquad-\Lambda{\ensuremath{\mbox{\boldmath$n$}}}{^{\mathrm T}}{\ensuremath{\mbox{\boldmath$\Omega$}}}\left({\ensuremath{\mbox{\boldmath$n$}}} {\wedge}{\ensuremath{\mathbb{S}}} {\ensuremath{\mbox{\boldmath$n$}}}\right) \\{\nonumber}&\qquad\qquad+({\ensuremath{\mbox{\boldmath$n$}}}{^{\mathrm T}}{\ensuremath{\mathbb{S}}} {\ensuremath{\mbox{\boldmath$n$}}}){\ensuremath{\mathbb{O}}} {\ensuremath{\mbox{\boldmath$n$}}} \\
    &\qquad\qquad+2 \Lambda ({\ensuremath{\mbox{\boldmath$n$}}}{^{\mathrm T}}{\ensuremath{\mathbb{S}}} {\ensuremath{\mbox{\boldmath$n$}}})({{\ensuremath{\mathbb{S}}}{\ensuremath{\mbox{\boldmath$n$}}} - {\ensuremath{\mbox{\boldmath$n$}}}({\ensuremath{\mbox{\boldmath$n$}}}{^{\mathrm T}}{\ensuremath{\mathbb{S}}}{\ensuremath{\mbox{\boldmath$n$}}})})\big] + \mathcal O({\ensuremath{\mbox{St}}}^2)\,. {\label{eqn:nequation}}\end{aligned}$$ [Eq. (\[eqn:nequation\])]{} is valid as an approximation to [Eq. (\[eqn:fulleqns\])]{} when particles are weakly inertial. The vector ${\ensuremath{\mbox{\boldmath$n$}}}$ is advected in an effective vector field ${\ensuremath{\mbox{\boldmath$F$}}}_0({\ensuremath{\mbox{\boldmath$n$}}}) + {\ensuremath{\mbox{St}}}{\ensuremath{\mbox{\boldmath$F$}}}_1({\ensuremath{\mbox{\boldmath$n$}}})$ that differs from Jeffery’s equation by a term linear in ${\ensuremath{\mbox{St}}}$. This approximation does not capture the formation of caustics that occur when particles with different orientational histories arrive at the same orientation but with different angular velocities [@gus13]. From [Eq. (\[eqn:nequation\])]{} we see that the geometrical factor multiplying the inertial correction is $1/\alpha^Y$. Therefore, in the following, we use the effective Stokes number ${\ensuremath{\widetilde {\ensuremath{\mbox{St}}}}}= {\ensuremath{\mbox{St}}}/\alpha^Y$. Solutions to the full equation of motion [Eq. (\[eqn:fulleqns\])]{} are compared to solutions of [Eq. (\[eqn:nequation\])]{} for a simple shear flow and varying $\lambda$ and ${\ensuremath{\widetilde {\ensuremath{\mbox{St}}}}}$. We observe that the approximation works well for the examples shown in Fig. \[fig:eqcheck\], although the deviation in panel (d) is discernible.

![\[fig:eqcheck\] Numerical solutions of [Eq. (\[eqn:fulleqns\])]{} (red), and [Eq. (\[eqn:nequation\])]{} (blue) for four different particles in a simple shear flow. (a) $\lambda=1/5$, ${\ensuremath{\widetilde {\ensuremath{\mbox{St}}}}}=0.1$; (b) $\lambda=5$, ${\ensuremath{\widetilde {\ensuremath{\mbox{St}}}}}=0.25$; (c) $\lambda=1/5$, ${\ensuremath{\widetilde {\ensuremath{\mbox{St}}}}}=0.25$; (d) $\lambda=5$, ${\ensuremath{\widetilde {\ensuremath{\mbox{St}}}}}=0.7$. Trajectories on the sphere are shown projected onto a disk with the area-preserving Lambert azimuthal stereographic projection $X=n_x\sqrt{2/(1+n_z)}$, $Y=n_y\sqrt{2/(1+n_z)}$ [@sny87]. The grey boundary marks the equator, the blue dot is the initial condition.](fig/comparison1.pdf)

Simple shear flow: stability analysis
=====================================

The simple shear flow, ${\ensuremath{\mbox{\boldmath$u$}}}({\ensuremath{\mbox{\boldmath$r$}}}, t)= s y \hat{{\ensuremath{\mbox{\boldmath$x$}}}}$, is interesting because it is the local linearisation of any flow with parallel streamlines, like pipe flows. It is also the flow in an ideal Couette rheometer used to measure the viscosity of fluids. The solutions of [Eq. (\[eqn:nequation\])]{} in the simple shear at ${\ensuremath{\mbox{St}}}=0$ are a one-parameter family of closed periodic orbits, the well-known Jeffery orbits [@Jeff22]. The effect of the $\mathcal O({\ensuremath{\mbox{St}}})$ terms in [Eq. (\[eqn:nequation\])]{} is a drift towards a stable periodic orbit. In this section we will characterise this *orbit drift* for general particle shapes $\lambda$ by computing the Floquet exponents of the limiting orbits.

![\[fig:coordinates\] Simple shear flow and coordinate system.](fig/coordinates.png)

The modified vector field on the right hand side of [Eq. (\[eqn:nequation\])]{} allows for two limiting orbits. First the equatorial orbit, where the ${\ensuremath{\mbox{\boldmath$n$}}}$-vector moves in the plane perpendicular to vorticity (usually referred to as *tumbling*), and second the polar orbit, where the ${\ensuremath{\mbox{\boldmath$n$}}}$-vector is fixed parallel to the vorticity (referred to as *log rolling*). We introduce a spherical coordinate system where $\theta$ is the polar angle, measured from the vorticity vector, and $\varphi$ is the azimuthal angle, measured from the flow direction, see Fig. \[fig:coordinates\]. Then the two limiting orbits correspond to $\theta=0$ and $\theta = \pi/2$. In spherical coordinates, [Eq. (\[eqn:nequation\])]{} takes the following form for a simple shear flow: $$\begin{aligned}
    \dot{\varphi} &= f(\theta, \varphi) = \frac{1}{2} (-1+\Lambda  \cos 2 \varphi )  \\
    &\quad + \frac{{\ensuremath{\widetilde {\ensuremath{\mbox{St}}}}}\Lambda}{8} 
     \left(2 \Lambda  \sin ^2\theta  \sin 4 \varphi +\sin 2 \varphi  ((\Lambda +1) \cos 2
   \theta +\Lambda -3)\right)\,,\\
    \dot{\theta} &= g(\theta, \varphi) = \Lambda  \sin \theta  \cos \theta  \sin \varphi  \cos \varphi   \\ 
    &\quad +  \frac{{\ensuremath{\widetilde {\ensuremath{\mbox{St}}}}}\Lambda}{2} \sin \theta  \cos \theta  \cos ^2(\varphi ) \left(4 \Lambda  \sin ^2\theta  \sin
   ^2\varphi -\Lambda +1\right)\,.\end{aligned}$$ These equations show that inertia induces a coupling between nutation and precession, in contrast with the inertia-free case where the $\varphi$-dynamics is independent of $\theta$. In the following we denote the limiting periodic orbits $\varphi_0(t)$ and $\theta_0=0$ or $\pi/2$. The solutions $\varphi_0(t)$ are in both cases periodic with the Jeffery period $T=4\pi/\sqrt{1-\Lambda^2}$. If we introduce arbitrarily small perturbations $\theta(t)=\theta_0+\delta \theta(t)$ and $\varphi(t)=\varphi_0(t)+\delta \varphi(t)$, then the perturbation $\delta \theta(t)$ evolves according to $$\begin{aligned}
    \frac{{\ensuremath{\mathrm{d}}}}{{\ensuremath{\mathrm{d}}}t} \delta \theta &=  \delta \theta \frac{\partial g}{\partial \theta}{\,\bigg|_{(\varphi_0(t), \theta_0)}} \!\!\!\!\!+\, \delta \varphi \frac{\partial g}{\partial \varphi}{\,\bigg|_{(\varphi_0(t), \theta_0)}} \,.\end{aligned}$$ We note that ${\partial g}/{\partial \varphi} = 0$ when evaluated at $\theta=0$ or $\theta=\pi/2$. This simplification yields $$\begin{aligned}
    \frac{{\ensuremath{\mathrm{d}}}}{{\ensuremath{\mathrm{d}}}t}\ln \delta \theta &= \frac{\partial g}{\partial \theta}{\,\bigg|_{(\varphi_0(t), \theta_0)}}{\label{eqn:lyap1}}\,.\end{aligned}$$ The Floquet exponent is obtained by time-averaging [Eq. (\[eqn:lyap1\])]{} along the unperturbed periodic orbit $(\varphi_0, \theta_0)$: $$\begin{aligned}
    \gamma_\theta &= \frac{1}{T}\int_0^{T}{\ensuremath{\mathrm{d}}}t \frac{\partial g}{\partial \theta}{\,\bigg|_{(\varphi_0(t), \theta_0)}} = -\frac{1}{T}\int_0^{2\pi}\frac{{\ensuremath{\mathrm{d}}}\varphi}{\dot \varphi} \frac{\partial g}{\partial \theta}{\,\bigg|_{(\varphi, \theta_0)}}\,.{\label{eqn:floquet}}\end{aligned}$$ The last equality follows from the fact that in one period, $\varphi$ travels monotonically through $\varphi=0$ to $-2\pi$. Upon evaluating [Eq. (\[eqn:floquet\])]{} we obtain for the polar, log rolling orbit $$\begin{aligned}
    \gamma_0 &= \frac{{\ensuremath{\widetilde {\ensuremath{\mbox{St}}}}}}{4}\Lambda(1-\Lambda) + \mathcal O({\ensuremath{\mbox{St}}}^2),\end{aligned}$$ and for the equatorial, tumbling orbit $$\begin{aligned}
    \gamma_{\pi/2} &= \frac{{\ensuremath{\widetilde {\ensuremath{\mbox{St}}}}}}{4}(1-\Lambda)(\sqrt{1-\Lambda^2}-\Lambda-1)+ \mathcal O({\ensuremath{\mbox{St}}}^2).\end{aligned}$$ A positive Floquet exponent $\gamma$ indicates an unstable orbit, and a negative an attracting, stable orbit. The exponents (divided by ${\ensuremath{\widetilde {\ensuremath{\mbox{St}}}}}$) are shown as a function of particle shape $\lambda$ in Fig. \[fig:lyapunovs\]. We observe that the stable limiting orbit for flat, disk-shaped particles is the polar, log rolling orbit, while rod-shaped particles are attracted to the equatorial, tumbling orbit. This is consistent with the intuition that the most stable rotation is that with the mass distributed far from the axis of rotation. Our result explains the orbit drift observed in [@Lundell2010] at small values of [$\mbox{St}$]{}. In Ref. [@Sub06] this drift (at small Stokes number and vanishing particle Reynolds number) was explained for nearly spherical particles. Here we generalise the result to arbitrary axisymmetric particles.

![\[fig:lyapunovs\] Floquet exponents $\gamma$ (divided by [$\widetilde {\ensuremath{\mbox{St}}}$]{}) as functions of aspect ratio $\lambda$. Polar orbit, $\gamma_0$, in blue, equatorial orbit, $\gamma_{\pi/2}$, in red.](fig/lyapunovs.pdf)

Competition between noise and inertia
=====================================

The indeterminacy of the Jeffery orbits poses a problem for predictions involving the orientational distribution of particles. Jeffery [@Jeff22] hypothesised that orbits corresponding to minimal energy dissipation should be preferred. Another approach, pursued by a number of authors [@Sch55; @Hin72; @Bre74], shows that thermal noise, after a long time, removes the dependence on initial conditions and produces a stationary distribution of the ${\ensuremath{\mbox{\boldmath$n$}}}$ vector. In particular, in [@Hin72] it is shown that the orientational distribution converges to a function independent of noise strength, as the noise strength is reduced. On the contrary, in the preceding sections we have shown that in the absence of noise a weakly inertial particle will rotate into a specific stable orbit. In this section we consider what effect weak particle inertia have on the orientational distribution.

The time scale for inertia to have an effect is equal to the inverse Floquet exponent, $s /\gamma \sim 1/{\ensuremath{\widetilde {\ensuremath{\mbox{St}}}}}$. The corresponding time scale for diffusion is $s \tau_D \sim {\ensuremath{\mbox{Pe}}}$, where the Péclet number ${\ensuremath{\mbox{Pe}}}=s/D$ measures the noise strength, or diffusion constant $D=k_B T/\widetilde{Y}^C$, relative to the shear strength [@Kim91]. Thus, when $\tau_D\gamma \sim {\ensuremath{\mbox{Pe}}}{\ensuremath{\widetilde {\ensuremath{\mbox{St}}}}}\ll 1$ we expect noise to determine the stationary distribution, as described by previous authors [@Bre74]. On the other hand, when ${\ensuremath{\mbox{Pe}}}{\ensuremath{\widetilde {\ensuremath{\mbox{St}}}}}\gg1$ we expect a distribution dominated by the inertial effects, with a peak around the stable orbit described in the previous section. In order to quantify these expectations we need to compute moments of the stationary distribution for various ${\ensuremath{\mbox{Pe}}}$, ${\ensuremath{\widetilde {\ensuremath{\mbox{St}}}}}$ and $\lambda$.

The Fokker-Planck equation describing the evolution of the orientational probability density $P({\ensuremath{\mbox{\boldmath$n$}}}, t)$ is $$\begin{aligned}
    &\partial_t P({\ensuremath{\mbox{\boldmath$n$}}}, t) = -\partial_{{\ensuremath{\mbox{\boldmath$n$}}}}\left[(F_0({\ensuremath{\mbox{\boldmath$n$}}}) + {\ensuremath{\mbox{St}}}F_1({\ensuremath{\mbox{\boldmath$n$}}}))P\right] + {\ensuremath{\mbox{Pe}}}^{-1}\partial^2_{{\ensuremath{\mbox{\boldmath$n$}}}}P \equiv \hat J P. {\label{eqn:fpe}} \end{aligned}$$ Here $\partial_{{\ensuremath{\mbox{\boldmath$n$}}}}$ is the gradient on the sphere defined as $\partial_{{\ensuremath{\mbox{\boldmath$n$}}}} \equiv ({{\ensuremath{\mathbb{I}}}}- {{\ensuremath{\mbox{\boldmath$n$}}} {\ensuremath{\mbox{\boldmath$n$}}}{^{\mathrm T}}})\nabla$. The boundary condition of [Eq. (\[eqn:fpe\])]{} is that the probability density $P$ must integrate to unity over the sphere. The Fokker-Planck equation can be solved numerically by expansion in spherical harmonics, the details of this procedure are described in \[app:fpe\].

For a discussion of the solutions of [Eq. (\[eqn:fpe\])]{} we compute the stationary average $\langle \sin^2\theta\rangle$, where $\theta$ is the polar angle from the vorticity vector of the simple shear flow (see Fig. \[fig:coordinates\]). Recall that weakly inertial rod-shaped particles exhibit a stable orbit at $\theta = \pi/2$, and disk-shaped particles at $\theta=0$. Moreover, values of $\langle \sin^2\theta\rangle$ have been tabulated (Table 6c in Ref. [@Bre74]) for different values of ${\ensuremath{\mbox{Pe}}}$ and $\lambda$. This allows us to validate our method at ${\ensuremath{\mbox{St}}}=0$.

In Fig. \[fig:sintheta1\] we show the angular average $\langle \sin^2\theta\rangle$ as a function of ${\ensuremath{\mbox{Pe}}}$ for different values of ${\ensuremath{\mbox{St}}}$, and for two particle shapes: prolate $\lambda=5$ and oblate $\lambda=1/5$. We note the following features. At strong noise, that is small ${\ensuremath{\mbox{Pe}}}$, the moment is independent of ${\ensuremath{\mbox{St}}}$ and approaches the uniformly random result $\langle \sin^2\theta\rangle=2/3$. Next we see that the curve for ${\ensuremath{\mbox{St}}}=0$ (red curve) levels out and converges at a finite value of ${\ensuremath{\mbox{Pe}}}$, as predicted in Ref. [@Hin72]. Finally we see that for finite ${\ensuremath{\mbox{St}}}$, the curves coincide with the ${\ensuremath{\mbox{St}}}=0$ curve up to ${\ensuremath{\mbox{Pe}}}{\ensuremath{\mbox{St}}}\approx1$, but then deviates and instead converges towards $\langle\sin^2\theta\rangle=1$ (prolate particle), or $\langle\sin^2\theta\rangle=0$ (oblate particle) when ${\ensuremath{\mbox{Pe}}}{\ensuremath{\mbox{St}}}\gg 1$.

The limit of ${\ensuremath{\mbox{St}}}\to0$ is singular: the stationary distribution at low noise and weak inertia differs substantially from the one predicted by a ${\ensuremath{\mbox{St}}}=0$ calculation, although the time to arrive at the steady state becomes longer for weaker inertia. The most important observation is that the stationary average at low noise strength but without inertia is the same for prolate and oblate particles, whereas the inertial correction has opposite directions for the two different shapes. We emphasise this observation in Fig. \[fig:sintheta2\], where we show the stationary average as function of ${\ensuremath{\mbox{Pe}}}{\ensuremath{\widetilde {\ensuremath{\mbox{St}}}}}$, for different values of ${\ensuremath{\widetilde {\ensuremath{\mbox{St}}}}}/{\ensuremath{\mbox{Pe}}}$. It turns out that when ${\ensuremath{\mbox{Pe}}}{\ensuremath{\widetilde {\ensuremath{\mbox{St}}}}}$ is large enough, then the average does not depend on ${\ensuremath{\widetilde {\ensuremath{\mbox{St}}}}}$ and ${\ensuremath{\mbox{Pe}}}$ separately, but only in the combination ${\ensuremath{\mbox{Pe}}}{\ensuremath{\widetilde {\ensuremath{\mbox{St}}}}}$. The stationary average is determined by the ratio of relaxation times $\tau_D\gamma \sim{\ensuremath{\mbox{Pe}}}{\ensuremath{\widetilde {\ensuremath{\mbox{St}}}}}$, and by the aspect ratio $\lambda$. For small values of ${\ensuremath{\mbox{Pe}}}{\ensuremath{\widetilde {\ensuremath{\mbox{St}}}}}$ the stationary distribution is instead determined by ${\ensuremath{\mbox{Pe}}}$, as shown in Fig. \[fig:sintheta1\].

![\[fig:sintheta1\] Stationary orientational average as function of ${\ensuremath{\mbox{Pe}}}$ for different ${\ensuremath{\widetilde {\ensuremath{\mbox{St}}}}}$ and $\lambda$. Solid lines are a prolate particle $\lambda=5$, dashed lines oblate $\lambda=1/5$. The red solid line is ${\ensuremath{\widetilde {\ensuremath{\mbox{St}}}}}=0$ result, then blue ${\ensuremath{\widetilde {\ensuremath{\mbox{St}}}}}= 10^{-4}$, green ${\ensuremath{\widetilde {\ensuremath{\mbox{St}}}}}= 10^{-3}$, purple ${\ensuremath{\widetilde {\ensuremath{\mbox{St}}}}}= 10^{-2}$ and orange ${\ensuremath{\widetilde {\ensuremath{\mbox{St}}}}}= 10^{-1}$. ](fig/sintheta1.pdf)

![\[fig:sintheta2\] Stationary orientational average as function of ${\ensuremath{\mbox{Pe}}}{\ensuremath{\widetilde {\ensuremath{\mbox{St}}}}}$ for different values of $k={\ensuremath{\widetilde {\ensuremath{\mbox{St}}}}}/{\ensuremath{\mbox{Pe}}}$. Solid lines are a prolate particle $\lambda=5$, dashed lines oblate $\lambda=1/5$. Then red $k = 10^{-9}$, blue $k = 10^{-8}$, green $k = 10^{-7}$, purple $k = 10^{-6}$ and orange $k = 10^{-5}$.](fig/sintheta2.pdf)

Conclusions {#sec:conc}
===========

The orientational dynamics of weakly inertial axisymmetric particles in a steady flow has been investigated. We have derived an approximate equation of motion for the unit axial vector along the particle symmetry axis, assuming that the Stokes number is small. The fact that we have obtained this equation of motion in explicit form (valid for any axisymmetric and weakly inertial particle in any steady linear viscous flow) has made it possible to address two important questions, both pertaining to the case of a simple shear flow.

First, we have shown that axisymmetric particles drift towards either the log-rolling (disk-shaped particles) or the tumbling orbit (rod-like particles). We have quantified this result for any particle aspect ratio by calculating the Floquet exponents of the periodic orbits to leading order in ${\ensuremath{\mbox{St}}}$. Our result generalises a result for nearly spherical particles in Ref. [@Sub06] to arbitrary axisymmetric particles.

Second, the approximate equation of motion derived in Sec. 2 has allowed us to formulate the Fokker-Planck equation describing the orientational distribution of axisymmetric, weakly inertial particles in a shear flow when rotational noise is present. Numerical solutions of the Fokker-Planck equation show that stationary orientational distributions are close to the inertia-free case when ${\ensuremath{\mbox{Pe}}}{\ensuremath{\widetilde {\ensuremath{\mbox{St}}}}}\ll1$, whereas they are determined by inertial effects, though small, when ${\ensuremath{\mbox{Pe}}}\gg1/{\ensuremath{\widetilde {\ensuremath{\mbox{St}}}}}\gg1$.

An important open question is to describe the effect of fluid inertia (at small but finite particle Reynolds number) upon the orientational motion for particles with arbitrary aspect ratios, in competition with the effects of particle inertia and rotational noise.

The calculations summarised in this article are restricted to axisymmetric particles. But our approach is readily generalised to asymmetric (triaxial) particles.

Acknowledgements {#acknowledgements .unnumbered}
================

Financial support by Vetenskapsrådet, by the Göran Gustafsson Foundation for Research in Natural Sciences and Medicine, and by the EU COST Action MP0806 on Particles in Turbulence’ are gratefully acknowledged.

[10]{} url \#1[`#1`]{}urlprefixhref \#1\#2[\#2]{} \#1[\#1]{}

C. Rolf, M. Kramer, C. Schiller, M. Hildebrandt, M. Riese, Lidar observation and model simulation of a volcanic-ash-induced cirrus cloud during the Eyjafjallajökull eruption, Atmos. Chem. Phys. Discuss. 12 (2012) 15675–15707.

H. Stommel, Trajectories of small bodies sinking slowly through convection cells., J. Marine Res. 8 (1949) 24–29.

J. S. Guasto, R. Rusconi, R. Stocker, Fluid mechanics of planktonic microorganisms, Annual Review of Fluid Mechanics 44 (1) (2012) 373–400.

F. Lundell, L. D. S[ö]{}derberg, P. H. Alfredsson, Fluid mechanics of papermaking, Annual Review of Fluid Mechanics 43 (1) (2011) 195–217.

G. B. Jeffery, The motion of ellipsoidal particles immersed in a viscous fluid, Proceedings of the Royal Society of London Series A 102 (715) (1922) 161–179.

H. Brenner, Rheology of a dilute suspension of axisymmetric Brownian particles, International Journal of Multiphase Flow 1 (2) (1974) 195–341.

C. J. Petrie, The rheology of fibre suspensions, Journal of Non-Newtonian Fluid Mechanics 87 (2–3) (1999) 369–402.

G. Gauthier, P. Gondret, M. Rabaud, Motions of anisotropic particles: Application to visualization of three-dimensional flows, Phys. Fluids 10 (9) (1998) 2147–2154.

M. Wilkinson, V. Bezuglyy, B. Mehlig, Fingerprints of random flows?, Physics of Fluids 21 (4) (2009) 043304.

M. Wilkinson, V. Bezuglyy, B. Mehlig, Emergent order in rheoscopic swirls, Journal of Fluid Mechanics 667 (2011) 158–187.

V. Bezuglyy, B. Mehlig, M. Wilkinson, Poincar[é]{} indices of rheoscopic visualisations, Europhysics Letters 89 (2010) 34003.

J. Einarsson, A. Johansson, S. K. Mahato, Y. Mishra, J. R. Angilella, D. Hanstorp, B. Mehlig, Periodic and aperiodic tumbling of microrods advected in a microchannel flow, [Acta Mechanica]{}, in press (2013).

A. Pumir, M. Wilkinson, Orientation statistics of small particles in turbulence, New Journal of Physics 13 (9) (2011) 093030.

S. Parsa, E. Calzavarini, F. Toschi, G. A. Voth, Rotation rate of rods in turbulent fluid flow, Physical Review Letters 109 (2012) 134501.

F. Lundell, A. Carlsson, Heavy ellipsoids in creeping shear flow: Transitions of the particle rotation rate and orbit shape, Phys. Rev. E 81 (2010) 016323.

F. Lundell, The effect of particle inertia on triaxial ellipsoids in creeping shear: From drift toward chaos to a single periodic solution, Phys. Fluids 23 (2011) 011704.

G. Subramanian, D. L. Koch, Inertial effects on the orientation of nearly spherical particles in simple shear flow, Journal of Fluid Mechanics 557 (2006) 257–296.

P. M. Lovalenti, J. F. Brady, The hydrodynamic force on a rigid particle undergoing arbitrary time-dependent motion at small reynolds number, Journal of Fluid Mechanics 256 (1993) 561–605.

C. Marchioli, M. Fantoni, A. Soldati, Orientation, distribution, and deposition of elongated, inertial fibers in turbulent channel flow, Physics of Fluids 22 (3) (2010) 033301.

S. Kim, S. Karrila, Microhydrodynamics: Principles And Selected Applications, Butterworth - Heinemann series in Chemical Engineering, 1991.

K. Gustavsson, J. Einarsson, B. Mehlig, Tumbling of small axisymmetric particles in random and turbulent flows, [arXiv]{} e-print 1305.1822 (2013).

J. P. Snyder, Map Projections: A Working Manual, Geological Survey (U.S.), 1987. [http://pubs.er.usgs.gov/publication/pp1395]{}.

H. A. Scheraga, [Non–Newtonian]{} viscosity of solutions of ellipsoidal particles, The Journal of Chemical Physics 23 (8) (1955) 1526–1532.

E. J. Hinch, L. G. Leal, The effect of Brownian motion on the rheological properties of a suspension of non-spherical particles, Journal of Fluid Mechanics 52 (04) (1972) 683–712.

G. Arfken, Mathematical Methods for Physicists, Second Edition, Academic Press, 1970.

J. J. Sakurai, S. F. Tuan (ed.), Modern quantum mechanics, Addison-Wesley Pub. 1994.

<https://github.com/jeinarsson/fpe_sphere>.

Solution of the Fokker-Planck equation {#app:fpe}
======================================

We start from the Fokker-Planck equation governing the evolution of the probability density $P({\ensuremath{\mbox{\boldmath$n$}}}, t)$ of finding a particle with orientation vector ${\ensuremath{\mbox{\boldmath$n$}}}$ at time $t$: $$\begin{aligned}
    \partial_t P({\ensuremath{\mbox{\boldmath$n$}}}, t) &= -\partial_{{\ensuremath{\mbox{\boldmath$n$}}}}\left[\dot{{\ensuremath{\mbox{\boldmath$n$}}}}P\right] + {\ensuremath{\mbox{Pe}}}^{-1}\partial^2_{{\ensuremath{\mbox{\boldmath$n$}}}}P \equiv \hat J P {\label{eqn:fpeapp}} \\
    \intertext{with the normalisation condition}
    \int_{S_2} \!\!\!{\ensuremath{\mathrm{d}}}{\ensuremath{\mbox{\boldmath$n$}}}\,P({\ensuremath{\mbox{\boldmath$n$}}}, t) &= 1.\end{aligned}$$ The differential operator $\partial_{{\ensuremath{\mbox{\boldmath$n$}}}}$ is the gradient on the sphere, defined by taking the usual gradient $\nabla$ in $\mathbb R^3$ projected onto the unit sphere, $\partial_{{\ensuremath{\mbox{\boldmath$n$}}}} \equiv ({{\ensuremath{\mathbb{I}}}}- {{\ensuremath{\mbox{\boldmath$n$}}} {\ensuremath{\mbox{\boldmath$n$}}}{^{\mathrm T}}})\nabla$. The drift term $\dot{{\ensuremath{\mbox{\boldmath$n$}}}}$ may be any polynomial function of ${\ensuremath{\mbox{\boldmath$n$}}}$, for example the equation of motion [Eq. (\[eqn:nequation\])]{}.

We approximate the solution of [Eq. (\[eqn:fpeapp\])]{} by an expansion in spherical harmonics [@Sch55], the eigenfunctions of the quantum-mechanical angular-momentum operators, which form a complete basis on $S_2$. We use bra-ket notation, $$\begin{aligned}
    P({\ensuremath{\mbox{\boldmath$n$}}}, t) &= \langle {\ensuremath{\mbox{\boldmath$n$}}}| P(t)\rangle =  \sum_{l=0}^{\infty}\sum_{m=-l}^l c_l^m(t) \langle {\ensuremath{\mbox{\boldmath$n$}}}{|l,m\rangle}{\label{eqn:basisexpansion}}
\intertext{where}
     \langle{\ensuremath{\mbox{\boldmath$n$}}}|l,m\rangle &= Y_l^m({\ensuremath{\mbox{\boldmath$n$}}}) = (-1)^m\sqrt{\frac{2l+1}{4\pi}\frac{(l-m)!}{(l+m)!}}P_l^m(\cos\theta)e^{im\varphi}.\end{aligned}$$ We use the standard spherical harmonics defined in for example [@Arf70] (p. 571). The functions $P_l^m$ are the associated Legendre polynomials. We call the time-dependent coefficients for each basis function $c_l^m(t)$, and the fact that $P$ is real-valued constrains the coefficients: $$\begin{aligned}
    c_l^{-m} = (-1)^m\overline{c_l^m},{\label{eqn:clmreal}}\end{aligned}$$ where the bar denotes complex conjugation. Then [Eq. (\[eqn:fpeapp\])]{} for the time evolution of the state $|P(t)\rangle$ reads $$\begin{aligned}
    \frac{{\ensuremath{\mathrm{d}}}}{{\ensuremath{\mathrm{d}}}t} |P(t)\rangle &= \hat J |P(t)\rangle\,.\end{aligned}$$ Inserting the expansion yields $$\begin{aligned}
    \sum_{l=0}^{\infty}\sum_{m=-l}^l \frac{{\ensuremath{\mathrm{d}}}}{{\ensuremath{\mathrm{d}}}t} c_l^m(t) {|l,m\rangle} &= \sum_{l=0}^{\infty}\sum_{m=-l}^l c_l^m(t) \hat J{|l,m\rangle}.
    \intertext{Multiplying with $\langle p,q|$, and using the orthogonality of the spherical harmonics we arrive at a system of coupled ordinary differential equations for the coefficients $c_l^m(t)$}
    \dot c_p^q(t)  &= \sum_{l=0}^{\infty}\sum_{m=-l}^l c_l^m(t) \langle p,q| \hat J{|l,m\rangle} {\label{eqn:clmeq1}},\\
    c_0^0 &= \frac{1}{\sqrt{4\pi}}{\label{eqn:clmeq1norm}}.\end{aligned}$$ In order to solve [Eq. (\[eqn:clmeq1\])]{} it remains to compute the matrix elements of the operator $\hat J$. This is achieved by expressing $\hat J$ as a combination of the angular momentum operators. How the angular momentum operators act on the spherical harmonics is well known, see for example Ref. [@Arf70].

The angular momentum operator $\hat{{\ensuremath{\mbox{\boldmath$L$}}}}$ is given by $$\begin{aligned}
    \hat{{\ensuremath{\mbox{\boldmath$L$}}}} &= -i \hat{{\ensuremath{\mbox{\boldmath$n$}}}}{\wedge}\nabla.
\intertext{In terms of $\hat{{\ensuremath{\mbox{\boldmath$L$}}}}$ we have}
    \partial_{{\ensuremath{\mbox{\boldmath$n$}}}}^2 &= ({{\ensuremath{\mathbb{I}}}}-\hat{{\ensuremath{\mbox{\boldmath$n$}}}}\hat{{\ensuremath{\mbox{\boldmath$n$}}}}{^{\mathrm T}})\nabla \cdot ({{\ensuremath{\mathbb{I}}}}-\hat{{\ensuremath{\mbox{\boldmath$n$}}}}\hat{{\ensuremath{\mbox{\boldmath$n$}}}}{^{\mathrm T}})\nabla = - \hat{{\ensuremath{\mbox{\boldmath$L$}}}}^2,
\intertext{and}
    \partial_{{\ensuremath{\mbox{\boldmath$n$}}}}&= ({{\ensuremath{\mathbb{I}}}}-\hat{{\ensuremath{\mbox{\boldmath$n$}}}}\hat{{\ensuremath{\mbox{\boldmath$n$}}}}{^{\mathrm T}})\nabla =  -i\hat{{\ensuremath{\mbox{\boldmath$n$}}}} {\wedge}\hat{{\ensuremath{\mbox{\boldmath$L$}}}}.\end{aligned}$$ The states ${|l,m\rangle}$ are the eigenfunctions of $\hat{{\ensuremath{\mbox{\boldmath$L$}}}}^2$ and $\hat{L}_3$ with eigenvalues $l(l+1)$ and $m$. Now, to evaluate the drift term $\partial_{{\ensuremath{\mbox{\boldmath$n$}}}}\dot{{\ensuremath{\mbox{\boldmath$n$}}}}$ we need to know how $\hat{L}_1$, $\hat{L}_2$, and $\hat{{\ensuremath{\mbox{\boldmath$n$}}}}$ act on ${|l,m\rangle}$. The first two are known through the use of ladder operators, defined by $$\begin{aligned}
    \hat{L}_{\pm} = \hat{L}_1 \pm i\hat{L}_2 \implies \hat{L}_1 = \frac{1}{2}(\hat{L}_+ + \hat{L}_-),\quad \hat{L}_2 = -\frac{i}{2}(\hat{L}_+ - \hat{L}_-)\end{aligned}$$ and their effect is $$\begin{aligned}
    \hat{L}_\pm {|l,m\rangle} &= \sqrt{(l\mp m)(1+l\pm m)}\,{|l,m\pm1\rangle}.\end{aligned}$$ Next, the drift term we consider is a polynomial in $\hat{n}_1, \hat{n}_2$ and $\hat{n}_3$ (the components of $\hat{{\ensuremath{\mbox{\boldmath$n$}}}}$), we therefore need to evaluate the action of a monomial $\hat{n}_1^\alpha \hat{n}_2^\beta \hat{n}_3^\gamma$ on ${|l,m\rangle}$. Any such monomial of order $k=\alpha+\beta+\gamma$ may be written as a linear combination of spherical tensor operators of up to order $k$: $$\begin{aligned}
    \hat{n}_1^\alpha \hat{n}_2^\beta \hat{n}_3^\gamma &= \sum_{l=0}^{k}\sum_{m=-l}^l a_l^m(\alpha,\beta,\gamma)\hat{Y}_l^m.\end{aligned}$$ The action of $\hat{Y}_p^q$ is computed with the Clebsch-Gordan coefficients (see for example Ref. [@Sak94], p.216) $$\begin{aligned}
    \hat{Y}_p^q {|l,m\rangle} &= \sum_{\Delta l=-p}^{p}K(p,q,l,m,l+\Delta l, m+q){|l+\Delta l,m+q\rangle},\\\end{aligned}$$ where $$\begin{aligned}
    K(l_1, m_1, l_2, m_2, l, m) &= \sqrt{\frac{(2l_1+1)(2l_2+1)}{4\pi (2l+1)}} \\
    &\quad \times \langle l_1 l_2; 0 0|l_1l_2;l0\rangle
    \langle l_1 l_2; m_1 m_2|l_1l_2;lm\rangle\,.\end{aligned}$$ The Clebsch-Gordan coefficients are denoted by $\langle l_1 l_2; m_1 m_2|l_1l_2;lm\rangle$ (see e.g Eq. (3.7.44) of Ref. [@Sak94]), and they are available in Mathematica by the function `ClebschGordan[l1,m1,l2,m2,l,m]`. Finally, we order the operators, so that as many terms as possible cancel before we actually begin evaluating the operators. In particular we want to reorder $\hat{n}_i$ against $\hat{L}_i$, and we make use of the commutator $$\begin{aligned}
    [\hat{n}_p, \hat{L}_q] &=\sum_{j=1}^3 i\varepsilon_{pqj}\hat{n}_j,\end{aligned}$$ where $\varepsilon_{ijk}$ is the Levi-Civita tensor.

Let us consider an example. In the special case of $\hat J = {\ensuremath{\mbox{Pe}}}^{-1} \partial_{{\ensuremath{\mbox{\boldmath$n$}}}}^2$ we see that $\hat J=-{\ensuremath{\mbox{Pe}}}^{-1} \hat{{\ensuremath{\mbox{\boldmath$L$}}}}^2$. It follows: $$\begin{aligned}
    \hat J {|l,m\rangle}&=-{\ensuremath{\mbox{Pe}}}^{-1}\hat{{\ensuremath{\mbox{\boldmath$L$}}}}^2 {|l,m\rangle} = -l(l+1){\ensuremath{\mbox{Pe}}}^{-1}{|l,m\rangle}.\end{aligned}$$ Knowing how $\hat J$ acts on ${|l,m\rangle}$ fully specifies [Eq. (\[eqn:clmeq1\])]{} which becomes $$\begin{aligned}
    \dot c_p^q(t)  &= -p(p+1){\ensuremath{\mbox{Pe}}}^{-1}c_p^q(t).\end{aligned}$$ This means that the diffusion operator exponentially suppresses all modes with $p>0$, and the lowest mode $p=0$ is determined by the normalisation condition [Eq. (\[eqn:clmeq1norm\])]{}. This is the solution of the diffusion equation on the sphere, with the uniform distribution as steady state.

Now take the for the drift term the standard Jeffery equation $\dot {{\ensuremath{\mbox{\boldmath$n$}}}} = {\ensuremath{\mathbb{O}}} {\ensuremath{\mbox{\boldmath$n$}}} + \Lambda ({\ensuremath{\mathbb{S}}} {\ensuremath{\mbox{\boldmath$n$}}} - {{\ensuremath{\mbox{\boldmath$n$}}} {\ensuremath{\mbox{\boldmath$n$}}}{^{\mathrm T}}}{\ensuremath{\mathbb{S}}} {\ensuremath{\mbox{\boldmath$n$}}}) = {\ensuremath{\mbox{\boldmath$F$}}}_0({\ensuremath{\mbox{\boldmath$n$}}})$, and call this operator $\hat J = -\partial_{{\ensuremath{\mbox{\boldmath$n$}}}}{\ensuremath{\mbox{\boldmath$F$}}}_0({\ensuremath{\mbox{\boldmath$n$}}}) + {\ensuremath{\mbox{Pe}}}^{-1}\partial^2_{{\ensuremath{\mbox{\boldmath$n$}}}}$. We find $$\begin{aligned}
    \hat J &= 
    i \sqrt{\frac{\pi }{30}} \Lambda  \hat{L}_-\hat{Y}_{2}^{-1}+i \sqrt{\frac{\pi }{30}} \hat{L}_-\hat{Y}_{2}^{1}-i \sqrt{\frac{\pi }{30}} 
   \Lambda  \hat{L}_+\hat{Y}_{2}^{1} \\ 
   &\quad-  i \sqrt{\frac{\pi }{30}} \hat{L}_+\hat{Y}_{2}^{-1} 
   -i \sqrt{\frac{2 \pi }{15}} \Lambda 
   \hat{L}_3\hat{Y}_{2}^{-2}-i \sqrt{\frac{2 \pi }{15}} \Lambda  \hat{L}_3\hat{Y}_{2}^{2} \\ 
   &\quad
   +\frac{2}{3} i \sqrt{\pi } 
   \hat{L}_3\hat{Y}_{0}^{0}-\frac{2}{3} i \sqrt{\frac{\pi }{5}} \hat{L}_3\hat{Y}_{2}^{0}-{\ensuremath{\mbox{Pe}}}^{-1} \hat{{\ensuremath{\mbox{\boldmath$L$}}}}^2.\end{aligned}$$ The fact that the operator $\hat J$ is real implies that the matrix elements must have the symmetry $$\begin{aligned}
    \langle p, q|\hat J |l,m\rangle = (-1)^q \overline{\langle p, -q|\hat J |l,m\rangle}.\end{aligned}$$ We can understand the reason for this because the system of differential equations in [Eq. (\[eqn:clmeq1\])]{} needs to preserve the condition that $P$ is real, as stated in [Eq. (\[eqn:clmreal\])]{}. It implies we have only to compute half of the matrix elements.

Upon including the correction [Eq. (\[eqn:nequation\])]{} due to weak particle inertia, the expression for the operator $\hat{J}$ becomes too lengthy to include here. However, we have implemented the identities and rules described in this appendix in a Mathematica notebook [@githubrepo]. Given an operator $\hat J$ expressed in terms of $\hat{{\ensuremath{\mbox{\boldmath$n$}}}}$ and $\hat{{\ensuremath{\mbox{\boldmath$L$}}}}$, it converts the expression into sums of angular momentum operators as exemplified above. We also provide an additional Mathematica notebook that, given the matrix elements, assembles a sparse matrix ${\ensuremath{\mathbb{J}}}$, up to a desired order $l_{\textrm{max}}$. This matrix can then be used to solve the truncated version of [Eq. (\[eqn:clmeq1\])]{}: $$\begin{aligned}
    \dot {{\ensuremath{\mbox{\boldmath$c$}}}} &= {\ensuremath{\mathbb{J}}} {\ensuremath{\mbox{\boldmath$c$}}},\end{aligned}$$ where ${\ensuremath{\mbox{\boldmath$c$}}}$ is a vector with elements $c_l^m$. Alternatively, the matrix can be solved for the stationary solution of ${\ensuremath{\mbox{\boldmath$c$}}}$. All results shown in this paper were computed using $l_{\mathrm{max}}=400$, which leads to very good convergence in all cases shown. Only when solutions approach delta peaks, as for example the limiting stable orbits in the large ${\ensuremath{\mbox{Pe}}}{\ensuremath{\mbox{St}}}$ case, the expansion procedure converges very slowly.
