2.1. The generalized reciprocal theorem

We use the generalized reciprocal theorem to obtain an expression for the angular velocity of a neutrally buoyant spheroid in a viscoelastic simple shear flow. The theorem relates two sets of velocity and stress fields – $(\boldsymbol{\sigma }^{(1)}, \boldsymbol{u}^{(1)} )$ pertains to the problem of interest above, a freely rotating spheroid of an arbitrary aspect ratio in simple shear flow of a second-order fluid, while $(\boldsymbol{\sigma }^{(2)}, \boldsymbol{u}^{(2)} )$ corresponds to the test problem, taken to be a spheroid rotating about an axis transverse to its symmetry axis, in an otherwise quiescent Newtonian fluid, at zero Reynolds number. In terms of the disturbance fields (indicated by primes) of the actual problem, the theorem yields the identity

(2.1) \begin{equation} \int _{S_{p}} \big (\boldsymbol{\sigma }^{(2)} \boldsymbol{\cdot }\boldsymbol{u}^{\prime (1)}-\boldsymbol{\sigma }^{\prime (1)} \boldsymbol{\cdot }\boldsymbol{u}^{(2)} \big)\boldsymbol{\cdot }\boldsymbol{n} \,{\rm d}S=De \int _{V} \boldsymbol{\sigma }^{\prime (1)}_{\textit{NN}} : \boldsymbol{\nabla }\boldsymbol{u}^{(2)} \,{\rm d}V, \end{equation}

where $\boldsymbol{n}$ is the unit outer normal to the spheroid surface ( $S_p$ ). The disturbance fields have been assumed to decay away sufficiently rapidly at large distances, so as to ensure the absence of any contributions from the surface integral at infinity. The latter is typically true in the small- ${\textit{De}}$ limit, owing to the ordered fluid constitutive equation being obtained from a truncation of the retarded motion expansion (see below), and the latter being essentially a regular perturbation expansion about the Newtonian limit (Bird et al. Reference Bird, Armstrong and Hassager1987).

The non-dimensional equations of motion and the continuity equation, for the problem of interest, are given by

(2.2) \begin{equation} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma }^{(1)}=0, \quad \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}^{(1)}=0, \end{equation}

with

(2.3) \begin{equation} \left . \begin{array}{ll} \boldsymbol{u}^{(1)}=\boldsymbol{\varOmega }^{(1)} \wedge \boldsymbol{x} \quad \text{for} \quad \boldsymbol{x} \in S_{P},\\ \boldsymbol{u}^{(1)} \rightarrow \boldsymbol{\varGamma } \boldsymbol{\cdot }\boldsymbol{x} \qquad \, \text{for} \quad \boldsymbol{x} \rightarrow \infty . \end{array}\!\! \right \} \end{equation}

The stress field in (2.2) is defined by

(2.4) \begin{equation} \boldsymbol{\sigma }^{(1)}=\boldsymbol{\sigma }_{N}^{(1)}+{\textit{De}}\, \boldsymbol{\sigma }_{\textit{NN}}^{(1)}, \end{equation}

where $\boldsymbol{\sigma }_{N}^{(1)}=-p^{(1)}\boldsymbol{I}+2 \boldsymbol{e}^{(1)}$ and $\boldsymbol{\sigma }_{\textit{NN}}^{(1)}$ are the Newtonian and non-Newtonian stress contributions, respectively. For ${\textit{De}} \ll 1$ , the retarded motion expansion can be employed to incorporate weak elastic effects in a perturbative manner. The first correction to Newtonian rheology corresponds then to a second-order fluid, and at this order, $\boldsymbol{\sigma }_{\textit{NN}}^{(1)}$ in (2.4) may be written as the sum of corotational ( $\boldsymbol{\sigma }_{\textit{NNC}}^{(1)}$ ) and quadratic ( $\boldsymbol{\sigma }_{\textit{NNQ}}^{(1)}$ ) contributions. Thus, $\boldsymbol{\sigma }^{(1)}_{\textit{NN}}= \boldsymbol{\sigma }^{(1)}_{\textit{NNC}} +\boldsymbol{\sigma }^{(1)}_{\textit{NNQ}}$ , with

(2.5) \begin{align} \boldsymbol{\sigma }^{(1)}_{\textit{NNC}}&=2 \epsilon \left ( \frac {\partial \boldsymbol{e}^{(1)}}{\partial t}+\boldsymbol{u}^{(1)} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{e}^{(1)} +\boldsymbol{w}^{(1)} \boldsymbol{\cdot }\boldsymbol{e}^{(1)}+\big(\boldsymbol{w}^{(1)} \boldsymbol{\cdot }\boldsymbol{e}^{(1)}\big)^{\dagger } \right )\!, \end{align}

(2.6) \begin{align} \boldsymbol{\sigma }^{(1)}_{\textit{NNQ}}&=4 (1+\epsilon ) \big ( \boldsymbol{e}^{(1)} \boldsymbol{\cdot }\boldsymbol{e}^{(1)} \big ), \end{align}

where $\boldsymbol{e}^{(1)}=(\boldsymbol{\nabla }\boldsymbol{u}^{(1)}+\boldsymbol{\nabla }\boldsymbol{u}^{(1)^{\dagger }})/2$ and $\boldsymbol{w}^{(1)}=(\boldsymbol{\nabla }\boldsymbol{u}^{(1)}-\boldsymbol{\nabla }\boldsymbol{u}^{(1)^{\dagger }})/2$ are the rate of strain and vorticity tensors, respectively, with $\epsilon$ being the second-order fluid parameter that, as mentioned in § 1, is related to the ratio of the normal stress differences (Subramanian & Koch Reference Subramanian and Koch2006, Reference Subramanian and Koch2007). The disturbance velocity field in (2.1) is defined by $\boldsymbol{u}^{\prime (1)}= \boldsymbol{u}^{(1)}-\boldsymbol{\varGamma } \boldsymbol{\cdot }\boldsymbol{x}$ , where $\boldsymbol{\varGamma } \boldsymbol{\cdot }\boldsymbol{x}$ is the ambient simple shear flow, with $\boldsymbol{\varGamma }$ being the transpose of the ambient velocity gradient tensor. On $S_{p}$ , $\boldsymbol{u}^{\prime (1)}= \boldsymbol{\varOmega }^{(1)} \wedge \boldsymbol{x}-\boldsymbol{\varGamma } \boldsymbol{\cdot }\boldsymbol{x}$ , with $\boldsymbol{\varOmega }^{(1)}$ being the spheroid angular velocity determined from a torque-free constraint ( ${\boldsymbol{\mathcal L}}^{(1)} = 0$ ).

In the above, $\boldsymbol{\varGamma }=\boldsymbol{1}^{\prime }_{x} \boldsymbol{1}^{\prime }_{y}$ in a space-fixed coordinate system with axes oriented along the flow ( $X^{\prime }$ ), gradient ( $Y^{\prime }$ ) and vorticity ( $Z^{\prime }$ ) directions of the ambient simple shear; $\boldsymbol{E}=(\boldsymbol{\varGamma } +\boldsymbol{\varGamma }^{\dagger })/2$ and $\boldsymbol{W}=(\boldsymbol{\varGamma }^{\dagger }-\boldsymbol{\varGamma })/2$ denote the (non-dimensional) ambient rate of strain and vorticity tensors, respectively. One may now write the governing equations for the actual problem as

(2.7) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma }^{ \prime (1)}=-\boldsymbol{\nabla }{\!p}^{\prime (1)} + {\nabla} ^{2} \boldsymbol{u}^{ \prime (1)} +{\textit{De}} \,\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma }_{\textit{NN}}^{ \prime (1)}&=0, \end{align}

(2.8) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}^{ \prime (1)}&=0, \end{align}

in terms of the disturbance fields. The disturbance stress field in (2.1) and (2.7) is defined as $\boldsymbol{\sigma }^{\prime (1)}= \boldsymbol{\sigma }^{ (1)}-\boldsymbol{\sigma }^{ \infty }$ , where $\boldsymbol{\sigma }^{\infty }=\boldsymbol{\sigma }_{N}^{ \infty } +{\textit{De}} \boldsymbol{\sigma }_{\textit{NN}}^{ \infty }$ ; $\boldsymbol{\sigma }_{N}^{ \infty }=-p^{\infty } \boldsymbol{I}+2 \boldsymbol{E}$ and $\boldsymbol{\sigma }_{\textit{NN}}^{ \infty }= \boldsymbol{\sigma }_{\textit{NNC}}^{ \infty } +\boldsymbol{\sigma }_{\textit{NNQ}}^{ \infty }$ , with $\boldsymbol{\sigma }_{\textit{NNC}}^{ \infty }= 2 \epsilon ( \boldsymbol{W} \boldsymbol{\cdot }\boldsymbol{E}+(\boldsymbol{W} \boldsymbol{\cdot }\boldsymbol{E})^{\dagger } )$ and $\boldsymbol{\sigma }_{\textit{NNQ}}^{ \infty }= 4 (1+\epsilon ) ( \boldsymbol{E} \boldsymbol{\cdot }\boldsymbol{E} )$ , are the Newtonian and non-Newtonian ambient stress components. The corotational and quadratic components of $\boldsymbol{\sigma }_{\textit{NN}}^{ \prime (1)}$ are therefore given by

(2.9) \begin{align} \boldsymbol{\sigma }^{\prime (1)}_{\textit{NNC}}&=\boldsymbol{\sigma }^{(1)}_{\textit{NNC}}-\boldsymbol{\sigma }^{\infty }_{\textit{NNC}}\nonumber \\ &= 2 \epsilon \Bigg (\frac {\partial \boldsymbol{e}^{\prime (1)}}{\partial t}+\boldsymbol{u}^{(1)} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{e}^{\prime (1)}+\boldsymbol{w}^{\prime (1)} \boldsymbol{\cdot }\boldsymbol{e}^{\prime (1)}+\big(\boldsymbol{w}^{\prime (1)} \boldsymbol{\cdot }\boldsymbol{e}^{\prime (1)}\big)^{\dagger }+\boldsymbol{W} \boldsymbol{\cdot }\boldsymbol{e}^{\prime (1)} \nonumber \\ &\quad +\big(\boldsymbol{W} \boldsymbol{\cdot }\boldsymbol{e}^{\prime (1)}\big)^{\dagger } +\boldsymbol{w}^{\prime (1)} \boldsymbol{\cdot }\boldsymbol{E}+\big(\boldsymbol{w}^{\prime (1)} \boldsymbol{\cdot }\boldsymbol{E}\big)^{\dagger } \Bigg), \end{align}

(2.10) \begin{align} \boldsymbol{\sigma }_{\textit{NNQ}}^{\prime (1)}&= \boldsymbol{\sigma }^{(1)}_{\textit{NNQ}}-\boldsymbol{\sigma }^{\infty }_{\textit{NNQ}}\nonumber \\ &=4 (1+\epsilon ) \big(\boldsymbol{e}^{\prime (1)} \boldsymbol{\cdot } \boldsymbol{e}^{\prime (1)}+\boldsymbol{e}^{\prime (1)} \boldsymbol{\cdot }\boldsymbol{E}+\boldsymbol{E} \boldsymbol{\cdot e}^{\prime (1)}\big). \end{align}

Note that $\epsilon =-1$ , with $\boldsymbol{\sigma }^{\prime (1)}_{\textit{NNQ}}=0$ , corresponds to the corotational limit, when viscoelasticity leaves the Stokesian velocity field unchanged (Bird et al. Reference Bird, Armstrong and Hassager1987).

The Stokesian test problem is defined by

(2.11) \begin{equation} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma }^{(2)}=0, \quad \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}^{(2)}=0, \end{equation}

with

(2.12) \begin{equation} \left . \begin{array}{ll} \boldsymbol{u}^{(2)}=\boldsymbol{\varOmega }^{(2)} \wedge \boldsymbol{x} \quad \text{for} \quad \boldsymbol{x} \in S_{P},\\ \boldsymbol{u}^{(2)} \rightarrow 0 \qquad \quad \ \ \ \text{for} \quad \boldsymbol{x} \rightarrow \infty , \end{array}\!\! \right \} \end{equation}

where $\boldsymbol{\varOmega }^{(2)}$ is the spheroid angular velocity constrained to be orthogonal to its symmetry axis ( $\boldsymbol p$ ). On using the boundary conditions for the actual and test problems, (2.1) takes the form

(2.13) \begin{equation} \begin{aligned} \boldsymbol{\varOmega }^{(1)} \boldsymbol{\cdot }\mathcal{\boldsymbol{L}}^{(2)} = \boldsymbol{\varGamma } \boldsymbol{:} \int _{Sp} \boldsymbol{x} \big(\boldsymbol{\sigma }^{(2)} \boldsymbol{\cdot }\boldsymbol{n} \big) \,{\rm d}S +{\textit{De}} \, \int _{V} \big(\boldsymbol{\sigma }_{\textit{NNC}}^{\prime (1)}+\boldsymbol{\sigma }_{\textit{NNQ}}^{\prime (1)}\big) : \boldsymbol{\nabla }\boldsymbol{u}^{(2)} {\rm d}V, \end{aligned} \end{equation}

where $\mathcal{\boldsymbol{L}}^{(2)}$ is the torque acting on the test spheroid. Since there exists a Giesekus pressure field $p_{G}^{\prime }$ in the corotational limit such that $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma }^{\prime (1)}_{\textit{NNC}} -\boldsymbol{\nabla }{\!p}_{G}^{\prime (1)}=0$ , with

(2.14) \begin{equation} p_{G}^{\prime (1)}=\epsilon \left ( \frac {\partial p^{\prime (1)}}{\partial t} +\boldsymbol{u}^{(1)} \boldsymbol{\cdot } \boldsymbol{\nabla }{\!p}^{\prime (1)} + \boldsymbol{e}^{\prime (1)} : \boldsymbol{e}^{\prime (1)}+ 2 \big(\boldsymbol{e}^{\prime (1)} : \boldsymbol{E}\big)\right )\!, \end{equation}

the corotational term in (2.13) may be rewritten in terms of a surface integral as

(2.15) \begin{equation} \int _{V} \boldsymbol{\sigma }_{\textit{NNC}}^{\prime (1)} : \boldsymbol{\nabla }\boldsymbol{u}^{(2)} {\rm d}V=\,\boldsymbol{\varOmega }^{(2)} \boldsymbol{\cdot }\int _{S}\boldsymbol{x} \times \big \{ \big (\boldsymbol{\sigma }_{\textit{NNC}}^{\prime (1)} - p^{\prime (1)}_{G} \boldsymbol{I} \big ) \boldsymbol{\cdot }\boldsymbol{n} \big \} {\rm d}S . \end{equation}

Further, noting that the test velocity field (2.11)–(2.12) is linear in $\boldsymbol{\varOmega }^{(2)}$ , one may write $\boldsymbol{u}^{(2)} = \boldsymbol{U}^{(2)} \boldsymbol{\cdot }\boldsymbol{\varOmega }^{(2)}$ , $\mathcal{\boldsymbol{L}}^{(2)} = \boldsymbol{L}^{(2)} \boldsymbol{\cdot }\boldsymbol{\varOmega }^{(2)}$ and $\boldsymbol{\sigma }^{(2)} = \boldsymbol{\varSigma }^{(2)} \boldsymbol{\cdot }\boldsymbol{\varOmega }^{(2)}$ , where $\boldsymbol{U}^{(2)}$ and $\boldsymbol{L}^{(2)}$ are second-order tensors, and $ \boldsymbol{\varSigma }^{(2)}$ a third-order tensor, dependent on position and the spheroid aspect ratio. All three tensors are known in closed form, see Dabade, Marath & Subramanian (Reference Dabade, Marath and Subramanian2016). Using (2.14) and (2.15), and the fact that $\boldsymbol{\varOmega }^{(2)}$ is arbitrary but for being transverse to $\boldsymbol{p}$ , (2.13) takes the final form

(2.16) \begin{equation} \begin{aligned} \boldsymbol{\varOmega }^{(1)} \boldsymbol{\cdot }\boldsymbol{L}^{(2)} &= \boldsymbol{\varGamma } \boldsymbol{:} \int _{Sp} \boldsymbol{x} \big (\boldsymbol{\varSigma }^{(2)} \boldsymbol{\cdot }\boldsymbol{n} \big ) \,{\rm d}S\\ &\quad+{\textit{De}}\,\left (\int _{S}\boldsymbol{x} \times \big \{ \big (\boldsymbol{\sigma }_{\textit{NNC}}^{\prime (1)} - p^{\prime (1)}_{G} \boldsymbol{I} \big ) \boldsymbol{\cdot }\boldsymbol{n} \big \} {\rm d}S + \int _{V} \boldsymbol{\sigma }_{\textit{NNQ}}^{\prime (1)} : \boldsymbol{\nabla }\boldsymbol{U}^{(2)} {\rm d}V \right )\!. \end{aligned} \end{equation}

The first term in (2.16) is the Stokesian torque that leads to the Jeffery orbits depicted in figure 1, while the second $O(De)$ term is the viscoelastic torque contribution. To $O(De)$ , the fields involved in the associated surface and volume integrals may be replaced by their Stokesian analogues (the resulting integrals are readily shown to be convergent), leading to additive contributions from the corotational and quadratic components. Although not analysed here, for the same reason, inertial and viscoelastic contributions to the spheroid angular velocity are also additive in the small $Re$ , ${\textit{De}}$ limit, so the inertioelastic trajectory topology is only a function of $El$ .

Figure 2.

The body-fixed ( $\textit{XYZ}$ ) and space-fixed ( $X^\prime Y^\prime Z^\prime$ ) coordinate systems for a spheroid in a simple shear flow; the $Y$ -axis being constrained to lie in the $X^\prime Y^\prime$ -plane at all times. The polar angle between the spheroid symmetry axis and ambient vorticity ( $\theta _j$ ), and the dihedral angle between the flow-vorticity and orientation-vorticity planes ( $\phi _{j}$ ) define the spheroid orientation in the space-fixed reference frame.

Figure 2. Long description

The diagram illustrates the body-fixed and space-fixed coordinate systems for a spheroid in a simple shear flow. The body-fixed coordinate system is represented by unprimed axes (X,Y, Z), while the space-fixed coordinate system is represented by primed axes (X', Y', Z'); the Y-axis being constrained to lie in the X' Y' -plane all the times. The polar angle (θj) is the angle between the spheroid’s symmetry axis and the ambient vorticity. The dihedral angle (φj) is the angle between the flow-vorticity and orientation-vorticity planes. These angles define the spheroid’s orientation in the space-fixed reference frame.

Following Dabade et al. (Reference Dabade, Marath and Subramanian2016), in § 2.2 below, we write down the disturbance velocity and pressure fields, that determine the various fields involved in (2.16), in terms of superpositions of the appropriate vector spheroidal harmonics. The latter are defined in a body-fixed coordinate system ( $\textit{XYZ}$ ) with the $Z$ -axis instantaneously aligned with the spheroid symmetry axis, and rotating with it such that the $Y$ -axis lies within the flow-gradient ( $X\prime Y\prime$ ) plane at all times, see figure 2. As a result, the angular velocity components of the coordinate system, orthogonal to $\boldsymbol{p}$ , coincide with those of the spheroid, both being given by $\boldsymbol{\varOmega }_{\textit{jeff}}\boldsymbol{\cdot }({\boldsymbol I} - {\boldsymbol p}{\boldsymbol p}) \equiv (\dot {\phi }_{j} \sin \theta _{j},\dot {\theta }_{j},0)$ ; $\dot {\theta }_j$ and $\dot {\phi }_j$ are defined below by (2.38). But, the spin about the $Z$ -axis ( $\dot {\phi }_{j} \cos \theta _{j}$ ) differs from the actual spheroid spin ( $({1}/{2})\boldsymbol \omega \boldsymbol{\cdot }\boldsymbol p = -({{1}/{2}})\cos \theta _j$ , $\boldsymbol \omega =-\boldsymbol{1}_{z\prime}$ being the ambient vorticity vector). Since the unsteady terms in (2.9) and (2.14) need to be evaluated in the aforementioned coordinate system, we rewrite the material-derivative combination in the said equations as

(2.17) \begin{align} \frac {\partial \boldsymbol{e}^{\prime (1)}}{\partial t} +\boldsymbol{u}^{(1)} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{e}^{\prime (1)} &= \left (\frac {\partial \boldsymbol{e}^{\prime (1)}}{\partial t}+ ( \boldsymbol{\varOmega }_{b} \times \boldsymbol{x} ) \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{e}^{\prime (1)} + \boldsymbol{\varOmega }_{\textit{bt}} \boldsymbol{\cdot }\boldsymbol{e}^{\prime (1)} + \big(\boldsymbol{\varOmega }_{\textit{bt}} \boldsymbol{\cdot }\boldsymbol{e}^{\prime (1)}\big)^{\dagger } \right ) \nonumber\\ &\quad +\boldsymbol{u}^{(1)} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{e}^{\prime (1)} - ( \boldsymbol{\varOmega }_{b} \times \boldsymbol{x} ) \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{e}^{\prime (1)}+ \boldsymbol{e}^{\prime (1)}\boldsymbol{\cdot }\boldsymbol{\varOmega }_{\textit{bt}}+\big(\boldsymbol{e}^{\prime (1)} \boldsymbol{\cdot }\boldsymbol{\varOmega }_{\textit{bt}}\big)^{\dagger }\nonumber\\ &= \left (\frac {\partial \boldsymbol{e}^{\prime (1)}}{\partial t}\right )_{r} +\big (\boldsymbol{u}^{(1)} - ( \boldsymbol{\varOmega }_{b} \times \boldsymbol{x} ) \big) \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{e}^{\prime (1)}+ \boldsymbol{e}^{\prime (1)}\boldsymbol{\cdot }\boldsymbol{\varOmega }_{\textit{bt}}\nonumber\\ &\quad +\big(\boldsymbol{e}^{\prime (1)} \boldsymbol{\cdot }\boldsymbol{\varOmega }_{\textit{bt}}\big)^{\dagger }\!, \end{align}

(2.18) \begin{align} \frac {\partial p^{\prime (1)}}{\partial t}+\boldsymbol{u}^{(1)} \boldsymbol{\cdot }\boldsymbol{\nabla }{\!p}^{\prime (1)} &= \left (\frac {\partial p^{\prime (1)}}{\partial t} + ( \boldsymbol{\varOmega }_{b} \times \boldsymbol{x} ) \boldsymbol{\cdot }\boldsymbol{\nabla }{\!p}^{\prime (1)} \right )+\boldsymbol{u}^{(1)} \boldsymbol{\cdot }\boldsymbol{\nabla }{\!p}^{\prime (1)} - ( \boldsymbol{\varOmega }_{b} \times \boldsymbol{x} ) \boldsymbol{\cdot }\boldsymbol{\nabla }{\!p}^{\prime (1)}\nonumber\\ &= \left (\frac {\partial p^{\prime (1)}}{\partial t} \right )_{r}+\big ( \boldsymbol{u}^{(1)} - ( \boldsymbol{\varOmega }_{b} \times \boldsymbol{x} ) \big ) \boldsymbol{\cdot }\boldsymbol{\nabla }{\!p}^{\prime (1)}, \end{align}

where $\boldsymbol{\varOmega }_{b}=(-\dot {\phi }_{j} \sin \theta _{j},\dot {\theta }_{j},\dot {\phi }_{j} \cos \theta _{j} )$ in (2.17)–(2.18) is the coordinate system angular velocity mentioned above. The terms with suffix $r$ in (2.17) and (2.18) denote the rates of change in the body-fixed frame. While the expression for the transformation of the scalar pressure field, to the said frame, is obvious, that for the tensor field ( $\boldsymbol{e}^{\prime (1)}$ ) involves a Jaumann-like contribution, $\boldsymbol{e}^{\prime (1)}\boldsymbol{\cdot }\boldsymbol{\varOmega }_{\textit{bt}}+(\boldsymbol{e}^{\prime (1)} \boldsymbol{\cdot }\boldsymbol{\varOmega }_{\textit{bt}})^{\dagger }$ , where $\boldsymbol \varOmega _{\textit{bt}}=({{\rm d} \boldsymbol{R(t)}}/{{\rm d}t}) \boldsymbol{\cdot }{\boldsymbol R}^{\dagger }(t)$ is the antisymmetric tensor corresponding to the pseudovector $\boldsymbol{\varOmega }_b$ , with

\begin{equation*} \boldsymbol{R}(t)= \left[\begin{matrix} \cos (\theta _{j}) \cos \phi _{j} & \cos \theta _{j} \sin \phi _{j} & -\sin \theta _{j} \\[3pt] -\sin \phi _{j}& \cos \phi _{j} & 0\\[3pt] \sin \theta _{j} \cos \phi _{j}& \sin \theta _{j} \sin \phi _{j} & \cos \theta _{j}\\ \end{matrix} \right] \end{equation*}

being the rotation tensor that connects the body-fixed and space-fixed frames. Using (2.17) and (2.18), one may write down the following final expressions for evaluation of the viscoelastic stress components and the Giesekus pressure field in (2.16):

(2.19) \begin{align} \boldsymbol{\sigma }^{\prime (1)}_{\textit{NNC}}&= 2 \epsilon \left ( \left (\frac {\partial \boldsymbol{e}^{\prime (1)}}{\partial t}\right )_{r} +\big (\boldsymbol{u}^{(1)} - ( \boldsymbol{\varOmega }_{b} \times \boldsymbol{x} ) \big) \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{e}^{\prime (1)}+ \boldsymbol{e}^{\prime (1)}\boldsymbol{\cdot }\boldsymbol{\varOmega }_{\textit{bt}}+\big(\boldsymbol{e}^{\prime (1)} \boldsymbol{\cdot }\boldsymbol{\varOmega }_{\textit{bt}}\big)^{\dagger }\right.\nonumber\\ &\quad\left.+\boldsymbol{w}^{\prime (1)} \boldsymbol{\cdot }\boldsymbol{e}^{\prime (1)}+\!\big(\!\boldsymbol{w}^{\prime (1)} \!\boldsymbol{\cdot }\boldsymbol{e}^{\prime (1)}\!\big)^{\dagger } \!+\boldsymbol{W} \!\boldsymbol{\cdot }\boldsymbol{e}^{\prime (1)}\!+\big(\!\boldsymbol{W} \boldsymbol{\cdot }\boldsymbol{e}^{\prime (1)}\!\big)^{\dagger }\!+\boldsymbol{w}^{\prime (1)} \boldsymbol{\cdot }\boldsymbol{E}\!+\!\big(\!\boldsymbol{w}^{\prime (1)} \boldsymbol{\cdot }\boldsymbol{E}\big)^{\dagger } \right )\!, \end{align}

(2.20) \begin{align} \boldsymbol{\sigma }_{\textit{NNQ}}^{\prime (1)}&=4 (1+\epsilon ) \big(\boldsymbol{e}^{\prime (1)} \boldsymbol{\cdot } \boldsymbol{e}^{\prime (1)}+\boldsymbol{e}^{\prime (1)} \boldsymbol{\cdot }\boldsymbol{E}+\boldsymbol{E} \boldsymbol{\cdot e}^{\prime (1)}\big ), \end{align}

(2.21) \begin{align} p_{G}^{\prime (1)}&=\epsilon \left ( \left ( \frac {\partial p^{\prime (1)}}{\partial t} \right )_{r}+\big ( \boldsymbol{u}^{(1)} - ( \boldsymbol{\varOmega }_{b} \times \boldsymbol{x} ) \big) \boldsymbol{\cdot }\boldsymbol{\nabla }{\!p}^{\prime (1)} + \boldsymbol{e}^{\prime (1)} : \boldsymbol{e}^{\prime (1)}+ 2 (\boldsymbol{e}^{\prime (1)} : \boldsymbol{E}) \right )\!. \end{align}

The time derivatives that appear in (2.19)–(2.21) are evaluated as $ ( {\partial }/{\partial \,t} )_{r}=\dot {\theta }_{j}( {\partial }/{\partial \theta _{j}})+\dot {\phi }_{j}( {\partial }/{\partial \phi _{j}})$ where, in keeping with the use of a Stokesian approximation for the viscoelastic integrand, $\dot {\theta }_j$ and $\dot {\phi }_j$ are evaluated using (2.38), which corresponds to the spheroid rotating along Jeffery orbits.

2.2. The Stokesian disturbance fields in spheroidal coordinates: Jeffery orbits at ${\textit{De}} = 0$

As already indicated, to $O(De)$ , the integrals in (2.16) may be evaluated using the Stokesian approximations for the disturbance fields in the actual problem. One therefore needs the Stokes disturbance fields associated with a torque-free spheroid in a simple shear flow ( $\boldsymbol{u}^{(1)}_{s}, p^{(1)}_s$ ), and that induced by a rotating spheroid in an otherwise quiescent fluid ( $\boldsymbol{u}^{(2)}, p^{(2)}$ ). The disturbance in the former case arises only due to the straining component of the ambient shear. Following Dabade et al. (Reference Dabade, Marath and Subramanian2016), one may decompose the ambient rate-of-strain tensor, given by $\boldsymbol{E}= \cos ^{2} \theta _{j} \sin \phi _{j} \cos \phi _{j} \boldsymbol{1}_{x}\boldsymbol{1}_{x} +{} ( {1}/{2})\cos \theta _{j} \cos 2 \phi _{j} (\boldsymbol{1}_{x}\boldsymbol{1}_{y}+\boldsymbol{1}_{y}\boldsymbol{1}_{x} )+ \sin \theta _{j} \cos \theta _{j} \sin \phi _{j} \cos \phi _{j} ( \boldsymbol{1}_{x}\boldsymbol{1}_{z} + \boldsymbol{1}_{z}\boldsymbol{1}_{x} )- \cos \phi _{j} {}\sin \phi _{j} \boldsymbol{1}_{y}\boldsymbol{1}_{y}+ ({1}/{2})\sin \theta _{j} \cos 2 \phi _{j} (\boldsymbol{1}_{y}\boldsymbol{1}_{z}+\boldsymbol{1}_{z}\boldsymbol{1}_{y} )+ \sin ^{2} \theta _{j} \sin \phi _{j} \cos \phi _{j} \boldsymbol{1}_{z}\boldsymbol{1}_{z}$ in the body-fixed frame, into the following five canonical components:

(2.22) \begin{equation} \begin{aligned} \boldsymbol{E} &=\begin{bmatrix} - \dfrac {1}{2} \sin ^{2} \theta _{j} \sin \phi _{j} \cos \phi _{j} & 0 & 0\\ 0 & -\dfrac {1}{2} \sin ^{2} \theta _{j} \sin \phi _{j} \cos \phi _{j} & 0\\ 0 & 0& \sin ^{2} \theta _{j} \sin \phi _{j} \cos \phi _{j}\\ \end{bmatrix}\\ &\quad+\begin{bmatrix} \dfrac {1}{2} (\cos ^{2} \theta _{j}+1) \sin \phi _{j} \cos \phi _{j} & 0 & 0\\ 0 & -\dfrac {1}{2} \left(\cos ^{2} \theta _{j} +1\right) \sin \phi _{j} \cos \phi _{j} & 0\\ 0 & 0& 0\\ \end{bmatrix}\\ &\quad+\begin{bmatrix} 0 & \dfrac {1}{2} \cos \theta _{j} \left(\cos ^{2} \phi _{j}-\sin ^{2} \phi _{j}\right) & 0\\ \dfrac {1}{2} \cos \theta _{j} \left(\cos ^{2} \phi _{j}-\sin ^{2} \phi _{j}\right) & 0 & 0\\ 0 & 0& 0\\ \end{bmatrix}\\ &\quad+\begin{bmatrix} 0 & 0 & \sin \theta _{j} \cos \theta _{j} \sin \phi _{j} \cos \phi _{j}\\ 0 & 0 & 0\\ \sin \theta _{j} \cos \theta _{j} \sin \phi _{j} \cos \phi _{j} & 0& 0\\ \end{bmatrix}\\ &\quad+\begin{bmatrix} 0 & 0 & 0\\ 0 & 0 & \dfrac {1}{2} \sin \theta _{j}\left(\cos ^{2} \phi _{j} -\sin ^{2} \phi _{j}\right)\\ 0 & \dfrac {1}{2} \sin \theta _{j}\left(\cos ^{2} \phi _{j} -\sin ^{2} \phi _{j}\right)& 0\\ \end{bmatrix}. \end{aligned} \end{equation}

The component matrices in (2.22) correspond, respectively, to (i) an axisymmetric extensional flow (uniaxial or biaxial, depending on the sign of $\sin \phi _{j}\cos \phi _{j}$ ) along the spheroidal symmetry axis ( ${\boldsymbol 1}_z$ ), with amplitude $\sin ^{2}\theta _{j}\sin \phi _{j}\cos \phi _{j}$ ; (ii,iii) a pair of planar extensions in the transverse ( $XY$ ) plane, one obtained from the other by a $45^{\circ }$ rotation about the symmetry axis, and with amplitudes $(\cos ^{2}\theta _{j}+1)\sin \phi _{j}\cos \phi _{j}/2$ and $\cos \theta _{j}(\cos ^{2}\phi _{j}-\sin ^{2}\phi _{j})/2$ ; and (iv,v) a pair of longitudinal extensions in the $XZ$ and $YZ$ planes, with amplitudes $\sin \theta _{j}\cos \theta _{j}\sin \phi _{j}\cos \phi _{j}$ and $\sin \theta _{j}(\cos ^{2}\phi _{j}-\sin ^{2}\phi _{j})/2$ . The disturbance velocity ( $\boldsymbol{u}_{1s}-\boldsymbol{u}_{5s}$ ) and pressure fields ( $p_{1s}-p_{5s}$ ), corresponding to each of the component flows, are given by

(2.23) \begin{align} &\boldsymbol{u}_{1s}=\frac {-{\rm d} \bar {\xi }_{0}}{Q^{1}_{1}(\xi _{0})-\xi _{0} Q^{1}_{2}(\xi _{0})} \big(\sin ^{2} \theta _{j} \sin \phi _{j} \cos \phi _{j} \big) \boldsymbol{S}_{2,0}^{(3)}, \end{align}

(2.24) \begin{align} &\boldsymbol{u}_{2s}=\frac {-{\rm d} \bar {\xi }_{0}}{3Q^{1}_{1}(\xi _{0})-\xi _{0} Q^{1}_{2}(\xi _{0})} \big(\sin \phi _{j} \cos \phi _{j}\big(1+\cos ^{2} \theta _{j}\big) \big) \big(\boldsymbol{S}_{2,2}^{(3)}+\boldsymbol{S}_{2,-2}^{(3)} \big), \end{align}

(2.25) \begin{align} &\boldsymbol{u}_{3s}=\frac {i {\rm d} \bar {\xi }_{0}}{3Q^{1}_{1}(\xi _{0})-\xi _{0} Q^{1}_{2}(\xi _{0})} \big(\cos \theta _{j} \cos 2 \phi _{j} \big) \big(\boldsymbol{S}_{2,2}^{(3)}-\boldsymbol{S}_{2,-2}^{(3)} \big), \end{align}

(2.26) \begin{align} &\boldsymbol{u}_{4s}=\frac {2 {\rm d} \xi _{0} \bar {\xi }_{0}}{Q^{1}_{2}(\xi _{0})\left(2 \xi _{0}^2-1\right)} \big(\sin \theta _{j} \cos \theta _{j} \sin \phi _{j} \cos \phi _{j} \big) \big(\boldsymbol{S}_{2,1}^{(3)}-\boldsymbol{S}_{2,-1}^{(3)} \big), \end{align}

(2.27) \begin{align} &\boldsymbol{u}_{5s}=\frac {-i {\rm d} \xi _{0} \bar {\xi }_{0}}{Q^{1}_{2}(\xi _{0})(2 \xi _{0}^2-1)} \big(\sin \theta _{j} \big( \cos ^{2} \phi _{j}-\sin ^{2} \phi _{j}\big) \big) \big(\boldsymbol{S}_{2,1}^{(3)}+\boldsymbol{S}_{2,-1}^{(3)} \big), \end{align}

(2.28) \begin{align} &p_{1s}=\frac {-2 {\rm d} \bar {\xi }_{0}}{Q^{1}_{1}(\xi _{0})-\xi _{0} Q^{1}_{2}(\xi _{0})} \big(\sin ^{2} \theta _{j} \sin \phi _{j} \cos \phi _{j}\big) D_{3}F_{1}^{0}, \end{align}

(2.29) \begin{align} &p_{2s}=\frac {-2{\rm d} \bar {\xi }_{0}}{3Q^{1}_{1}(\xi _{0})-\xi _{0} Q^{1}_{2}(\xi _{0})} \big(\sin \phi _{j} \cos \phi _{j}\big(1+\cos ^{2} \theta _{j}\big)\big) \big(D_{3} F_{1}^{2}+D_{3} F_{1}^{-2}\big), \end{align}

(2.30) \begin{align} &p_{3s}=\frac { 2i {\rm d} \bar {\xi }_{0}}{3Q^{1}_{1}(\xi _{0})-\xi _{0} Q^{1}_{2}(\xi _{0})} \big(\cos \theta _{j} \cos 2 \phi _{j}\big) \big(D_{3} F_{1}^{2}-D_{3} F_{1}^{-2}\big), \end{align}

(2.31) \begin{align} &p_{4s}=\frac {4 {\rm d} \xi _{0} \bar {\xi }_{0}}{Q^{1}_{2}(\xi _{0})\left(2 \xi _{0}^2-1\right)} \big(\sin \theta _{j} \cos \theta _{j} \sin \phi _{j} \cos \phi _{j}\big) \big(D_{3} F_{1}^{1}-D_{3} F_{1}^{-1}\big), \end{align}

(2.32) \begin{align} &p_{5s}=\frac {- 2 i {\rm d} \xi _{0} \bar {\xi }_{0}}{Q^{1}_{2}(\xi _{0})\left(2 \xi _{0}^2-1\right)} \big(\sin \theta _{j} \big( \cos ^{2} \phi _{j}-\sin ^{2} \phi _{j}\big) \big) \big(D_{3} F_{1}^{1}+D_{3} F_{1}^{-1}\big). \end{align}

The test velocity and pressure fields, corresponding to transverse rotations about the $X$ and $Y$ axes, are given by

(2.33) \begin{align} \boldsymbol{u}^{(2)}_{\textit{sx}}&=\frac {i {\rm d} (2 \xi _{0}^{2}-1)}{2 \xi _{0} Q^{0}_{1}(\xi _{0})-\bar {\xi }_{0} Q^{1}_{1}(\xi _{0})}\big(\boldsymbol{S}^{(2)}_{1,1}-\boldsymbol{S}^{(2)}_{1,-1}\big) \nonumber\\&\quad +\frac {i {\rm d} (\xi _{0} Q^{1}_{1}(\xi _{0})+2\bar {\xi }_{0} Q^{0}_{1}(\xi _{0}) ) }{Q^{1}_{2}(\xi _{0})(2 \xi _{0} Q^{0}_{1}(\xi _{0})-\bar {\xi }_{0} Q^{1}_{1}(\xi _{0}))}\big(\boldsymbol{S}^{(3)}_{2,1}+\boldsymbol{S}^{(3)}_{2,-1}\big), \end{align}

(2.34) \begin{align} \boldsymbol{u}^{(2)}_{\textit{sy}}&=\frac {{\rm d} (2 \xi _{0}^{2}-1)}{2 \xi _{0} Q^{0}_{1}(\xi _{0})-\bar {\xi }_{0} Q^{1}_{1}(\xi _{0})}\big(\boldsymbol{S}^{(2)}_{1,1}+\boldsymbol{S}^{(2)}_{1,-1}\big) \nonumber\\&\quad +\frac { {\rm d} (\xi _{0} Q^{1}_{1}(\xi _{0})+2\bar {\xi }_{0} Q^{0}_{1}(\xi _{0}) ) }{Q^{1}_{2}(\xi _{0})(2 \xi _{0} Q^{0}_{1}(\xi _{0})-\bar {\xi }_{0} Q^{1}_{1}(\xi _{0}))}\big(\boldsymbol{S}^{(3)}_{2,1}-\boldsymbol{S}^{(3)}_{2,-1}\big), \end{align}

(2.35) \begin{align} p^{(2)}_{\textit{sx}}&=\frac {2 i {\rm d} (\xi _{0} Q^{1}_{1}(\xi _{0})+2\bar {\xi }_{0} Q^{0}_{1}(\xi _{0}) ) }{Q^{1}_{2}(\xi _{0})\left(2 \xi _{0} Q^{0}_{1}(\xi _{0})-\bar {\xi }_{0} Q^{1}_{1}(\xi _{0})\right)}\big(D^{1}_{2}+D^{-1}_{2}\big), \end{align}

(2.36) \begin{align} p^{(2)}_{\textit{sy}}&=\frac { 2 {\rm d} (\xi _{0} Q^{1}_{1}(\xi _{0})+2\bar {\xi }_{0} Q^{0}_{1}(\xi _{0}) ) }{Q^{1}_{2}(\xi _{0})\left(2 \xi _{0} Q^{0}_{1}(\xi _{0})-\bar {\xi }_{0} Q^{1}_{1}(\xi _{0})\right)}\big(D^{1}_{2}-D^{-1}_{2}\big). \end{align}

The $\boldsymbol S^{(q)}_{t,s}$ ( $q = 1-3$ ) in the above expressions are vector spheroidal harmonics, with $D_t^s$ denoting the pressure contribution associated with $\boldsymbol S^{(3)}_{t,s}$ . The detailed expressions for these harmonics in spheroidal coordinates $(\xi ,\eta ,\phi )$ , and their connections to the vector spherical harmonics in the classical Lamb expansion (chapter $4$ in Kim & Karrila (Reference Kim and Karrila1991)) have been discussed in earlier papers (Dabade et al. Reference Dabade, Marath and Subramanian2015, Reference Dabade, Marath and Subramanian2016; Marath & Subramanian Reference Marath and Subramanian2018). Here, $\xi _0$ in (2.23)–(2.36) is the coordinate label for the spheroid, with $d$ being half the interfoci distance (focal radius) for the prolate (oblate) spheroid; $\bar {\xi }_0 =\sqrt {\xi _0^2-1}$ . It must be said that the above expressions are based on the pioneering results obtained by Kushch and coworkers (Kushch Reference Kushch1995, Reference Kushch1997, Reference Kushch1998; Kushch & Sangani Reference Kushch and Sangani2000; Kushch & Sevostianov Reference Kushch and Sevostianov2004) in the context of linear elasticity.

In the Newtonian limit, corresponding to ${\textit{De}}=0$ , (2.16) reduces to

(2.37) \begin{equation} \boldsymbol{\varOmega }_{\textit{jeff}} \boldsymbol{\cdot }\boldsymbol{L}^{(2)} = \boldsymbol{\varGamma } \boldsymbol{:} \int _{Sp} \boldsymbol{x} \big (\boldsymbol{\varSigma }^{(2)} \boldsymbol{\cdot }\boldsymbol{n} \big ) \,{\rm d}S=\boldsymbol{\varOmega }^{(1)} \boldsymbol{\cdot }\boldsymbol{L}^{(2)}, \end{equation}

where $\boldsymbol{\varOmega }_{\textit{jeff}}$ is the Jeffery angular velocity. On evaluating the surface integral on the right-hand side using the known expression for $\boldsymbol{\varSigma }^{(2)}$ in spheroidal coordinates (Dabade Reference Dabade2009; Dabade et al. Reference Dabade, Marath and Subramanian2016), one obtains the components of $\boldsymbol{\varOmega }_{\textit{jeff}}$ , orthogonal to $\boldsymbol{p}(\equiv \boldsymbol{1}_{z}) =\sin \theta _{j} \cos \phi _{j} \boldsymbol{1}_{x^{\prime }} +\sin \theta _{j} \sin \phi _{j} \boldsymbol{1}_{y^{\prime }} +\cos \theta _{j} \boldsymbol{1}_{z^{\prime }}$ , as

(2.38) \begin{equation} \dot {\phi }_{j}=-\frac {1}{2}+\frac {(\kappa ^2 -1)}{2 (\kappa ^2+1)} \cos 2 \phi _{j}, \qquad \dot {\theta }_{j} =\frac {(\kappa ^2 -1)}{4 (\kappa ^2+1)} \sin 2 \phi _{j} \sin 2 \theta _{j}. \end{equation}

The above equations lead to the Jeffery orbits in figure 1. The orbits are spherical ellipses; the ellipse orientation and shape as a function of $\kappa$ was described in § 1. Note that, in his original paper, Jeffery (Reference Jeffery1922) used ellipsoidal harmonics to first obtain the angular velocity of a triaxial ellipsoid, and then considered the simpler expression that results for a spheroid.

In the absence of inertia, viscoelasticity, Brownian motion or interparticle interactions, a spheroid is confined to the Jeffery orbit, determined by its initial orientation, for all time. This is readily seen in $(C,\tau )$ coordinates defined by

(2.39) \begin{equation} C=\tan \theta _{j} \frac {\sqrt {\kappa ^2 \sin ^2 \phi _{j}+\cos ^2 \phi _{j}}}{\kappa }; \quad \tau =\tan ^{-1}\left (\frac {1}{\kappa \tan \phi _{j} } \right )\!. \end{equation}

Here, $C\,\in \,[0,\infty )$ is the orbit constant, with $C=0$ and $C=\infty$ corresponding to the spinning (log-rolling) and tumbling modes, respectively. All finite non-zero $C$ correspond to three-dimensional precessional motions termed the kayaking modes; $\tau$ is the phase along a given orbit. Equation (2.38), written in $(C,\tau )$ coordinates, take the simple forms

(2.40) \begin{equation} \frac {{\rm d} C}{{\rm d} t}=0, \qquad \frac {{\rm d} \tau }{{\rm d} t}=\frac {\kappa }{\kappa ^2+1}, \end{equation}

where, unlike $\phi _{j}$ , $\tau$ is seen to change at a uniform albeit $\kappa$ -dependent rate. The period of rotation, $T_{\textit{jeff}}= 2\pi (\kappa ^2+1)/\kappa$ , is independent of $C$ . It is worth mentioning that, for inelastic fluids with a shear-rate-dependent viscosity, although reversibility constraints still lead to a spheroid moving along closed orbits, $T_{\textit{jeff}}$ is nevertheless a function of $C$ . For shear-thinning fluids, for instance, one expects the greater disturbance induced by a tumbling spheroid to lead to a smaller effective viscosity, and thereby faster rotation, compared with the kayaking modes. That $T_{\textit{jeff}}$ is a increasing function of $C$ for power-law fluids with $n \lt 1$ ( $n$ being the power-law index) has been confirmed in recent calculations (Abtahi & Elfring Reference Abtahi and Elfring2019).

In the context of the results presented below, $\kappa =\xi _{0}/\bar {\xi }_0$ for a prolate spheroid, $\xi _{0}$ being the coordinate label mentioned above in the context of the expressions for the disturbance fields. In what follows, we will use both $\kappa$ and $\xi _0$ to refer to the spheroid aspect ratio; $\xi _0 \rightarrow \infty \,(\kappa \rightarrow 1)$ and $\xi _0 \rightarrow 1\,(\kappa \rightarrow \infty )$ correspond to the spherical and slender-fibre limits. One may rewrite (2.38) in terms of $\xi _0$ , in which case the orbit equations for an oblate spheroid are obtained, from those for a prolate one, by using the transformation $\xi _{0} \leftrightarrow i \bar {\xi _{0}}$ , $d \leftrightarrow -i d$ in the dimensional form of the equations. This transformation is used below, for non-zero ${\textit{De}}$ , to obtain the oblate orbital drift from the prolate one; the aspect ratio on the oblate side is given by $\kappa = \bar {\xi }_0/\xi _0$ .

2.3. The rotation rate in presence of weak viscoelastic effects: ${\textit{De}}\ll 1$

Linearity of the Stokes equations translating to a linear dependence of the spheroid rotation rate $\dot {\boldsymbol p}$ on the ambient velocity gradient, the unit length constraint ( $\boldsymbol{p}\boldsymbol{\cdot }\dot {\boldsymbol p}=0$ ) and the requirement of an affine response in the limit of solid body rotation ( $\dot {\boldsymbol p} = \boldsymbol{p} \boldsymbol{\cdot }{\boldsymbol W}$ ), immediately leads to the general form $\dot {\boldsymbol p} = \boldsymbol{p} \boldsymbol{\cdot }\boldsymbol{W} + B(\kappa )\left \{ \boldsymbol{E} \boldsymbol{\cdot }\boldsymbol{p} - (\boldsymbol{E}:\boldsymbol{p}\!\boldsymbol{p})\boldsymbol{p}\right \}$ , for ${\textit{De}} = 0$ , in an ambient linear flow. The detailed calculations show that the Bretherton constant, $B(\kappa ) ={\kappa ^2-1}/{\kappa ^2+1}$ ; the components of the resulting equation in spherical polar coordinates, for simple shear flow, are given in (2.38). The regular nature of the viscoelastic correction, and the quadratic dependence of the second-order fluid stress on the velocity gradient, imply that the $O(De)$ correction to the above rate of change must be a quadratic function of $\boldsymbol{E}$ and $\boldsymbol{W}$ . The other two constraints above then lead to the following form for the rotation rate, to $O(De)$ in an ambient linear flow:

(2.41) \begin{align} \begin{aligned} \boldsymbol{\dot {p}}&=\boldsymbol{p} \boldsymbol{\cdot }\boldsymbol{W} + B(\kappa ) \left \{ \boldsymbol{E} \boldsymbol{\cdot }\boldsymbol{p} - (\boldsymbol{E}:\boldsymbol{p}\!\boldsymbol{p})\boldsymbol{p} \right \} + \textit{De} \big \{\beta _{1}(\xi _0,\epsilon )(\boldsymbol{E:pp})(\boldsymbol{I-pp}) \boldsymbol{\cdot }(\boldsymbol{E \boldsymbol{\cdot }p})\\ &\quad+\beta _{2}(\xi _0,\epsilon )(\boldsymbol{E:pp})\boldsymbol{W\boldsymbol{\cdot }p}+\beta _{3}(\xi _0,\epsilon )\boldsymbol{(I-pp)\boldsymbol{\cdot }W \boldsymbol{\cdot }(E\boldsymbol{\cdot }p)}\\ &\quad+\beta _{4}(\xi _0,\epsilon )(\boldsymbol{I-pp})\boldsymbol{\cdot E \boldsymbol{\cdot }}(\boldsymbol{E \boldsymbol{\cdot }p}) \big \}. \end{aligned} \end{align}

The purpose of a detailed calculation is to determine the four independent functions of $\xi _0$ and $\epsilon$ ; the dependence on $\epsilon$ is linear, mirroring the dependence in the original constitutive equation – see (2.5) and (2.6). The quadratic form in (2.41) is also valid for weak inertial effects, with $Re$ now being the proportionality factor, and the $\beta _i$ only being functions of $\xi _0$ . The closed-form expressions for these functions are implicit in the results of Dabade et al. (Reference Dabade, Marath and Subramanian2016) and Marath & Subramanian (Reference Marath and Subramanian2018), and have been tabulated in Appendix A.2.

Returning to the viscoelastic correction and the case of simple shear flow, the components of the $O(De)$ correction in (2.41), in terms of $\theta _j$ and $\phi _j$ , assume the following forms:

(2.42) \begin{align} (\dot {\phi }_{j})_{\textit{De}}&= \cos \phi _{j} \sin \phi _{j} \big (G_{1}(\xi _{0},\epsilon )+ G_{2}(\xi _{0},\epsilon )\cos 2 \theta _{j} + G_{3}(\xi _{0},\epsilon )\cos 2 \phi _{j}\nonumber\\ &\quad +G_{4}(\xi _{0},\epsilon ) \cos 2\theta _{j} \cos 2\phi _{j} \big ), \end{align}

(2.43) \begin{align} (\dot {\theta }_{j})_{\textit{De}}&= \cos \theta _{j} \sin \theta _{j} \big ( F_{1}(\xi _{0},\epsilon ) +F_{2}(\xi _{0},\epsilon ) \cos 2 \phi _{j} +F_{3}(\xi _{0},\epsilon ) \cos 2 \theta _{j} \nonumber\\ &\quad + F_{4}(\xi _{0},\epsilon ) \cos 4\phi _{j} +F_{5}(\xi _{0},\epsilon ) \cos \bigl(2\theta _{j}-4\phi _{j}\bigr)+ F_{6}(\xi _{0},\epsilon ) \cos \bigl(2\theta _{j}+4\phi _{j}\bigr) \big ), \end{align}

where the $F_{i}$ and $G_{i}$ are related to the $\beta _{i}$ above through the following relations: $F_{1}=( \beta _{1} +4 \beta _{4})/16$ , $F_{2}=-\beta _{3}/4$ , $F_{3}=-\beta _{1}/16$ , $F_{3}=F_{4}$ , $F_{5}=F_{6}=-F_{3}/2$ , $G_{1}=(\beta _{2}+2\beta _{3})/4$ , $G_{2}=-\beta _{2}/4$ and $G_{3}=-G_{4}=\beta _{1}/4$ . These functions have been obtained in closed form by analytically evaluating the surface integral involving $p\prime_G$ , and the volume integral involving $\boldsymbol{\sigma }^{\prime (1)}_{\textit{NNQ}}$ , in (2.16). The detailed expressions are given in Appendix A.1, where the corotational and quadratic contributions to each of these functions have been identified. Interestingly, the final expressions imply $\beta _2 = \beta _3 = 0$ , so that the original invariant form, (2.41), is independent of $\boldsymbol{W}$ , this being consistent with the general expression proposed by Brunn (Reference Brunn1977b ) (and evaluated by him for a tridumbbell). Thus, one may write (2.42)–(2.43) in the following simpler form:

(2.44) \begin{align} (\dot {\phi }_{j})_{\textit{De}}&= \frac {1}{2}G_{3}(\xi _{0},\epsilon ) \sin ^2\theta _j \sin 4 \phi _j, \end{align}

(2.45) \begin{align} (\dot {\theta }_{j})_{\textit{De}}&= \cos \theta _j \sin \theta _j \big ( F_{1}(\xi _{0},\epsilon )+F_{3}(\xi _{0}, \epsilon )\big ( 1+2 \sin ^2 \theta _j (\cos 4 \phi _j -1) \big ) \big ). \end{align}

Note that the above is not true for the inertial case, which ends up involving both $\boldsymbol{E}$ and $\boldsymbol{W}$ (Dabade et al. Reference Dabade, Marath and Subramanian2016; Marath & Subramanian Reference Marath and Subramanian2018). We also note that there is no contribution from the nonlinear terms of the corotational component in (2.19) and (2.21); in contrast, both linear and nonlinear terms of the quadratic component in (2.20) contribute to (2.42) and (2.43).

3.1. Multiple time scales analysis of spheroid rotation for ${\textit{De}} \ll 1$

Provided that the drift and Jeffery-rotation time scales are well separated, one may proceed via a multiple scales analysis with $t_1 = \dot {\gamma }t$ and $t_{2} = De\dot {\gamma }t$ as the fast and slow time variables (Bender & Orszag Reference Bender and Orszag1999; Subramanian & Brady Reference Subramanian and Brady2004). Expanding the time derivative as $d/{\rm d}t = d/{\rm d}\,t_1+{\textit{De}}\,d/{\rm d}t_2$ , one formally writes (3.1) in the form

(3.2) \begin{align} \frac {{\rm d}C}{{\rm d}t_1} +{\textit{De}} \frac {{\rm d}C}{{\rm d}t_2} &= {\textit{De}} \mathcal{H}(C,\tau ;\kappa ,\epsilon ), \end{align}

(3.3) \begin{align} &= {\textit{De}} \mathcal{H}\left(C,\frac {\kappa }{\kappa ^2+1}t_1;\kappa ,\epsilon \right)\!, \end{align}

where the expression for $\tau$ from (2.39) has been used in going from (3.2) to (3.3), keeping in mind that this represents a fast-time dependence. The function $\mathcal H$ is obtained by first expressing $\theta _j$ and $\phi _j$ in terms of $C$ and $\tau$ (this corresponds to inverting the $(C,\tau )$ -relations given by (2.39)), substituting these in (2.42) and (2.43), which may then be substituted in (3.1). Equation (3.3) suggests writing $C$ in terms of a two time scale expansion of the form (Marath & Subramanian Reference Marath and Subramanian2018; Banerjee & Subramanian Reference Banerjee and Subramanian2021) $C =C_{0}(t_1,t_2)+{\textit{De}}\, C_1(t_1,t_2)+O(De^2)$ , which gives

(3.4) \begin{equation} \frac {{\rm d} C_{0}}{{\rm d} t_1} = 0, \end{equation}

at leading order, so $C_{0} \equiv C_0(t_2)$ represents the slow drift due to viscoelasticity. At $O(De)$ , one obtains

(3.5) \begin{equation} \frac {{\rm d} C_{1}}{{\rm d} t_{1}} + \frac {{\rm d} C_{0}}{{\rm d} t_{2}} = \mathcal{H} \left (C_0, \frac {\kappa t_{1}}{\kappa ^2+1}; \kappa , \epsilon \right )\!, \end{equation}

where $C_1$ accounts for the small amplitude fast fluctuations about the aforementioned drift. To determine $C_{0}$ , we integrate (3.5), with $t_2$ fixed, over $t_1\,\in \,[0,T_{\textit{jeff}}]$ , corresponding to a single fast-time period (a single Jeffery rotation). In doing so, $C_1$ is constrained to be a periodic function of $t_{1}$ , which is equivalent to eliminating secular growth on the fast time scale. This leads to the first term in (3.5) integrating out to zero, yielding

(3.6) \begin{equation} \begin{aligned} \frac {{\rm d} C_{0}}{{\rm d} t_{2}} &= \frac {1}{T_{\textit{jeff}}} \int _{0}^{T_{\textit{jeff}}} \mathcal{H} \left (C_0, \frac {\kappa t_{1}}{\kappa ^2+1}; \kappa , \epsilon \right ) {\rm d} t_{1} \\ &= \frac {1}{2 \pi } \int _{0}^{2\pi } \mathcal{H} (C_0, \tau ; \kappa , \epsilon ) {\rm d}\tau . \end{aligned} \end{equation}

Using (3.6) in (3.5) gives the following expression for $C_{1}$ :

(3.7) \begin{equation} C_{1} = \frac {T_{\textit{jeff}}}{2 \pi } \left ( \int _{0}^{\tau } \mathcal{H} (C_0,\tau ^{\prime }; \kappa , \epsilon ) {\rm d} \tau ^{\prime } - \frac {\tau }{2\pi }\int _{0}^{2\pi } \mathcal{H} (C_0, \tau ; \kappa , \epsilon ) {\rm d}\tau \right )\!, \end{equation}

where $C_0$ is assumed constant during the $\tau$ -integration. Together, (3.6) and (3.7) provide a reduced one-dimensional description of the effect of viscoelasticity. The autonomous ordinary differential equation governing $C_0$ may be solved numerically to obtain the slow drift of the orbit constant with time. This slow variation is then used in (3.7), after the $\tau$ -integration, to obtain the fast orbit-constant fluctuations that accompany the drift.

Figures 3(a) and 3(b) compare the orbit constant prediction ( $C_0 +{\textit{De}} C_1)$ , from solving (3.6) and (3.7), with that obtained from directly solving the equations for $\theta _j$ and $\phi _j$ , to $O(De)$ , and calculating $C$ using (2.39). This is done for ${\textit{De}}=0.02$ and $0.06$ for prolate spheroid with $\kappa =4$ and $\epsilon = -2$ , and for ${\textit{De}}=0.05$ and $0.1$ for prolate with $\kappa =4$ and $\epsilon =-0.01$ . In all cases, the spheroid starts from the initial orientation $(\theta _{j},\phi _j) \equiv (23\pi /48,0)$ , which corresponds to $C(0) =3.814$ for $\kappa =4$ . The main plots show $C$ decreasing to zero (figure 3 a), and increasing to infinity (figure 3 b), on a time scale of $O(De^{-1}\dot {\gamma }^{-1})$ , for the two values of $\epsilon$ considered. This is consistent with the organization of the small- ${\textit{De}}$ orientation dynamics regimes on the $\kappa -\epsilon$ plane, shown later in § 3.3, where the spinning mode is the only stable mode for $\epsilon = -2$ , in the interval $1 \lt \kappa \leqslant 9.324$ , with the tumbling mode alone being stable for $\epsilon = -0.01$ , in the interval $1 \lt \kappa \leqslant 4.050$ . The inset views in figures 3(a) and 3(b) show that the fluctuations in $C$ are also captured accurately by the multiple scales analysis, for small ${\textit{De}}$ ; as expected, the amplitude of these fluctuations scales as ${\textit{De}}$ .

Figure 3.

The orbital coordinate $C$ as a function of time $t$ , starting from $(\theta _{j},\phi _{j})=(23 \pi /48,0)$ , for a spheroid with aspect ratio $\kappa =4$ for (a) $\epsilon =-2$ and (b) $\epsilon =-0.01$ . The calculations are based on (i) $C$ obtained by numerically integrating the governing equations of $\theta _{j}$ and $\phi _{j}$ (solid line), (ii) $C_0$ obtained from solving (3.6) (dash–dotted line) and (iii) the drift together with the fluctuations, $C_0 +{\textit{De}} C_1$ (dashed line).

Figure 3. Long description

Two line graphs depict the orbital coordinate as a function of time for a spheroid with aspect ratio 4 for different calculations. The x-axis represents the dimensionless time, and the y-axis represents the normalized orbital coordinate. The graphs compare three different methods of calculation: numerically integrating the primitive governing equations, solving the leading order drift equation, and considering the drift with fluctuations. Each method is represented by a different line style: solid, dash-dotted, and dashed. The graphs show how the orbital coordinate evolves over time for each method, highlighting differences in the results. The inset in each graph provides a zoomed-in view of a specific time range to show more detailed differences between the methods. All values are approximated.

3.2. Effect of viscoelasticity based on the orbital drift interpretation

Although we have solved (3.6) for the orbit constant evolution, characterizing the small- ${\textit{De}}$ trajectory topology does not necessarily require solving the said ordinary differential equation. The right-hand side of (3.6) may be directly interpreted as the average rate of change of $C$ (scaled by ${\textit{De}}$ ) over a single Jeffery period, in the neighbourhood of a Jeffery orbit with orbit constant $C$ , and is therefore denoted as $ {\Delta C(C;\kappa ,\epsilon )}/{De\,T_{\textit{jeff}}}$ . Here, $\Delta C$ may also be regarded as the pitch of a tightly spiralling small- ${\textit{De}}$ trajectory, measured in units of $C$ . From the right-hand side of (3.6), one obtains

(3.8) \begin{align} \Delta C &= {\textit{De}} \int _{0}^{T_{\textit{jeff}}} \mathcal{H} \left (C_0, \frac {\kappa t_{1}}{\kappa ^2+1}; \kappa , \epsilon \right ) {\rm d} t_{1} \nonumber \\ &=\textit{De}\, C \frac {\kappa ^2+1}{\kappa } \big \{ \big ( I_{1} \,F_{1}(\xi _{0},\,\epsilon ) + I_{2} \,F_{2}(\xi _{0},\,\epsilon )+ I_{3}\, F_{3}(\xi _{0},\,\epsilon ) + I_{4} \, F_{4}(\xi _{0},\,\epsilon )\nonumber \\ &\quad + I_{5} \,F_{5}(\xi _{0},\,\epsilon ) + I_{6}\, F_{6}(\xi _{0},\,\epsilon ) \big ) + \big ( J_{1} \,G_{1}(\xi _{0},\,\epsilon ) + J_{2}\,G_{2}(\xi _{0},\,\epsilon )+ J_{3}\, G_{3}(\xi _{0},\,\epsilon ) \nonumber \\ &\quad + J_{4} \, G_{4}(\xi _{0},\,\epsilon ) \big )\big \}, \end{align}

where the $I_1-I_6$ and $J_1-J_4$ are functions of $C$ and $\kappa$ , that result from integrating the trigonometric functions in (2.42) and (2.43) over $\tau$ . These expressions have been given in Dabade et al. (Reference Dabade, Marath and Subramanian2016), and are tabulated in Appendix A.3 for convenience.

Figure 4.

The drift due to elasticity as characterized by the normalized change in the orbit constant in a single Jeffery period, ${\textit{De}}^{-1} \Delta C /(C^2+1)$ , plotted as a function of $C/(C + 1)$ , for an oblate spheroid; $C/(C + 1)=0$ and $C/(C+ 1)=1$ correspond to the spinning and tumbling modes, respectively. Panels (a) and (c), with $\Delta C\lt 0$ and $\Delta C\gt 0$ $\forall \,C$ , respectively, correspond to the stable-spinning and stable-tumbling regimes. The corresponding unit-sphere trajectories are shown in (b) and (d) with ${\textit{De}}=0.1$ and $0.3$ ; red arrows indicate the direction of drift.

Figure 4. Long description

The image consists of four panels. Panels (a) and (c) display line graphs plotting the normalized change in the orbit constant as a function of C/(C+1) for an oblate spheroid. Panel (a) shows a blue curve for the stable-spinning regime with κ = 0.3 and ε = −0.1, while panel (c) shows a red curve for the stable-tumbling regime with κ = 0.3 and ε = −0.7. Panels (b) and (d) depict 3D trajectories on a unit sphere for the corresponding regimes shown in panels (a) and (c), respectively. The trajectories are illustrated with black contour lines, and red arrows indicate the direction of drift. The x, y, and z axes are labeled, with C = 0 and C = ∞ marked on the sphere.

Figure 5.

The drift due to elasticity as characterized by the normalized change in the orbit constant in a single Jeffery period, ${\textit{De}}^{-1} \Delta C/(C^2+1)$ , plotted as a function of $C/(C + 1)$ , for a prolate spheroid; a larger value, ${\textit{De}}=5.0$ , is chosen for clarity. (a) Here $\Delta C$ has a negative-slope zero crossing, corresponding to the stable kayaking regime. (c) Here $\Delta C$ has a positive-slope zero crossing, corresponding to the unstable kayaking regime. Unit-sphere trajectories are shown in (b) and (d), where the red curve is the intermediate kayaking orbit and blue arrows indicate the direction of drift. In (a) and (c), the blue dot marks the zero crossing of $\Delta C$ .

Figure 5. Long description

The image consists of four subfigures labeled (a), (b), (c), and (d). Subfigure (a) presents a graph plotting the normalized change in the orbit constant in a single Jeffery period as a function of C/(C+1) for a prolate spheroid. The graph has a negative-slope zero crossing, corresponding to the stable kayaking regime. The green curve represents the data, and the blue dot marks the zero crossing. Subfigure (b) shows a unit-sphere trajectory with a red curve indicating the intermediate kayaking orbit and blue arrows showing the direction of drift. Subfigure (c) presents another graph plotting the normalized change in the orbit constant in a single Jeffery period as a function of C/(C+1) for a prolate spheroid. This graph has a positive-slope zero crossing, corresponding to the unstable kayaking regime. The yellow curve represents the data, and the blue dot marks the zero crossing. Subfigure (d) shows a unit-sphere trajectory with a red curve indicating the intermediate kayaking orbit and blue arrows showing the direction of drift. The trajectories in subfigures (b) and (d) illustrate the different orientations and drift directions of the prolate spheroid under varying conditions.

Plots of ${\textit{De}}^{-1} \Delta C/(C^2+1)$ versus $C/(C+1)$ may now be used to characterize the nature of spiralling trajectories in presence of weak viscoelasticity. As $C/(C+1)$ varies from $0$ to $1$ , one moves from the spinning to the tumbling mode; further, $\Delta C$ is normalized by $C^2+1$ , so that $\Delta C/(C^2+1) = 0$ for $C = 0$ and for $C \rightarrow \infty$ – as pointed out in § 1, in the context of the theoretical effort of Brunn (Reference Brunn1977b ), symmetry dictates that spinning and tumbling modes are always allowed equilibria, and it is only their stability that depends on the relevant physical parameters. Our detailed calculations for different $\xi _0$ (or $\kappa$ ) and $\epsilon$ reveal four different types of ${\textit{De}}^{-1} \Delta C/(C^2+1)$ versus $C/(C+1)$ plots, examples of which are illustrated in figures 4 and 5. In figure 4(a) corresponding to $(\kappa ,\epsilon ) \equiv (0.3,-0.1)$ , $\Delta C$ is seen to be negative for all finite non-zero $C$ , implying that the spheroid drifts towards the spinning mode starting from any initial orientation; the associated unit-sphere trajectory in figure 4(b) spirals in towards the vorticity axis. In contrast, in figure 4(c), corresponding to $(\kappa ,\epsilon ) \equiv (0.3,-0.7)$ , $\Delta C$ is positive for all $C$ , implying that a spheroid for these parameters drifts towards the tumbling mode starting from any initial orientation; the associated unit-sphere trajectory in figure 4(d) spirals out towards the flow-gradient plane. The $\Delta C$ versus $C$ curves in figures 5(a) and 5(c) exhibit a zero crossing, corresponding to a limit cycle equilibrium intermediate between the spinning and tumbling modes; for small ${\textit{De}}$ , the limit cycle closely approximates a finite- $C$ Jeffery orbit. When the slope at the zero-crossing is negative, as in figure 5(a), the limit cycle is stable, and almost all initial orientations (except those corresponding to the tumbling and spinning modes) approach a finite- $C$ kayaking motion for sufficiently long times – figure 5(b) shows a pair of spiralling trajectories approaching the limit cycle (red curve) for long times. When the slope at the zero-crossing is positive, as in figure 5(c), the limit cycle is an unstable equilibrium, and acts to divide each unit hemisphere into distinct basins of attraction for the tumbling and spinning modes. All initial orientations in one basin approach the spinning mode, while the ones in the other approach the tumbling mode. Figure 5(d) depicts a pair of trajectories starting from initial orientations on either side of the unstable limit cycle (red curve), and spiralling away from it. For future reference, the scenarios in figures 4(a,c), 5(a) and 5(c) will be referred to as Regimes $1$ , $2$ , $3$ and $4$ , respectively.

3.3. Viscoelastic-orientation dynamics on the $\kappa$ – $\epsilon$ plane

In this section, we will show ${\textit{De}}^{-1} \Delta C/(C^2+1)$ versus $C/(C+1)$ plots over a wide range of spheroid aspect ratios, in order to illustrate transitions in the orientation dynamics that occur with changing $\kappa$ , at a fixed $\epsilon$ . This leads naturally to an organization of the different orientation-dynamics regimes, identified in figures 4 and 5, on the $\kappa -\epsilon$ plane. To do the above, there is the need to scale ${\textit{De}}^{-1} \Delta C/(C^2+1)$ appropriately, so it remains finite both in the near-sphere limit, and for extreme aspect ratios (slender fibres and flat disks). Barring exceptional values of $\epsilon$ , ${\textit{De}}^{-1} \Delta C$ scales $O(\xi _{0}^{-2})$ in the near-sphere limit, with the sign of this contribution depending on whether one approaches the sphere from the prolate or the oblate side. Here ${\textit{De}}^{-1} \Delta C$ is $O(1/\sqrt {\xi _0-1})$ in the slender-fibre and flat-disk limits. The divergence for $\xi _0 \rightarrow 1$ arises because the near-meridional nature of the Jeffery orbits, in the said limits, leads to them being squeezed together in an increasingly small neighbourhood of the flow or gradient axis. As a result, a fixed $O(De)$ angular displacement in these neighbourhoods amounts to an increasingly large drift when measured in units of $C$ . The divergent scaling above already points to a breakdown of the weak-drift approximation when $\Delta C \sim De/\sqrt {\xi _0-1} \sim O(1)$ ; as already mentioned, this breakdown and the associated rotation arrest is examined in § 4. It is worth adding here that, for the inertial case, $\Delta C \sim Re.O(1/\sqrt {\xi _0-1}\ln (\xi _0-1))$ in the slender-fibre limit (Subramanian & Koch Reference Subramanian and Koch2005; Dabade et al. Reference Dabade, Marath and Subramanian2016), with the logarithmic smallness arising due to different asymptotic regions that contribute dominantly to fibre rotation in the inertial (outer region with $r \sim O(1)$ in units of the fibre length) and viscoelastic cases (matching region with $\sqrt {\xi _0-1} \ll r \ll 1$ ) – these aspects are clarified by the analysis in Appendix B. With the above scalings in mind, we replace ${\textit{De}}^{-1} \Delta C /(C^2+1)$ by ${\textit{De}}^{-1} \Delta C \xi _0^{ {3}/{2}}\sqrt {\xi _0-1}/(C^2+1)$ in the orbital drift plots below.

While the generic $\xi _0$ -scalings of $\Delta C$ have been clarified above, the exceptional scenarios where these scalings do not hold are also important, since they correspond to boundaries demarcating the different orientation dynamics regimes. To examine these instances, we consider the different limiting forms of $\Delta C$ , as defined by (3.8). In the near-sphere limit, one finds $\lim _{\xi _0\rightarrow \infty } De^{-1} \Delta C = \pm \pi C (8 \epsilon +3)/(14 \xi _{0}^{2})$ , with the plus (minus) sign corresponding to a prolate (oblate) spheroid. This identifies a critical value of $\epsilon \,(=-3/8)$ at which ${\textit{De}}^{-1} \Delta C \sim O(\xi _0^{-4})$ in the near-sphere limit, and that marks a bifurcation in the orientation dynamics for $\kappa \rightarrow 1^\pm$ – this bifurcation mediates a transition from Regime $2$ to $1$ ( $1$ to $2$ ) as $\epsilon$ decreases below $-3/8$ on the prolate (oblate) side. In the slender-fibre limit, one finds $\lim _{\xi _0\rightarrow 1} De^{-1} \Delta C = {\pi C (2 \epsilon +1)}/{2 \sqrt {2} \sqrt {\xi _{0}-1}}$ , implying that the point $(\kappa ,\epsilon ) \equiv (\infty ,-1/2)$ now marks a bifurcation in the orientation dynamics, with $\lim _{\xi _0\rightarrow 1} De^{-1} \Delta C \sim O(1/\sqrt {\xi _0-1}\ln (\xi _0-1))$ at $\epsilon =-1/2$ . Inferring the change in the orientation dynamics regime in the slender-fibre limit is more subtle, however, since the limiting expression for $\Delta C$ above is not uniformly valid. A separate asymptote exists for small $C$ , and is given by ${\textit{De}}^{-1} \Delta C = -\sqrt {2}\pi C (1+\epsilon ) \sqrt {\xi _0-1}$ . For $\epsilon \gt -1/2$ , this small- $C$ asymptote has a sign (negative) opposite to the asymptote for $C \sim O(1)$ above, implying a zero crossing with a positive slope in between, with a positive slope, and therefore, an unstable intermediate limit cycle. For $-1 \lt \epsilon \lt -1/2$ , both asymptotes have the same (negative) sign, and the spinning mode alone is stable. Finally, for $\epsilon \lt -1$ , the small- $C$ asymptote reverses sign, and one must then have a stable intermediate limit cycle. Thus, in the limit $\kappa \rightarrow \infty$ , one must transition from Regime $4$ , through Regime  $1$ , to Regime $3$ , with decreasing $\epsilon$ . In the flat-disk limit, one finds $\lim _{\xi _0\rightarrow 1} De^{-1} \Delta C \sim -( {\pi [2 (\sqrt {C^2+1}-1 ) (6 \epsilon +1)+C^2 ]})/({6 \sqrt {2} C \sqrt {\xi _0-1}})$ which is uniformly valid for all $C$ . For $C \rightarrow 0$ , this asymptote reduces to $-( {\pi }/{\sqrt {2}\sqrt {\xi _0-1}})[C(\epsilon +1/3) -({C^3}/{4})(\epsilon + 1/6)]+ O(C^5)$ , which shows that the point $(\kappa ,\epsilon ) \equiv (0,-1/3)$ corresponds to a bifurcation in the orientation dynamics. For $\epsilon \gt -1/3$ , the spinning mode alone is stable, while for $\epsilon$ less than but close to $-1/3$ , one obtains $\Delta C = 0$ for $C^* \approx 2[-6(\epsilon +1/3)]^{ {1}/{2}}$ from the small- $C$ asymptote, which corresponds to a stable limit cycle emerging from the vorticity axis; one may, in fact, obtain the location of the stable limit cycle for arbitrary $\epsilon$ , $C^* = 2\sqrt {2}[(3 \epsilon +1) (6 \epsilon +1)]^{{1}/{2}}$ , using the aforementioned flat-disk asymptote. Thus, in the limit $\kappa \rightarrow 0$ , one must transition from Regime $1$ to Regime $3$ as $\epsilon$ decreases across $-1/3$ . Consideration of the next term in the flat-disk asymptote for $\Delta C$ shows that, for any small but finite $\kappa$ , there must be a second transition from Regime $3$ to Regime $2$ across a large but negative $\epsilon$ of $O(1/\sqrt {\xi _0-1})$ .

The above analysis helps identify the limiting forms of the bifurcation loci, that separate the different orientation-dynamics regimes on the $\kappa -\epsilon$ plane. A pair of loci on the prolate side must asymptote to the vertical lines $\epsilon = -1$ and $\epsilon = -1/2$ in the limit $\kappa \rightarrow \infty$ . On the oblate side, one of the loci must asymptote to the line $\epsilon = -1/3$ , with the other assuming the limiting form $\epsilon \propto -\kappa ^{-1}$ , in the limit $\kappa \rightarrow 0$ . Further, a critical point $(\kappa ,\epsilon ) \equiv (1,-3/8)$ on the sphere-line is expected to act as an anchor point for one or more of the aforementioned curves, and control the bifurcations in the near-sphere limit. An additional benefit of the arguments above is a broad identification of $\epsilon$ -intervals corresponding to distinct bifurcation sequences for the spheroid orientation dynamics. These intervals are $-\infty \lt \epsilon \lt -1$ , $-1 \lt \epsilon \lt -1/2$ and $-1/2 \lt \epsilon \lt \infty$ on the prolate side. Neglecting the locus that diverges to infinity for $\kappa \rightarrow 0$ , the intervals are $-\infty \lt \epsilon \lt -1/3$ and $-1/3 \lt \epsilon \lt \infty$ on the oblate side. Figures 6 and 7 below show the ${\textit{De}}^{-1} \Delta C\xi _0^{ {3}/{2}}\sqrt {\xi _0-1}/(C^2+1)$ versus $C/(C+1)$ plots for $\epsilon$ values in each of the aforesaid intervals.

Figure 6.

The scaled orbital drift, ${\textit{De}}^{-1} \Delta C \xi _{0}^{3/2}\sqrt {\xi _{0}-1}/(C^2+1)$ , plotted as a function of $C/(C + 1)$ , for prolate spheroids of a varying aspect-ratio $\kappa$ with $\epsilon$ fixed; the black dashed and dash–dotted curves denote the near-sphere and slender-fibre asymptotes, respectively. Here (a–b) $\epsilon =-2$ , (c–d) $\epsilon =-0.6$ and (e–f) $\epsilon =-0.1$ . The inset views in (b) and (f) highlight the emergence of a zero-crossing from the spinning end ( $C=0$ ) above a threshold $\kappa$ ; the zero-crossing corresponds to a stable and unstable limit cycle in the respective cases. Black arrows, either straight or curved, indicate the movement of the drift curves with increasing $\kappa$ .

Figure 6. Long description

The image contains six line graphs labeled (a) through (f), each showing the scaled orbital drift as a function of C/(C+1) for prolate spheroids with varying aspect ratios. The black dashed and dash-dotted curves denote the near-sphere and slender-fibre asymptotes, respectively. Each graph represents different values of the aspect ratio, with curves color-coded and labeled accordingly. The inset views in graphs (b) and (f) highlight the emergence of a zero-crossing from the spinning end above a threshold, corresponding to stable and unstable limit cycles. Black arrows indicate the movement of the drift curves with increasing aspect ratio. The graphs illustrate how the drift curves evolve with changes in the aspect ratio, providing insights into the orientation dynamics of prolate spheroids in viscoelastic fluids.

Figure 6(a–f) illustrate the orbital drift curves for a prolate spheroid, over the entire range of aspect ratios, for $\epsilon = -2$ (figures 6 a and 6 b), $-0.6$ (figures 6 c and 6 d) and $-0.1$ (figures 6 e and 6 f); note that $\epsilon = -0.6$ falls in the range typical of polymer solutions. The black dashed and dash–dotted curves in each of these figures denote the limiting near-sphere and slender-fibre asymptotes, respectively. The drift curves in figure 6(a) conform to the stable spinning regime. The scaled orbital drift decreases in magnitude starting from the near-sphere asymptote, attaining its minimal limit for $\kappa \approx 2$ . The movement of the drift curve reverses thereafter. As shown in figure 6(b), the drift increases in magnitude for $\kappa$ larger than $2$ , overshooting the slender-fibre asymptote for $\kappa \gtrsim 15$ , before turning around and converging back to this asymptote in the limit $\kappa \rightarrow \infty$ . The convergence is only logarithmic, with the difference between the asymptote and the large- $\kappa$ curves decreasing as $O(\ln \kappa )^{-1}$ , as also evident from the perceptible difference between the two for $\kappa = 10^4$ . The inset in figure 6(b) shows a magnified view of the region near $C = 0$ , highlighting the emergence of a negative-slope zero-crossing from $C = 0$ at $\kappa \approx 9.324$ . The zero-crossing moves towards larger $C$ with increasing $\kappa$ , but turns around thereafter, heading back towards $C = 0$ , and eventually coalescing with it in the limit $\kappa = \infty$ . Thus, for $\epsilon = -2$ , a prolate spheroid transitions from Regime $1$ to Regime $3$ with increasing $\kappa$ . This transition is marked by a stable limit cycle emerging from the vorticity axis at $\kappa \approx 9.324$ , attaining a maximum size at $\kappa \approx 1000$ , and collapsing back onto the vorticity axis for $\kappa \rightarrow \infty$ .

Figures 6(c) and 6(d) depict the orbital drift for $\epsilon =-0.6$ . In this case, the prolate spheroid drifts to the spinning mode regardless of $\kappa$ , with there being no bifurcations. Interestingly, the movement of the scaled drift-curve from the near-sphere to the slender-fibre asymptote continues to have a non-monotonic character. With increasing $\kappa$ , the drift first decreases away from the near-sphere asymptote until $\kappa \approx 2$ , then increases and overshoots the near-sphere asymptote for $\kappa \approx 10$ , before turning around and converging at a logarithmic rate towards the slender-fibre asymptote. Figures 6(e) and 6(f) illustrate the drift of a prolate spheroid for $\epsilon =-0.1$ . The drift curves are now positive for all $C$ to begin with. Mirroring the scenario for $\epsilon = -2$ , the drift curve in figure 6(e) decreases relative to the near-sphere asymptote, attaining a minimum for $\kappa \approx 2$ . It increases monotonically thereafter, gradually approaching the slender-fibre asymptote. The inset view in figure 6(f) shows the emergence of a positive-slope zero-crossing (an unstable limit cycle) from $C=0$ , at $\kappa \approx 3.695$ , that eventually collapses back into $C = 0$ for $\kappa \rightarrow \infty$ . Thus, a prolate spheroid, for $\epsilon = -0.1$ , transitions from Regime $2$ to Regime $4$ with increasing $\kappa$ .

Figure 7(a–d) illustrate the orbital drift curves for an oblate spheroid, over the entire range of aspect ratios, for $\epsilon = -0.6$ (figures 7 a and 7 b) and $-0.1$ (figures 7 c and 7 d). In figure 7(a), the drift curves conform to the stable tumbling regime to begin with. The portion of the drift curve for $C/(C+1) \lesssim 0.6$ again exhibits a non-monotonic movement, decreasing from the near-sphere asymptote to begin with, but heading back up towards the flat-disk asymptote for $\kappa \lesssim 0.8$ . While the movement of the drift curve from the near-sphere to the flat-disk asymptote has a monotonic character for larger $C/(C+1)$ , there is nevertheless a bifurcation at $\kappa \approx 0.28$ due to the emergence of a negative-slope zero-crossing from the tumbling end ( $C= \infty$ ), and that asymptotes to a finite $C\,(=C^*)$ in the flat-disk limit. The expression for $C^*$ , for $\kappa \rightarrow 0$ , was obtained above, and equals $4$ for the chosen $\epsilon$ . Thus, an oblate spheroid, for $\epsilon = -0.6$ , transitions from Regime $2$ to Regime  $3$ with decreasing $\kappa$ . Figures 7(c) and 7(d) depict the drift curves for $\epsilon =-0.1$ . While all of the drift curves conform to the stable spinning regime in this case, the movement with decreasing $\kappa$ has a non-monotonic character, with the drift attaining a maximum at $\kappa \approx 0.25$ .

Figure 7.

The scaled orbital drift, ${\textit{De}}^{-1} \Delta C \xi _{0}^{3/2}\sqrt {\xi _{0}-1}/(C^2+1)$ , plotted as a function of $C/(C + 1)$ , for oblate spheroids of a varying aspect-ratio $\kappa$ with $\epsilon$ fixed. The black dashed and dash–dotted curves denote the near-sphere and flat-disk asymptotes, respectively; the latter exhibits a zero-crossing at $C^* = 2\sqrt {2}[(3 \epsilon +1) (6 \epsilon +1)]^{({1}/{2})}$ for $\epsilon \lt -1/3$ . Here (a,b) $\epsilon =-0.6$ and (c,d) $\epsilon =-0.1$ . Black arrows, either straight or curved, indicate the movement of the drift curves with decreasing $\kappa$ .

Figure 7. Long description

The image contains four graphs labeled (a), (b), (c), and (d), each showing the scaled orbital drift, De^{−1} \Delta C ξ0^{3/2} (ξ0−1)^{1/2} / (1+C^2), plotted as a function of C/(C+1) for oblate spheroids with varying aspect ratios κ. The black dashed and dash-dotted curves denote the near-sphere and flat-disk asymptotes, respectively. The latter exhibits a zero-crossing at C/(C+1) = 0.8 for κ = 0.1. The graphs illustrate how the drift curves move with decreasing κ, indicated by black arrows. Each graph shows different ranges of κ values, with (a) and (b) focusing on \epsilon=-0.6  and (c) and (d) on \epsilon=-0.1. The curves demonstrate the relationship between the aspect ratio of the spheroids and their orbital drift, highlighting the differences in behavior as the aspect ratio changes.

Having characterized the orientation dynamics bifurcations as a function of varying $\kappa$ , for specific $\epsilon$ values, figure 8 shows the different small- ${\textit{De}}$ orientation dynamics regimes on the $\kappa -\epsilon$ plane, for an arbitrary aspect ratio spheroid, in an ambient simple shear flow. It is worth reiterating that ${\textit{De}}$ here is not fixed; instead, in order for the orbital-drift interpretation to remain valid, ${\textit{De}}$ is constrained to be increasingly small for $\kappa \rightarrow 0$ ( ${\textit{De}} \ll \kappa$ ) and $\infty$ ( ${\textit{De}} \ll \kappa ^{-1}$ ). The boundaries, that demarcate the different regimes in figure 8, have a more complicated form than that inferred using arguments based on the limiting forms of $\Delta C$ above. One of the prolate-side boundaries starts from its large- $\kappa$ asymptotic form, $\epsilon = -1$ , but turns around for finite $\kappa$ , eventually heading to negative infinity parallel to the $\epsilon$ -axis, with $\kappa \rightarrow 6.963$ in his limit. The second prolate-side boundary continues down to the sphere line ( $\kappa =1$ ), intersecting it at $\epsilon = -3/8$ . It then continues further into the oblate side, eventually asymptoting to the curve $\epsilon = (2 \pi \kappa ^{-1})/ (-1168-24 \pi +117 \pi ^2 )$ in the limit $\kappa \rightarrow 0$ . The oblate-side boundary that starts off along $\epsilon = -1/3$ in the limit $\kappa \rightarrow 0$ , continues up and intersects the sphere-line at $\epsilon = -3/8$ . It then turns around on the prolate side, heading off to positive infinity parallel to the $\epsilon$ -axis, with $\kappa$ again approaching $6.963$ in this limit. The full forms of the equations governing the two bounding curves in figure 8, in the form $\epsilon \equiv f(\xi _0)$ , are given in Appendix A.4. The two bounding curves, that intersect the sphere line at $\epsilon = -3/8$ , also intersect in a second point, $(\epsilon ,\kappa ) \equiv (-0.288,2.485)$ , on the prolate side. The unit-sphere exhibits a closed Jeffery-orbit topology at the latter point, even for small but finite ${\textit{De}}$ ; the orbits are (trivially) circles at $(-3/8,1)$ . The portions of the boundaries between $(-3/8,1)$ and $(-0.288,2.485)$ enclose an island where the orientation dynamics conforms to Regime $3$ . Accounting for this island, the transitions in the small- ${\textit{De}}$ orientation dynamics, for $\kappa$ increasing from zero to infinity, may be organized across the following six distinct $\epsilon$ -intervals:

1. (i) $-\infty \lt \epsilon \lt -1$ , stable kayaking $\rightarrow$ stable tumbling $\rightarrow$ stable spinning $\rightarrow$ stable kayaking;

2. (ii) $-1\lt \epsilon \lt -1/2$ , stable kayaking $\rightarrow$ stable tumbling $\rightarrow$ stable spinning;

3. (iii) $-1/2\lt \epsilon \lt -3/8$ , stable kayaking $\rightarrow$ stable tumbling $\rightarrow$ stable spinning $\rightarrow$ unstable kayaking sandwiched between stable spinning and tumbling;

4. (iv) $-3/8 \lt \epsilon \lt -1/3$ , stable kayaking $\rightarrow$ stable spinning $\rightarrow$ stable tumbling $\rightarrow$ stable kayaking $\rightarrow$ stable spinning $\rightarrow$ unstable kayaking sandwiched between stable spinning and tumbling;

5. (v) $-1/3 \lt \epsilon \lt -0.288$ , stable spinning $\rightarrow$ stable tumbling $\rightarrow$ stable kayaking $\rightarrow$ stable spinning $\rightarrow$ unstable limit cycle;

6. (vi) $-0.288 \lt \epsilon \lt \infty$ , stable spinning $\rightarrow$ stable tumbling $\rightarrow$ unstable kayaking sandwiched between stable spinning and tumbling.

Figure 8.

The orientation dynamics regimes on the $\kappa$ – $\epsilon$ plane, for small but finite ${\textit{De}}$ ; the four different regimes are numbered as follows: $1$ , stable spinning; $2$ , stable tumbling; $3$ , stable kayaking; $4$ ,unstable kayaking sandwiched between stable spinning and tumbling modes. The horizontal black line is $\kappa =1$ , with regions above and below corresponding to prolate and oblate spheroids, respectively; the horizontal dashed line, $\kappa = 6.963$ , is the infinite- $\epsilon$ asymptote. The cyan vertical band denotes the range $\epsilon \in (-0.7,-0.5)$ relevant to polymer solutions. Vertical dashed lines serve to demarcate the $\epsilon$ -intervals corresponding to distinct bifurcation sequences for the orientation dynamics, as a function of $\kappa$ .

Figure 8. Long description

A line graph illustrates the orientation dynamics regimes on the plane for small but finite Deborah number. The graph is divided into four distinct regimes, each numbered and color-coded: stable spinning, stable tumbling, stable kayaking, and unstable kayaking. The horizontal black line represents a unit aspect ratio, with regions above and below corresponding to prolate and oblate spheroids, respectively. The horizontal dashed line marks the infinite-\epsilon asymptote. A cyan vertical band denotes the range relevant to polymer solutions. Vertical dashed lines demarcate the epsilon-intervals corresponding to distinct bifurcation sequences for the orientation dynamics as a function of epsilon. The graph includes various labeled intervals and points, such as negative infinity to 6.963, negative 1 to infinity, negative 0.288 to 2.485, and others, indicating different dynamic behaviors. All values are approximate.

Figure 9.

The scaled orbital drift, ${\textit{De}}^{-1} \Delta C \xi _{0}^{3/2}\sqrt {\xi _{0}-1}/(C^2+1)$ , plotted as a function of $C/(C + 1)$ , for spheroids with $\kappa$ increasing from zero to infinity; $\epsilon =-0.36$ . Panels (a) and (b) pertain to oblate spheroids, and (d)–(f) to prolate spheroids; (c) shows the transition of the drift curve from the oblate to the prolate side. Black dashed and dash–dotted curves denote the near-sphere and flat-disk asymptotes in (a) and (b), and the near-sphere and slender-fibre asymptotes in (c) and (f). Black arrows in all figures denote the direction of movement of the drift curve with increasing $\kappa$ .

Figure 9. Long description

The image contains six panels (a) through (f) of line graphs depicting the scaled orbital drift, De^{−1} \Delta C  ξ0^{3/2} (ξ0−1)^{1/2} / (1 + C^2), as a function of C/(C + 1) for spheroids with increasing aspect ratios (κ). Panels (a) and (b) pertain to oblate spheroids, panels (d), (e), and (f) to prolate spheroids, and panel (c) shows the transition of the drift curve from the oblate to the prolate side. Each panel contains multiple lines representing different values of κ, ranging from 0.001 to 10000. Black dashed and dash-dotted curves denote the near-sphere and flat-disk asymptotes in (a) and (b), and the near-sphere and slender-fibre asymptotes in (c) and (f). Black arrows in all figures denote the direction of movement of the drift curve with increasing κ. The graphs illustrate how the drift behavior changes with the aspect ratio of the spheroids, showing different trends and intersections for oblate and prolate shapes. All values are approximated.

The Regime $3$ island in figure 8 implies that a vertical path corresponding to any $\epsilon$ in the interval $(-3/8,-1/3)$ exhibits a more complex sequence of bifurcations than those identified earlier in figures 6 and 7, based on the limiting forms of $\Delta C$ ; this sequence corresponds to item (iv) in the list above. The scaled orbital drift curves for $\epsilon = -0.36\,\in \,(-3/8,-1/3)$ are shown in figure 9(a–f). One starts off in figure 9(a) with the drift curve exhibiting a negative-slope zero crossing at $C^* =0.862$ for $\kappa = 0$ . The zero-crossing migrates towards the spinning end, coalescing with it at $\kappa \approx 0.265$ . The drift in figure 9(b) decreases in magnitude towards the near-sphere asymptote, while conforming to a stable spinning regime. Figure 9(c) shows the expected transition to a stable tumbling regime as $\kappa$ increases across unity. A negative-slope zero-crossing emerges from the tumbling end at $\kappa \approx 1.098$ , and as shown in figure 9(d), moves towards the spinning end with increasing $\kappa$ , coalescing with it at $\kappa \approx 1.499$ ; this movement also leads to $\Delta C$ reversing its sign (from positive to negative). Figure 9(e) shows that, with a further increase in $\kappa$ , the drift magnitude exhibits a non-monotonic variation, attaining a maximum at $\kappa \approx 4$ . In figure 9(f), a positive-slope zero crossing is seen to emerge from the spinning end at $\kappa \approx 8.154$ , again coalescing with it in the limit $\kappa \rightarrow \infty$ . With increasing $\kappa$ , therefore, figure 9(a–f) describe the following sequence of regimes: (i) $0\lt \kappa \lt 0.265$ , stable kayaking; (ii) $0.265\lt \kappa \lt 1$ , stable spinning; (iii) $1\lt \kappa \lt 1.098$ , stable tumbling; (iv) $1.098\lt \kappa \lt 1.499$ , stable kayaking; (v) $1.499\lt \kappa \lt 8.154$ , stable spinning; (vi) $8.154\lt \kappa \lt \infty$ , unstable kayaking.