2.1. Using a generalized reciprocal theorem to obtain ${\boldsymbol {G}^{PIST}}$ in terms of $\boldsymbol {\varPi }$

In its original form, $\boldsymbol {G}^{PIST}=\int _{\boldsymbol {r}_{p}} {\rm d}S\,\boldsymbol {r}\times \boldsymbol {\tau }^{P}\boldsymbol{\cdot} \boldsymbol {n}$ is a function of $\boldsymbol {\tau }^{P}$ which in turn is driven by $\boldsymbol {\varPi }$ through (2.8a,b). In this section, using the balance of the divergence of the polymeric stress ($\boldsymbol {\varPi }$) and $\boldsymbol {\tau }^{P}$ from (2.8a,b) and using a generalized reciprocal theorem we derive an expression to obtain $\boldsymbol {G}^{PIST}$ directly from $\boldsymbol {\varPi }$ without the need to compute $\boldsymbol {\tau }^{P}$.

The auxiliary or comparison problem in a generalized reciprocal theorem must be chosen based on the surface integral one is interested in evaluating. The effect of the torque is to rotate the particle. Hence, we consider the Stokes flow around a particle rotating in a quiescent Newtonian fluid. We define the following auxiliary Stokes problem for a ‘velocity’ field consisting of a rank-two pseudotensor $\boldsymbol {b}$, a ‘pressure’ field that is a pseudovector, $\boldsymbol {q}$, and a ‘fluid stress’ field that is a rank-three pseudotensor $\boldsymbol {B}$,

(2.18a–c)\begin{equation} \frac{\partial B_{ijk}}{\partial r_i}=0,\quad \frac{\partial b_{ij}}{\partial r_i}=0,\quad B_{ijk}=-\delta_{ij}q_k+\frac{\partial b_{jk}}{\partial r_i}+\frac{\partial b_{ik}}{\partial r_j},\end{equation}

with boundary condition

(2.19)\begin{equation} b_{ij}=\epsilon_{ijk}r_k, \quad \text{on particle surface}, \quad \text{and}, \quad b_{ij}=0, \quad \text{as } |\boldsymbol{r}|\rightarrow\infty, \end{equation}

where $\epsilon _{ijk}$ is the permutation tensor. Here $\boldsymbol {b}\boldsymbol{\cdot}\boldsymbol {\omega }_{auxiliary}$ is the velocity field around a particle rotating with an angular velocity $\boldsymbol {\omega }_{auxiliary}$ in a quiescent Newtonian fluid. Using the definitions of $\boldsymbol {\tau }^{P}$ and $\boldsymbol {B}$ in terms of the respective velocities, $\boldsymbol {u}^{P}$ and $\boldsymbol {b}$, incompressibility of the velocities in (2.8a,b) and (2.18a–c), $\boldsymbol {\nabla }\boldsymbol{\cdot}\boldsymbol {B}=0$ from (2.18a–c) and the symmetry of $B_{lki}$ and ${\tau }^{P}_{lk}$ about the $l$ and $k$ indices we obtain

(2.20)\begin{equation} {{\tau}^{P}_{lk}\frac{\partial b_{ki}}{\partial r_l}=\frac{\partial B_{lki} {u}^{P}_k}{\partial r_l}.}\end{equation}

Using the balance of the divergence of $\boldsymbol {\varPi }$ and $\boldsymbol {\tau }^{P}$ from (2.8a,b), the volume integral of (2.20) in a region bounded by the particle surface, $\boldsymbol {r}_{p}$, and a far-field spherical surface at $|\boldsymbol {r}|\rightarrow \infty$ is

(2.21)\begin{equation} {\int_{Fluid} {\rm d}V\,\frac{\partial}{\partial r_l}[{\tau}^{P}_{lk}b_{ki}-B_{lki}{u}^{P}_k]=- \int_{Fluid} {\rm d}V\,b_{ki}\frac{\partial{\varPi}_{lk}}{\partial r_l}.} \end{equation}

Using the divergence theorem, the left-hand side of the above equation can be written as

(2.22)\begin{align} \int_{Fluid} {\rm d}V\,\frac{\partial}{\partial r_l}[{\tau}^{P}_{lk}b_{ki}-B_{lki}{u}^{P}_k]& = \int_{|\boldsymbol{r}|\rightarrow\infty} {\rm d}S\,n_l[{\tau}^{P}_{lk}b_{ki}-B_{lki}{u}^{P}_k]\nonumber\\ &\quad -\int_{\boldsymbol{r}_{p}}{\rm d}S\,n_l[{\tau}^{P}_{lk}b_{ki}- B_{lki}{u}^{P}_k], \end{align}

where the surface normal $n_l$ points into the fluid (away from particle region) on the particle surface. The fluid velocity due to a particle that exerts a force (force dipole) on the fluid decays as $1/r$ ($1/{r}^2$) in the far field. Hence, the velocities $\boldsymbol {u}^{P}$ and $\boldsymbol {b}$ scale as $1/r$ and $1/{r}^2$, respectively, and the stresses $\boldsymbol {\tau }^{P}$ and $\boldsymbol {B}$ scale as $1/{r}^2$ and $1/{r}^3$. Therefore, the first surface integral in (2.22) vanishes. Using the boundary conditions for $\boldsymbol {u}^{P}$ and $\boldsymbol {b}$ from (2.9) and (2.19) in the second surface integral in (2.22) we obtain

(2.23)\begin{equation} \int_{Fluid} {\rm d}V\,\frac{\partial}{\partial r_l}({\tau}^{P}_{lk}b_{ki}-B_{lki}{u}^{P}_k)= \int_{\boldsymbol{r}_{p}}{\rm d}S\,n_l{\tau}^{P}_{lk}\epsilon_{kim}r_m=-{G}^{PIST}_i.\end{equation}

Therefore, from (2.21) and (2.23), for a particle of any shape,

(2.24)\begin{equation} \boldsymbol{G}^{PIST}=\int_{Fluid} {\rm d}V(\boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{\varPi}}) \boldsymbol{\cdot}\boldsymbol{b}.\end{equation}

This completes the derivation expressing $\boldsymbol {G}^{PIST}$ in terms of the polymer stress $\boldsymbol {\varPi }$ and a quasisteady two-tensor field $\boldsymbol {b}$ that is dependent on the particle shape and is a solution of the auxiliary Stokes problem defined by (2.18a–c) and (2.19).

In a linear imposed flow such as a simple shear, the undisturbed polymer stress, $\boldsymbol {\varPi }^U$, i.e. the polymer stress without the particle, is spatially constant. Hence,

(2.25)\begin{equation} \boldsymbol{G}^{PIST}=\int_{Fluid} {\rm d}V(\boldsymbol{\nabla}\boldsymbol{\cdot}( \boldsymbol{\varPi}-\boldsymbol{\varPi}^U))\boldsymbol{\cdot}\boldsymbol{b}.\end{equation}

Using the chain rule and divergence theorem,

(2.26)\begin{align} {G}^{PIST}_i&=\int_{|\boldsymbol{r}|\rightarrow\infty} {\rm d}S({\varPi}_{lk}-{\varPi}^{U}_{lk}) b_{ki} n_l-\int_{\boldsymbol{r}_{p}} {\rm d}S({\varPi}_{lk}-{\varPi}^{U}_{lk}) b_{ki}n_l\nonumber\\ &\quad -\int_{Fluid} {\rm d}V({\varPi}_{lk}-{\varPi}^{U}_{lk})\frac{\partial b_{ki}}{\partial r_l}. \end{align}

The disturbance of the polymer stress created by the particle, ${\varPi }_{lk}-{\varPi }^{U}_{lk}$, also decays as $1/{r}^2$ in the far-field since it is forced by the disturbance to the far-field velocity gradients. In the case of Oldroyd-B fluids or other dumbbell models such as FENE-P and Giesekus (Bird et al. Reference Bird, Armstrong and Hassager1987), this far-field scaling of ${\varPi }_{lk}-{\varPi }^{U}_{lk}$ is ascertained by linearizing the polymer constitutive equation (for example (2.4) and (2.5) for the Oldroyd-B model) about the far-field velocity and polymer conformation to obtain a governing equation for ${\varPi }_{lk}-{\varPi }^{U}_{lk}$. Therefore, the first surface integral in the above equation vanishes, and using the surface boundary condition of (2.19) leads to

(2.27)\begin{equation} {G}^{PIST}_i=-\int_{\boldsymbol{r}_{p}} {\rm d}S({\varPi}_{lk}-{\varPi}^{U}_{lk}) \epsilon_{kif}r_fn_l-\int_{Fluid}{\rm d}V({\varPi}_{lk}-{\varPi}^{U}_{lk}) \frac{\partial b_{ki}}{\partial r_l}, \end{equation}

and using ${G}^{Elastic}_i=\int _{\boldsymbol {r}_{p}} {\rm d}S\,{\varPi }_{lk} \epsilon _{kif}r_fn_l$,

(2.28)\begin{equation} {G}^{PIST}_i+{G}^{Elastic}_i=\int_{\boldsymbol{r}_{p}} {\rm d}S\,{\varPi}^{U}_{lk} \epsilon_{kif}r_fn_l-\int_{Fluid} {\rm d}V({\varPi}_{lk}-{\varPi}^{U}_{lk})\frac{\partial b_{ki}}{\partial r_l}.\end{equation}

The first term on the right-hand side is the torque on a particle about its centre of mass due to (constant) undisturbed polymer stress acting on its surface. It is zero by symmetry for a constant density fore-and-aft and axisymmetric particle such as a slender prolate spheroid, and hence for such a particle,

(2.29)\begin{equation} {\boldsymbol{G}^{PIST}}+{\boldsymbol{G}^{Elastic}}=-\int_{Fluid} {\rm d}V (\boldsymbol{\varPi}-\boldsymbol{\varPi}^{U}):\boldsymbol{\nabla}\boldsymbol{b}\end{equation}

and

(2.30)\begin{equation} {\boldsymbol{G}^{Elastic}}=-\int_{Fluid}{\rm d}V\,\boldsymbol{\nabla}\boldsymbol{\cdot} ((\boldsymbol{\varPi}-\boldsymbol{\varPi}^{U})\boldsymbol{\cdot}\boldsymbol{b}). \end{equation}

Equation (2.28) (or (2.29) for a fore-and-aft and axisymmetric particle) represents the total torque acting on the particle due to the presence of the polymers (including effects from both the polymeric stress and the polymer-induced solvent stress), as a function of polymer stress, $\boldsymbol {\varPi }$, and its undisturbed value, $\boldsymbol {\varPi }^U$. For the Oldroyd-B fluid in a simple shear flow with 1, 2 and 3 as flow, gradient and vorticity directions, respectively, in Cartesian coordinates, ${\varPi }^{U}_{ij}=c(2De\delta _{i1}\delta _{j1}+(\delta _{i2}\delta _{j1}+\delta _{i1}\delta _{j2}))$.

4.1. Flow of Newtonian fluid around a slender fibre and particle rotation rates

Due to the relative velocity between the particle's centreline and the imposed flow, the flow disturbance at $O(1/\log (\kappa ))$ and higher orders in $1/\log (\kappa )$ can be considered to be generated by a line of point forces or force Stokeslets located at the particle centreline. This flow is considered by the general SBT of Cox (Reference Cox1970) up to $O(1/\log (\kappa )^2)$. An equivalent theory by Batchelor (Reference Batchelor1970) is valid at all orders in $1/\log (\kappa )$ for a prolate spheroid. For quantitative accuracy, it is advantageous to use the Stokeslet distribution, defined as $\boldsymbol {h}^{(0)}$ below, from the SBT of Batchelor (Reference Batchelor1970) instead of Cox (Reference Cox1970). For a prolate spheroidal particle fixed in the flow–vorticity plane of a simple shear flow or rotating about its centreline in a quiescent fluid, the force Stokeslet from these two theories is zero. The flow driven by the local velocity gradients, which acts at $O(1/\kappa ^2)$, is dominant in these cases. Considering the velocity gradients in the far-field relative to the particle, Cox (Reference Cox1971) provides the solution for this flow. It is a combination of flows driven by $O(1/\kappa ^2)$ force Stokeslets ($\boldsymbol {f}^{(0)}$) and doublets ($\boldsymbol {\mathcal {G}}^{(0)}$).

The Newtonian or the leading-order (in $c$) outer flow generated by a fibre freely rotating with an angular velocity $\boldsymbol {\omega }^{(0)}_p$ in an imposed shear flow of Newtonian fluid with 1, 2 and 3 as the flow, gradient and vorticity directions, when in orientation $\boldsymbol {p}$,

(4.2)\begin{equation} \boldsymbol{p}=\left[\,p_1\quad p_2\quad p_3\right]^{\rm T}, \end{equation}

close to the flow–vorticity plane is obtained by superposition of flows from Batchelor (Reference Batchelor1970), and Cox (Reference Cox1971) with the assumption $p_2\ll 1$ and $\kappa \gg 1$. The outer flow in the presence of the fibre is

(4.3)\begin{equation} \boldsymbol{u}_{out}^{{M}{(0)}}(\boldsymbol{r};\boldsymbol{\omega}^{(0)}_p,\boldsymbol{p})=\boldsymbol{r}\boldsymbol{\cdot} (\boldsymbol{\nabla}\boldsymbol{u})_\infty+{\boldsymbol{u}_{out}^{{M}{(0)^{\prime}}}} (\boldsymbol{r};\boldsymbol{\omega}^{(0)}_p,\boldsymbol{p}).\end{equation}

The disturbance velocity field,

(4.4)\begin{align} {\boldsymbol{u}_{out}^{{M}{(0)^{\prime}}}}(\boldsymbol{r};\boldsymbol{\omega}^{(0)}_p,\boldsymbol{p})&= \frac{1}{8{\rm \pi}}\int_{-1}^1{\rm d}\lambda\left\{\frac{\boldsymbol{h}^{(0)}(\lambda;\boldsymbol{\omega}^{(0)}_p, \boldsymbol{p})}{\log(2\kappa)-1.5}+\frac{\boldsymbol{f}^{(0)}(\lambda;\boldsymbol{p})}{\kappa^2}\right\}\nonumber\\ &\quad \boldsymbol{\cdot}\left\{\frac{\boldsymbol{I}}{|\boldsymbol{r}-\lambda\boldsymbol{p}|}+\frac{(\boldsymbol{r}-\lambda \boldsymbol{p})(\boldsymbol{r}-\lambda\boldsymbol{p})}{|\boldsymbol{r}-\lambda\boldsymbol{p}|^3}\right\}\nonumber\\ &\quad +\frac{1}{8{\rm \pi}}\int_{-1}^1{\rm d}\lambda\,\frac{1}{\kappa^2}\boldsymbol{\mathcal{G}}^{(0)}( \lambda;\boldsymbol{\omega}^{(0)}_p,\boldsymbol{p}):\boldsymbol{\nabla} \left\{\frac{\boldsymbol{I}}{|\boldsymbol{r}-\lambda\boldsymbol{p}|}+\frac{(\boldsymbol{r}-\lambda \boldsymbol{p})(\boldsymbol{r}-\lambda\boldsymbol{p})}{|\boldsymbol{r}-\lambda\boldsymbol{p}|^3}\right\}, \end{align}

is the sum of flows due to the point forces distributions $\boldsymbol {h}^{(0)}$ and $\boldsymbol {f}^{(0)}$, and the force doublet distribution $\boldsymbol {\mathcal {G}}^{(0)}$ discussed above. The different source distributions are

(4.5)\begin{align} \left.\begin{gathered} \boldsymbol{h^{(0)}}(\lambda;\boldsymbol{\omega}^{(0)}_p,\boldsymbol{p})=-4 {\rm \pi}\lambda [\boldsymbol{p}\boldsymbol{\cdot}(\boldsymbol{\nabla} \boldsymbol{u})_\infty-\boldsymbol{\omega}^{(0)}_p\times\boldsymbol{p}]\boldsymbol{\cdot} \left(\frac{1}{2}\boldsymbol{pp}+(\boldsymbol{I}-\boldsymbol{p}\boldsymbol{p})\frac{\log(2\kappa)-1.5}{ \log(2\kappa)-0.5}\right),\\ \boldsymbol{f}^{(0)}(\lambda;\boldsymbol{p})=-4{\rm \pi}\lambda \left[0\quad p_1\quad 0\right]^{\rm T} \left(1-\frac{1}{\log(2\kappa)-0.5}\right)+O(\,p_2),\\ \boldsymbol{\mathcal{G}}^{(0)}(\lambda;\boldsymbol{\omega}^{(0)}_p,\boldsymbol{p})=2{\rm \pi} (1-\lambda^2)\left(\begin{bmatrix} 0 & 1+\dfrac{p_3^2}{2} & 0\\ \dfrac{p_3^2}{2} & 0 & -\dfrac{1}{2}p_1p_3\\ 0 & -\dfrac{1}{2}p_1p_3 & 0 \end{bmatrix}+(\boldsymbol{\omega}^{(0)}_p\boldsymbol{\cdot} \boldsymbol{p})\boldsymbol{\epsilon}\boldsymbol{\cdot}\boldsymbol{p}\right)+O(\,p_2). \end{gathered}\right\} \end{align}

A torque-free spheroid with aspect ratio, $\kappa$ rotating in a simple shear flow has the exact result for the temporal evolution of the orientation vector (Jeffery Reference Jeffery1922; Kim & Karrila Reference Kim and Karrila2013),

(4.6)\begin{equation} \dot{\boldsymbol{p}}^{(0)}=\boldsymbol{\omega}^{(0)}_p\times\boldsymbol{p}=\boldsymbol{\omega}_\infty \times\boldsymbol{p}+\frac{\kappa^2+1}{\kappa^2-1}(\boldsymbol{E}_\infty\boldsymbol{\cdot} \boldsymbol{p})\boldsymbol{\cdot}(\boldsymbol{I}-\boldsymbol{pp}), \end{equation}

where $\boldsymbol {\omega }_\infty$ and $\boldsymbol {E}_\infty$ are the vorticity vector and strain rate tensor of the imposed simple shear,

(4.7a,b)\begin{equation} \boldsymbol{\omega}_\infty=-\tfrac{1}{2}\left[ 0\quad 0\quad 1 \right]^{\rm T},\quad \boldsymbol{E}_\infty=\tfrac{1}{2}\begin{bmatrix}0 & 1 & 0\\1 & 0 & 0\\0 & 0 & 0\end{bmatrix}. \end{equation}

For a slender spheroid close to the flow–vorticity plane,

(4.8)\begin{equation} \boldsymbol{\omega}^{(0)}_p\boldsymbol{\cdot}\boldsymbol{p}=-\frac{p_3}{2}+O(\,p_2).\end{equation}

Using $p_3=1$ in the above equation, we obtain the log-rolling velocity (when the particle is oriented with the vorticity axis) of the particle in a Newtonian fluid and find that in this orientation, a slender particle rotates at the angular velocity of the fluid, i.e. half of the shear rate. $\boldsymbol {\omega }^{(0)}_p\boldsymbol{\cdot}\boldsymbol {p}$ from (4.8) and $\boldsymbol {\omega }^{(0)}_p\times \boldsymbol {p}$ from (4.6) are used to obtain $\boldsymbol {\mathcal {G}}^{(0)}$ and $\boldsymbol {h^{(0)}}$, respectively, in (4.5).

The flow due to $\boldsymbol {h}^{(0)}$ is taken from Batchelor (Reference Batchelor1970). It is obtained by an expansion in $1/\log (\kappa )$ for a slender particle. But for a slender prolate spheroid, this expansion terminates such that the flow due to $\boldsymbol {h}^{(0)}$ captures the disturbance created by the prolate spheroid at all orders in $1/\log (\kappa )$ when expressed as in (4.4) and (4.6). The flow generated at the next order in $\kappa$ arises at $O(\kappa ^{-2})$ and is generated by $\boldsymbol {f}^{(0)}$ and $\boldsymbol {\mathcal {G}}^{(0)}$. It is taken from Cox (Reference Cox1971) where only the $O(\kappa ^{-2})$ flow is available. Thus, the disturbance velocity in the outer region given by (4.4) has an overall error of $O(\kappa ^{-3})$. If only $\boldsymbol {h}^{(0)}$ is considered, the error is of $O(\kappa ^{-2})$. The flow generated by $\boldsymbol {h}^{(0)}$ is proportional to $p_2$ and therefore, no flow is produced by $\boldsymbol {h}^{(0)}$ when the particle is aligned in the flow–vorticity plane. In this plane $\boldsymbol {f}^{(0)}$ and $\boldsymbol {\mathcal {G}}^{(0)}$ capture the (highest) $O(\kappa ^{-2})$ disturbance created by the particle. Accounting for the flow generated by $\boldsymbol {f}^{(0)}$ and $\boldsymbol {\mathcal {G}}^{(0)}$ allows us to consider the influence of elasticity within the flow–vorticity plane.

As mentioned just before (4.1) we expect most of the changes in the particle's rotation rate due to elasticity to arise when the particle is near the flow–vorticity plane, i.e. when $|p_2|\le O(1/\kappa )$. Therefore, in $\boldsymbol {h}^{(0)}$ we only consider the flow at $O(\,p_2)$, i.e. a flow of $O(\,p_2/\log (\kappa ))$ with an error of $O(\,p_2^2/\log (\kappa ))$. Since the primary purpose of using the $\boldsymbol {f}^{(0)}$ and $\boldsymbol {\mathcal {G}}^{(0)}$ flow is to capture the finite effect of elasticity when the particle is in the flow–vorticity plane, we use $\boldsymbol {f}^{(0)}$ and $\boldsymbol {\mathcal {G}}^{(0)}$ with $p_2=0$ (as expressed in (4.5)) which leads to a flow of $O(1/\kappa ^2)$ with an error of $O(\,p_2/\kappa ^2)$. When $|p_2|/\log (\kappa )$ is more than $1/\kappa ^2$, $\boldsymbol {h}^{(0)}$ driven flow dominates. In the $|p_2|\le O(1/\kappa )$ regime considered, the errors in $\boldsymbol {h}^{(0)}$ driven flow, i.e. $O(\,p_2^2/\log (\kappa ))$, are always smaller than the flow generated by $\boldsymbol {f}^{(0)}$ and $\boldsymbol {\mathcal {G}}^{(0)}$, i.e. $O(1/\kappa ^2)$. As $p_2$ approaches zero, the $\boldsymbol {h}^{(0)}$ driven flow and the associated errors fall rapidly to zero making $\boldsymbol {f}^{(0)}$ and $\boldsymbol {\mathcal {G}}^{(0)}$ driven flow at $O(1/\kappa ^2)$ the dominating one. Hence, the error terms arising from either of the $\boldsymbol {h}^{(0)}$ or $\boldsymbol {f}^{(0)}$ and $\boldsymbol {\mathcal {G}}^{(0)}$ driven flow are always lower than the actual flow when $|p_2|\le O(1/\kappa )$. For $|p_2|>O(1/\kappa )$, the Newtonian rotation rate of the fibre dominates over the changes due to elasticity, and we consider the exact Jeffery rotation to account for it.

We have considered spheroidal particles in our theory. However, in experiments with slender particles (Gauthier et al. Reference Gauthier, Goldsmith and Mason1971; Bartram et al. Reference Bartram, Goldsmith and Mason1975; Iso et al. Reference Iso, Cohen and Koch1996a,Reference Iso, Koch and Cohenb) it is convenient to fabricate cylindrical particles. The forces generated by the blunt ends of a slender cylinder lead to an additional torque of $O(1/\kappa ^2)$ (Cox Reference Cox1971) rendering the $O(1/\kappa ^2)$ flow generated by $\boldsymbol {f}^{(0)}$ and $\boldsymbol {\mathcal {G}}^{(0)}$ inaccurate. Instead of being valid at all orders in $1/\log (\kappa )$, the $\boldsymbol {h}^{(0)}$ generated flow may be considered up to order $1/\log (\kappa )^3$ (Batchelor Reference Batchelor1970). Taking the $O(\,p_2)$ terms with this velocity field can still allow us to consider the flow accurately up to $O(\,p_2/\log (\kappa ))$. Thus, the upcoming theoretical development for the orientation dynamics of a slender prolate spheroid leading up to (4.46) can still be used for a slender cylinder while ignoring the $\boldsymbol {f}^{(0)}$ and $\boldsymbol {\mathcal {G}}^{(0)}$ flow and using only the appropriate terms up to $1/\log (\kappa )^3$ instead of the factor $1/(2\log (2\kappa )-3)$ appearing in the following text. Most of the features described in § 5 while analysing the theoretical prediction of the influence of viscoelasticity on the orientation of a slender spheroid will be qualitatively valid for a slender cylinder or a general slender particle.

The Newtonian flow at $O(c)$ is due to a fibre rotating with the perturbed angular velocity, $\boldsymbol {\omega }^{(1)}_p$, in a quiescent fluid, i.e. (3.18a,b)–(3.19). In the outer region, the fluid velocity is

(4.9) \begin{align} {{\boldsymbol{u}_{out}^{M}}^{(1)}}(\boldsymbol{r};\boldsymbol{\omega}^{(1)}_p,\boldsymbol{p})& =\frac{1}{8{\rm \pi}}\int_{-1}^1{\rm d}\lambda\,\frac{\boldsymbol{h}^{(1)}(\lambda;\boldsymbol{\omega}^{(1)}_p, \boldsymbol{p})}{\log(2\kappa)-0.5}\boldsymbol{\cdot} \left\{\frac{\boldsymbol{I}}{|\boldsymbol{r}-\lambda\boldsymbol{p}|}+\frac{(\boldsymbol{r}-\lambda \boldsymbol{p})(\boldsymbol{r}-\lambda\boldsymbol{p})}{|\boldsymbol{r}-\lambda\boldsymbol{p}|^3}\right\}\nonumber\\ &\quad +\frac{1}{8{\rm \pi}}\int_{-1}^1{\rm d}\lambda\,\frac{1}{\kappa^2}\boldsymbol{\mathcal{G}}^{(1)}( \lambda;\boldsymbol{\omega}^{(1)}_p,\boldsymbol{p}):\boldsymbol{\nabla} \left\{\frac{\boldsymbol{I}}{|\boldsymbol{r}-\lambda\boldsymbol{p}|}+\frac{(\boldsymbol{r}-\lambda \boldsymbol{p})(\boldsymbol{r}-\lambda\boldsymbol{p})}{|\boldsymbol{r}-\lambda\boldsymbol{p}|^3}\right\}, \end{align}

where

(4.10a,b)\begin{align} \boldsymbol{h^{(1)}}(\lambda;\boldsymbol{\omega}^{(1)}_p,\boldsymbol{p})=4{\rm \pi}\lambda \boldsymbol{\omega}^{(1)}_p\times\boldsymbol{p},\quad \boldsymbol{\mathcal{G}}^{(1)}(\lambda;\boldsymbol{\omega}^{(1)}_p,\boldsymbol{p})=-2{\rm \pi} (1-\lambda^2)(\boldsymbol{\omega}^{(1)}_p\boldsymbol{\cdot} \boldsymbol{p})\boldsymbol{\epsilon}\boldsymbol{\cdot}\boldsymbol{p}.\end{align}

The torque generated by this flow is

(4.11)\begin{align} &{\boldsymbol{G}^{MIST}}^{(1)}(\boldsymbol{\omega}^{(1)}_p,\boldsymbol{p})\nonumber\\ &\quad =\int_{-1}^1{\rm d}\lambda\left(\frac{1}{\kappa^2}\boldsymbol{\boldsymbol{\epsilon}}:\boldsymbol{\mathcal{G}}^{(1)}- \frac{\lambda\boldsymbol{p}\times\boldsymbol{h}^{(1)}}{\log(2\kappa)-0.5}\right)=\frac{8{\rm \pi}}{3} \left(\frac{(\boldsymbol{\omega}^{(1)}_p\boldsymbol{\cdot}\boldsymbol{p}) \boldsymbol{p}}{\kappa^2}-\frac{\boldsymbol{\omega}^{(1)}_p\boldsymbol{\cdot} (\boldsymbol{I}-\boldsymbol{pp})}{\log(2\kappa)-0.5}\right). \end{align}

Taking the cross product of this equation with orientation vector $\boldsymbol {p}$ leads to the $O(c)$ rotation rate

(4.12)\begin{equation} \dot{\boldsymbol{p}}^{(1)}=\boldsymbol{\omega}^{(1)}_p\times\boldsymbol{p}=-\frac{3}{8{\rm \pi}}({\log(2\kappa)-0.5}) ({\boldsymbol{G}^{MIST}}^{(1)}(\boldsymbol{\omega}^{(1)}_p,\boldsymbol{p})\times\boldsymbol{p}). \end{equation}

Using the torque balance at $O(c)$ from (3.17) we obtain

(4.13)\begin{equation} \dot{\boldsymbol{p}}^{(1)}=\dot{\boldsymbol{p}}^{(1)}_{PIST}+\dot{\boldsymbol{p}}^{(1)}_{Elastic},\end{equation}

where

(4.14)\begin{equation} \dot{\boldsymbol{p}}^{(1)}_{PIST}=\frac{3}{8{\rm \pi}}({\log(2\kappa)-0.5})({\boldsymbol{G}^{PIST}}^{(1)}\times\boldsymbol{p}) \end{equation}

is the effect of polymer-induced solvent stresses on the particle rotation rate, and

(4.15)\begin{equation} \dot{\boldsymbol{p}}^{(1)}_{Elastic}=\frac{3}{8{\rm \pi}}({\log(2\kappa)-0.5})({\boldsymbol{G}^{Elastic}}^{(1)} \times\boldsymbol{p})\end{equation}

is the effect of the elastic stress on the particle rotation rate.

4.2. Rotation due to polymer-induced solvent stress

From (2.4), (3.14), (3.16) and (4.14), the rotation rate due to the polymer-induced solvent stress is

(4.16)\begin{equation} \dot{\boldsymbol{p}}^{(1)}_{PIST}=-\frac{3}{8{\rm \pi} De}({\log(2\kappa)-0.5})\boldsymbol{p}\times\int_{Fluid} {\rm d}V(\boldsymbol{\nabla}\boldsymbol{\cdot}({\boldsymbol{\varLambda}^{(0)}}-{\boldsymbol{\varLambda}^U})) \boldsymbol{\cdot}\boldsymbol{b}.\end{equation}

We remind the reader that here $\boldsymbol {b}$ is the auxiliary ‘velocity’ field used in the reciprocal theorem for deriving the polymer-induced solvent torque ((2.24) or (3.16)). Here $\boldsymbol {b}\boldsymbol{\cdot}\boldsymbol {\omega }$ corresponds to the fluid velocity around a fibre rotating with an angular velocity $\boldsymbol {\omega }$ in a quiescent fluid. Therefore,

(4.17)\begin{equation} \boldsymbol{b}=\boldsymbol{\nabla}_{\boldsymbol{\omega}^{(1)}_p} {\boldsymbol{u}^{M}}^{(1)}, \end{equation}

where ${\boldsymbol {u}^{M}}^{(1)}$ is the solution to (3.18a,b) and (3.19). The volume integral in (4.16) is approximated from the outer region of slender body theory, where the particle is replaced with a line of point forces and doublets. We find the integral from the outer region to converge, and the contribution from the inner region is expected to be small due to the smaller volume of that region. Therefore,

(4.18)\begin{equation} \dot{\boldsymbol{p}}^{(1)}_{PIST}\approx-\frac{3}{8{\rm \pi} De}({\log(2\kappa)-0.5})\boldsymbol{p}\times\int {\rm d}^3\boldsymbol{r}\,\boldsymbol{\nabla}\boldsymbol{\cdot}({\boldsymbol{\varLambda}_{out}^{(0)}}-{ \boldsymbol{\varLambda}^U})\boldsymbol{\cdot}\boldsymbol{b}_{out},\end{equation}

where the volume integral is taken over the entire space, ${\boldsymbol {\varLambda }_{out}^{(0)}}$ is the polymer conformation in the outer region and from (4.9), (4.10a,b) and (4.17),

(4.19)\begin{equation} \boldsymbol{b}_{out}=\frac{1}{8{\rm \pi}}\int_{-1}^1{\rm d}\lambda\, \frac{4{\rm \pi}\lambda\boldsymbol{p}\boldsymbol{\cdot}\boldsymbol{\epsilon}}{\log(2\kappa)-0.5} \boldsymbol{\cdot}\left(\frac{\boldsymbol{I}}{|\boldsymbol{r}-\lambda\boldsymbol{p}|}+\frac{(\boldsymbol{r}-\lambda \boldsymbol{p})(\boldsymbol{r}-\lambda \boldsymbol{p})}{|\boldsymbol{r}-\lambda\boldsymbol{p}|^3}\right)+O(\kappa^{-2}). \end{equation}

It is straightforward to show

(4.20)\begin{equation} \dot{\boldsymbol{p}}^{(1)}_{PIST}\approx\tfrac{3}{2}(\boldsymbol{I}-\boldsymbol{p}\boldsymbol{p}) \boldsymbol{\cdot}\int_{-1}^1{\rm d}\lambda\,\lambda\boldsymbol{u}_{out}^{polymer}(\lambda\boldsymbol{p}), \end{equation}

where

(4.21)\begin{align} \boldsymbol{u}_{out}^{polymer}(\lambda\boldsymbol{p})=\frac{1}{8{\rm \pi}}\int {\rm d}^3\boldsymbol{r}\,\frac{1}{De}\boldsymbol{\nabla}\boldsymbol{\cdot}({\boldsymbol{\varLambda}_{out}^{(0)}}-{ \boldsymbol{\varLambda}^U})\boldsymbol{\cdot}\left(\frac{\boldsymbol{I}}{|\boldsymbol{r}-\lambda \boldsymbol{p}|}+\frac{(\boldsymbol{r}-\lambda\boldsymbol{p})(\boldsymbol{r}-\lambda \boldsymbol{p})}{|\boldsymbol{r}-\lambda\boldsymbol{p}|^3}\right) \end{align}

is the velocity field evaluated on the particle centreline that is produced by polymeric force, $({1}/{De})\boldsymbol {\nabla }\boldsymbol{\cdot}({\boldsymbol {\varLambda }_{out}^{(0)}}-{\boldsymbol {\varLambda }^U})$, acting in the outer region. This form of the rotation rate, (4.20) and (4.21), was considered by Harlen & Koch (Reference Harlen and Koch1993). However, they did not account for all the relevant terms in calculating ${\boldsymbol {\varLambda }_{out}^{(0)}}$ as we show below.

As discussed in § 3 and shown by (3.13) and (4.1), the polymer conformation, ${\boldsymbol {\varLambda }^{(0)}}$ depends on the leading-order Newtonian velocity. In the outer region, the velocity disturbance created by the particle is $O(\max [\,p_2/\log (\kappa ),1/\kappa ^2])$ smaller than the velocity of the imposed simple shear ((4.3), (4.4) and (4.5)). Thus, we linearize (4.1) in the outer region about the imposed flow field, $\boldsymbol {r}\boldsymbol{\cdot}(\boldsymbol {\nabla } \boldsymbol {u})_\infty$, and obtain the governing equation for the disturbance in the polymer conformation from its undisturbed value,

(4.22)\begin{equation} {\boldsymbol{\varLambda}_{out}^{(0)^{\prime}}}=\boldsymbol{\varLambda}_{out}^{(0)}-\boldsymbol{\varLambda}_{out}^U, \end{equation}

to be

(4.23) \begin{align} \left(\boldsymbol{r}\boldsymbol{\cdot}(\boldsymbol{\nabla} \boldsymbol{u})_\infty\boldsymbol{\cdot}\boldsymbol{\nabla}+\frac{1}{De}\right){\boldsymbol{\varLambda}_{out}^{(0)^{\prime}}}&- (\boldsymbol{\nabla} \boldsymbol{u})_\infty^T\boldsymbol{\cdot}{\boldsymbol{\varLambda}_{out}^{(0)^{\prime}}}-{\boldsymbol{\varLambda}_{out}^{(0)^{\prime}}} \boldsymbol{\cdot}(\boldsymbol{\nabla}\boldsymbol{u})_\infty\nonumber\\ &=(\boldsymbol{\nabla}{\boldsymbol{u}_{out}^{{M}{(0)^{\prime}}}}^{\rm T}\boldsymbol{\cdot} \boldsymbol{\varLambda}^U+\boldsymbol{\varLambda}^U\boldsymbol{\cdot}\boldsymbol{\nabla} {\boldsymbol{u}_{out}^{{M}{(0)^{\prime}}}}. \end{align}

It is more convenient to solve (4.23) in Fourier space,

(4.24)\begin{align} \left(-k_1\frac{\partial}{\partial k_2}+\frac{1}{De}\right){\hat{\boldsymbol{\varLambda}}_{out}^{(0)^{{\prime}}}}-(\boldsymbol{\nabla} \boldsymbol{u})_\infty^T\boldsymbol{\cdot}{\hat{\boldsymbol{\varLambda}}_{out}^{(0)^{{\prime}}}} & -{\hat{\boldsymbol{\varLambda}}_ {out}^{(0)^{{\prime}}}}\boldsymbol{\cdot}(\boldsymbol{\nabla}\boldsymbol{u})_\infty\nonumber\\ &={\rm i}({{\hat{\boldsymbol{u}}_{out}^{{M}^{(0)^{{\prime}}}}}}\boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{\varLambda}^U+\boldsymbol{\varLambda}^U\boldsymbol{\cdot}\boldsymbol{k}{{\hat{\boldsymbol{u}}_{out}^{{M}^{(0)^{{\prime}}}}}}), \end{align}

where ${\hat {\boldsymbol {\varLambda }}_{out}^{(0)^{\prime}}}=\mathcal {F}({{\boldsymbol {\varLambda }}_{out}^{(0)^{\prime}}})$ and ${{\hat {\boldsymbol {u}}_{out}^{M^{(0)^{\prime}}}}}=\mathcal {F}({{{\boldsymbol {u}}_{out}^{M^{(0)^{\prime}}}}})$ are the Fourier transforms of the disturbance of polymer conformation and fluid velocity in the outer region. The rotation rate in (4.20) is expressed as

(4.25)\begin{equation} \dot{\boldsymbol{p}}^{(1)}_{PIST}\approx\tfrac{3}{2}(\boldsymbol{I}-\boldsymbol{p}\boldsymbol{p}) \boldsymbol{\cdot}\int_{-1}^1{\rm d}\lambda\int {\rm d}^3\boldsymbol{k}\,\lambda\exp(2{\rm \pi} {\rm i}\lambda\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p})\hat{\boldsymbol{u}}_{out}^{polymer},\end{equation}

where from the convolution theorem,

(4.26)\begin{equation} \hat{\boldsymbol{u}}_{out}^{polymer}={\rm i}\frac{1}{De}\boldsymbol{k}\boldsymbol{\cdot} {\hat{\boldsymbol{\varLambda}}_{out}^{(0)^{{\prime}}}}(\boldsymbol{k};\boldsymbol{p})\boldsymbol{\cdot} \hat{\boldsymbol{J}}(\boldsymbol{k}),\quad\text{with } \hat{\boldsymbol{J}}(\boldsymbol{k})=\mathcal{F}\left(\frac{1}{8{\rm \pi}} \left(\frac{\boldsymbol{I}}{r}-\frac{\boldsymbol{rr}}{r^3}\right)\right)= \frac{1}{k^2}\left(\boldsymbol{I}-\frac{\boldsymbol{kk}}{k^2}\right). \end{equation}

We solve the polymer constitutive equation (4.24) in Fourier space using the method of characteristics to obtain ${\hat {\boldsymbol {\varLambda }}_{out}^{(0)^{\prime}}}$. The characteristics are the streamlines of the imposed shear flow (in the Fourier space) and so only integration along the $k_2$ direction is needed. The limits of the integral are $\infty \textrm {sgn}(k_1)$ and $k_2$ and the boundary condition is ${\hat {\boldsymbol {\varLambda }}_{out}^{(0)^{\prime}}}(\infty \textrm {sgn}(k_1))=0$. We use computer algebra for this calculation. Here ${\hat {\boldsymbol {\varLambda }}_{out}^{(0)^{\prime}}}$, hence obtained, allows us to evaluate the integral in (4.25). The expression used for the disturbance in the Newtonian velocity field, ${\boldsymbol {u}_{out}^{{M}{(0)^{\prime}}}}$, has errors of $O(\textrm {max}[\,p_2^2/\log (\kappa ),p_2/\kappa ^2,1/\kappa ^3])$. Due to the linearization of the constitutive equation, the error in evaluation of rotation rate is $O(\epsilon _{\dot {p}})$, where

(4.27)\begin{equation} \epsilon_{\dot{p}}=\text{max}\left[\frac{p_2^2}{\log(\kappa)},\frac{p_2}{\kappa^2},\frac{1}{\kappa^3}\right].\end{equation}

Due to fortuitous cancelling of terms, a simple expression for the second (gradient direction) component of the rotation rate is obtained,

(4.28)\begin{equation} \dot{{p}}^{(1)}_{2,PIST}=-\frac{De p_1^2 p_2}{2\log(2\kappa)-3}+O(\epsilon_{\dot{p}}),\end{equation}

valid at all $De$. Equivalently, from the outer region integral of (3.16), i.e.

(4.29)\begin{equation} {\boldsymbol{G}^{PIST}}^{(1)}\approx \frac{1}{De}\int {\rm d}^3\boldsymbol{r}\,\boldsymbol{\nabla}\boldsymbol{\cdot}({\boldsymbol{\varLambda}_{out}^{(0)}}- {\boldsymbol{\varLambda}^U})\boldsymbol{\cdot}\boldsymbol{b}_{out}, \end{equation}

we obtain the first and third components of the polymer-induced solvent torque,

(4.30)\begin{equation} \frac{{{G}_1^{PIST}}^{(1)}}{p_3}=-\frac{{{G}_3^{PIST}}^{(1)}}{p_1}=\frac{4De^2 p_1^2 p_2}{3(\log(2\kappa)-0.5)(\log(2\kappa)-1.5)}+O\left(\frac{\epsilon_{\dot{p}}}{\log(\kappa)}\right). \end{equation}

Since, $\dot {\boldsymbol {p}}^{(1)}_{PIST}={3}/{(8{\rm \pi} )}({\log (2\kappa )-0.5})({\boldsymbol {G}^{PIST}}^{(1)} \times \boldsymbol {p})$ (4.14) we may obtain $\dot {{p}}^{(1)}_{2,PIST}=-De{p_1^2 p_2}/({2\log (2\kappa )-3})$. However, the equations for determining the torque component in the gradient direction, ${{G}_2^{PIST}}^{(1)}$ or rotation rates in flow and vorticity direction, $\dot {{p}}^{(1)}_{1,PIST}$ and $\dot {{p}}^{(3)}_{1,PIST}$ are not tractable for a general $De$. Therefore, to obtain these components we consider small and large $De$ limits separately in §§ 4.2.1 and 4.2.2.

We find the particle rotation rate due to polymer-induced solvent stress in gradient direction, $\dot {{p}}^{(1)}_{2, PIST}$ to be dependent on the first normal stress difference of the Oldroyd-B fluid. This is because $\dot {\boldsymbol {p}}^{(1)}_{PIST}$ in (4.25) is directly proportional to ${\hat {\boldsymbol {\varLambda }}_{out}^{(0)^{\prime}}}$ that is in turn dependent on $\boldsymbol {\varLambda }^U$ via (4.24). This first normal stress difference dependence is further confirmed by repeating the above calculation by artificially setting ${\varLambda }^U_{ij}=De^2\delta _{i1}\delta _{j1}$, such that the only non-zero non-Newtonian property of the fluid is a finite first normal stress difference equivalent to that of an Oldroyd-B fluid. Similarly, we find the two remaining rotation rate components in the limit of large $De$ calculation of § 4.2.1 to arise from the first normal stress difference of the Oldroyd-B fluid. The rotation rate for a second-order fluid ($O(De)$ rotation rate in the small $De$ calculation of § 4.2.2) is also due to the first normal stress difference. An Oldroyd-B fluid has no second normal stress difference. In addition to its appropriate modelling of the simple shear flow of a dilute polymeric liquid, the Oldroyd-B model the advantages over other models such as FENE-P or Giesekus (Bird et al. Reference Bird, Armstrong and Hassager1987) of simplicity and hence better analytical tractability. As mentioned in § 1, the second normal stress difference for most polymeric fluids is much smaller than the first normal stress difference. Hence, the effect of the first normal stress difference as ascertained from the Oldroyd-B model is likely to be the most important contribution in determining the influence of polymers on the orientational dynamics of a prolate spheroid in simple shear flow.

4.2.1. Large $De$

In the large $De$ regime, the relaxation of the disturbance in polymer conformation in the outer region is much slower than its convection and stretching by the imposed velocity field. Therefore, when $De\gg 1$, (4.24) simplifies to

(4.31)\begin{align} -k_1\frac{\partial}{\partial k_2}{\hat{\boldsymbol{\varLambda}}_{out}^{(0)^{\prime}}}-(\boldsymbol{\nabla} \boldsymbol{u})_\infty^T\boldsymbol{\cdot}{\hat{\boldsymbol{\varLambda}}_{out}^{(0)^{\prime}}}- {\hat{\boldsymbol{\varLambda}}_{out}^{(0)^{\prime}}}\boldsymbol{\cdot}(\boldsymbol{\nabla} \boldsymbol{u})_\infty={\rm i}({{\hat{\boldsymbol{u}}_{out}^{{M}^{(0)^{\prime}}}}} \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{\varLambda}^U+\boldsymbol{\varLambda}^U\boldsymbol{\cdot} \boldsymbol{k}{{\hat{\boldsymbol{u}}_{out}^{{M}^{(0)^{\prime}}}}}).\end{align}

This equation is solved in a similar way as (4.24) described earlier. Using ${\hat {\boldsymbol {\varLambda }}_{out}^{(0)^{\prime}}}$ obtained from this equation we find the particle rotation rate due to polymer-induced solvent stresses at large $De$ to be

(4.32)\begin{align} \lim_{De\gg 1}\dot{\boldsymbol{p}}^{(1)}_{PIST}=De\begin{bmatrix} -p_1p_3^2 \alpha_{De\gg 1}\\- \dfrac{p_1^2 p_2}{2\log(2\kappa)-3}\\ p_1^2p_3\alpha_{De\gg 1} \end{bmatrix},\quad \text{with } \alpha_{De\gg 1}=\frac{-4}{\kappa^2}+\frac{3{\rm \pi}}{4\log(2\kappa)-6}p_2\tan^{-1}\frac{p_1}{|p_3|}.\end{align}

The large $De$ approximation can be viewed as the leading $O(De)$ term in an expansion in powers of $1/De$. Thus, the error due to the expansion in $1/De$ is $O(1)$ in the first and third components of $\lim _{De\gg 1}\dot {\boldsymbol {p}}^{(1)}_{PIST}$. There is no error due to expansion in $1/De$ for $\dot {{p}}^{(1)}_{2,PIST}$ as same value is obtained via $De\gg 1$ approximation as that for a general $De$ in (4.28). As mentioned earlier, $\dot {{p}}^{(1)}_{2,PIST}$ has an error of $O(\epsilon _{\dot {p}})$ (4.27). Here $\dot {{p}}^{(1)}_{1,PIST}$ and $\dot {{p}}^{(1)}_{3,PIST}$ have an error of $O(De\epsilon _{\dot {p}})$. To the best of our knowledge, the only previous theoretical study concerning the rotation of a slender particle in a viscoelastic fluid at high $De$ was conducted by Harlen & Koch (Reference Harlen and Koch1993). However, they did not account for the stretching of the conformation disturbance, ${\hat {\boldsymbol {\varLambda }}_{out}^{(0)^{\prime}}}$, by the mean velocity gradients, i.e. they omitted the $(\boldsymbol {\nabla } \boldsymbol {u})_\infty ^T\boldsymbol{\cdot}{\hat {\boldsymbol {\varLambda }}_{out}^{(0)^{\prime}}}+ {\hat {\boldsymbol {\varLambda }}_{out}^{(0)^{\prime}}}\boldsymbol{\cdot}(\boldsymbol {\nabla } \boldsymbol {u})_\infty$ term in (4.31).

4.2.2. Small $De$

When $De\ll 1$, a solution of (4.23) or (4.24) is obtained via a regular perturbation expansion of ${\hat {\boldsymbol {\varLambda }}_{out}^{(0)^{\prime}}}$ in the powers of $De$

(4.33)\begin{equation} {{\hat{\boldsymbol{\varLambda}}_{out}^{(0)^{\prime}}}}=\varSigma_{n=0}^\infty De^n {{\hat{\boldsymbol{\varLambda}}_{out}^{(0)^{{\prime}^{(n)}}}}}.\end{equation}

The first two terms in this expansion are

(4.34)\begin{align} \left.\begin{gathered} {{\hat{\boldsymbol{\varLambda}}_{out}^{(0)^{{\prime}^{(0)}}}}}=\boldsymbol{0},\\ {{\hat{\boldsymbol{\varLambda}}_{out}^{(0)^{\prime^{(1)}}}}}={\rm i}(\boldsymbol{k}{\hat{\boldsymbol{u}}^{M^{(0)}}} +{\hat{\boldsymbol{u}}^{M^{(0)}}}\boldsymbol{k}). \end{gathered}\right\} \end{align}

The polymer conformation at higher order in $De$ is obtained from the following equations:

(4.35)\begin{equation} \left.\begin{aligned} {\hat{\boldsymbol{\varLambda}}_{out}^{(0^{\prime})^{(2)}}} & =(\boldsymbol{\nabla}\boldsymbol{u})_\infty^T \boldsymbol{\cdot}{\hat{\boldsymbol{\varLambda}}_{out}^{(0^{\prime})^{(1)}}}+{\hat{\boldsymbol{\varLambda}}_{out}^{(0^{\prime})^{(1)}}} \boldsymbol{\cdot}(\boldsymbol{\nabla}\boldsymbol{u})_\infty\\ & \quad +{\rm i}({\hat{\boldsymbol{u}}^{M^{(0)}}}\boldsymbol{k}\boldsymbol{\cdot} {\boldsymbol{\varLambda}^{U^{(1)}}}+{\boldsymbol{\varLambda}^{U^{(1)}}}\boldsymbol{\cdot} \boldsymbol{k}{\hat{\boldsymbol{u}}^{M^{(0)}}})+k_1\frac{\partial}{\partial k_2}{\hat{\boldsymbol{\varLambda}}_{out}^{(0^{\prime})^{(1)}}},\\ {\hat{\boldsymbol{\varLambda}}_{out}^{(0^{\prime})^{(3)}}} & = (\boldsymbol{\nabla}\boldsymbol{u})_\infty^T\boldsymbol{\cdot}{\hat{\boldsymbol{\varLambda}}_{out}^{(0^{\prime})^{(2)}}} + {\hat{\boldsymbol{\varLambda}}_{out}^{(0^{\prime})^{(2)}}}\boldsymbol{\cdot}(\boldsymbol{\nabla}\boldsymbol{u})_\infty\\ & \quad +{\rm i}({\hat{\boldsymbol{u}}^{M^{(0)}}}\boldsymbol{k}\boldsymbol{\cdot}{\boldsymbol{\varLambda}^{U^{(2)}}}+ {\boldsymbol{\varLambda}^{U^{(2)}}}\boldsymbol{\cdot}\boldsymbol{k}{\hat{\boldsymbol{u}}^{M^{(0)}}})+ k_1\frac{\partial}{\partial k_2}{\hat{\boldsymbol{\varLambda}}_{out}^{(0^{\prime})^{(2)}}},\\ {\hat{\boldsymbol{\varLambda}}_{out}^{(0^{\prime})^{(n)}}} & =( \boldsymbol{\nabla}\boldsymbol{u})_\infty^T\boldsymbol{\cdot}{\hat{\boldsymbol{\varLambda}}_{out}^{(0^{\prime})^{(n-1)}}} +{\hat{\boldsymbol{\varLambda}}_{out}^{(0^{\prime})^{(n-1)}}}\boldsymbol{\cdot} (\boldsymbol{\nabla}\boldsymbol{u})_\infty+k_1\frac{\partial}{\partial k_2}{\hat{\boldsymbol{\varLambda}}_{out}^{(0^{\prime})^{ (n-1)}}},\quad n\ge 4. \end{aligned}\right\}\end{equation}

Therefore, in the limit of low $De$, $\dot {\boldsymbol {p}}^{(1)}_{PIST}$ can be obtained up to arbitrary order in $De$,

(4.36)\begin{equation} \lim_{De\ll 1}\dot{\boldsymbol{p}}^{(1)}_{PIST}=-\frac{3}{8{\rm \pi} }({\log(2\kappa)-0.5})\boldsymbol{p}\times\varSigma_{n=1}^\infty De^{n-1} \int_{Fluid} {\rm d}V(\boldsymbol{\nabla}\boldsymbol{\cdot}({{{\boldsymbol{\varLambda}}_{out}^{(0)^{\prime}}}}^{(n)})) \boldsymbol{\cdot}\boldsymbol{b}. \end{equation}

We will discuss the contribution of the $n=1$ term, i.e. $O(1)$ term in $De$ expansion of $\dot {\boldsymbol {p}}^{(1)}_{PIST}$ in the above equation in § 4.4, along with the similar term in the expansion of $\dot {\boldsymbol {p}}^{(1)}_{Elastic}$. Higher-order contributions to $\dot {\boldsymbol {p}}^{(1)}_{PIST}$ can be evaluated in Fourier space. From (4.33)–(4.35), using the expansion of ${{\boldsymbol {\varLambda }}_{out}^{(0)}}'$ up to $O(De^{8})$ (even higher orders can be easily evaluated), we obtain

(4.37) \begin{align} \lim_{De\ll 1} \dot{\boldsymbol{p}}^{(1)}_{PIST}&=-\frac{3}{8{\rm \pi} }({\log(2\kappa)-0.5}) \boldsymbol{p}\times \int_{Fluid} {\rm d}V(\boldsymbol{\nabla}\boldsymbol{\cdot}({{{\boldsymbol{\varLambda}}_{out}^{(0)}}'}^{(1)})) \boldsymbol{\cdot}\boldsymbol{b}\nonumber\\ &\quad +De\begin{bmatrix} -p_1p_3^2 \alpha_{De\ll 1}\\ -\dfrac{p_1^2 p_2}{2\log(2\kappa)-3}\\ p_1^2p_3\alpha_{De\ll 1} \end{bmatrix}, \end{align}

where

(4.38)\begin{align} \alpha_{De\ll 1}=-\frac{4}{\kappa^2}+\frac{ p_1p_2 }{8\log(2\kappa)-12} \left(5 De-\frac{3}{2}(2+p_3^2)De^3+\frac{13}{10}(8+p_3^2(4+3p_3^2))De^5\right).\end{align}

As expected, we obtain the same expression for $\dot {{p}}^{(1)}_{2,PIST}$ from this perturbation expansion in $De$ valid at small $De$ as was obtained directly for a general $De$ in (4.28). The error in the rotation rate in (4.37) is $O(De\epsilon _{\dot {p}})$ for all the components ($\epsilon _{\dot {p}}$ is given in (4.27)).

4.3. Rotation due to elastic stress

The elastic torque, ${\boldsymbol {G}^{Elastic}}^{(1)}$, is due to the polymer stress on the particle surface ((3.14) and (3.16)). This is evaluated through the polymer constitutive equations on the particle surface written in the frame of reference moving with the particle surface. Due to the absence of polymer convection relative to the particle on the latter's surface, the quasisteady constitutive equation (4.1) simplifies to a set of coupled algebraic equations,

(4.39)\begin{equation} (\boldsymbol{\nabla}{\boldsymbol{u}_{in}^{M}}^{(0)}(\boldsymbol{r}_p; \boldsymbol{\omega}^{(0)}_p,\boldsymbol{p}))^{\rm T}\boldsymbol{\cdot} \boldsymbol{\varLambda}^{(0)}_p+\boldsymbol{\varLambda}^{(0)}_p \boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{u}_{in}^{M}}^{(0)}( \boldsymbol{r}_p;\boldsymbol{\omega}^{(0)}_p,\boldsymbol{p})+ \frac{1}{De}(\boldsymbol{I}-\boldsymbol{\varLambda}^{(0)}_p)=0, \end{equation}

at each point on the surface. Here $\boldsymbol {r}_p$ is the position vector of a point on the surface. $\boldsymbol {\nabla }{\boldsymbol {u}_{in}^{M}}^{(0)}(\boldsymbol {r}_p;\boldsymbol {\omega }^{(0)}_p,\boldsymbol {p})$ and $\boldsymbol {\varLambda }^{(0)}_p$ represent the surface velocity gradient of the inner velocity field and the surface polymer conformation, respectively, at $\boldsymbol {r}_p$. The inner velocity field is obtained from Cox (Reference Cox1970, Reference Cox1971). These velocity fields include the effects of point force Stokeslets and doublets that correspond to the outer velocity field accurate up to $O(\,p_2/\log (\kappa )^2)$ and $O(1/\kappa ^2)$. Solving (4.39) via an asymptotic expansion in $1/\kappa$ and ignoring higher-order terms leads to the elastic torque

(4.40)\begin{equation} {\boldsymbol{G}}^{(1)}_{Elastic}=\frac{8{\rm \pi} c}{3\kappa^2}\begin{bmatrix} p_1 p_3\\ 0\\-p_1^2 \end{bmatrix}+O(\max[\,p_2\kappa^{-2},\kappa^{-3}]),\end{equation}

and the corresponding rotation rate

(4.41)\begin{align} \dot{\boldsymbol{p}}^{(1)}_{Elastic}&=-\frac{3}{8{\rm \pi}}({\log(2\kappa)-0.5}) \boldsymbol{p}\times\int {\rm d}S\left(\boldsymbol{r}\times\frac{1}{De}(\boldsymbol{\varLambda}^{(0)}_p- \boldsymbol{I})\boldsymbol{\cdot}\boldsymbol{n}\right)\nonumber\\ &=\frac{{\log(2\kappa)-0.5}}{\kappa^2}p_1\begin{bmatrix} 0\\ 1\\0 \end{bmatrix}+O(\log(\kappa)\max[\,p_2\kappa^{-2},\kappa^{-3}]). \end{align}

Therefore, for all $De$, in limits of large $\kappa$ and small $p_2$, the $O(c)$ elastic torque is independent of the polymer relaxation time, $De$.

4.4. Net rotation due to the polymers

Combining the results of the previous two sections, we obtain the net change in the fibre's rotation rate due to polymers.

For $De\ll 1$, we can obtain the $O(1)$ term in the $De$ expansion of $\dot {\boldsymbol {p}}^{(1)}_{PIST}$ and $\dot {\boldsymbol {p}}^{(1)}_{Elastic}$ for all particle orientations, without resorting to the outer region approximation for the former. The first two terms in the polymer configuration, everywhere (inner and outer region) is similar to that mentioned earlier in (4.34), i.e.

(4.42)\begin{equation} {{{{\boldsymbol{\varLambda}}^{(0)}}}=\boldsymbol{I}+De(\boldsymbol{\nabla} {\boldsymbol{u}^{M}}^{(0)}+(\boldsymbol{\nabla} {\boldsymbol{u}^{M}}^{(0)})^{\rm T})+O(De^2).}\end{equation}

From the combined effect of polymers from (3.16), (4.13), (4.14) and (4.15), along with using the value of ${{{\boldsymbol {\varLambda }}^{(0)}}}$ from (4.42), divergence theorem, the leading-order momentum equation (i.e. $\boldsymbol {\nabla }{p^{P}}^{(0)}=\boldsymbol {\nabla }\boldsymbol{\cdot}[\boldsymbol {\nabla } {\boldsymbol {u}^{P}}^{(1)}+(\boldsymbol {\nabla } {\boldsymbol {u}^{P}}^{(1)})^\textrm {T}]$) and the auxiliary Stokes problem ((2.18a–c) and (2.19)), we obtain

(4.43) \begin{align} \lim_{De\ll 1}{\dot{\boldsymbol{p}}^{(1)}}_{Elastic}&+{\dot{\boldsymbol{p}}^{(1)}}_{PIST}\nonumber\\ &=\int_{\boldsymbol{r}_{p}} {\rm d}S\,\boldsymbol{r}\times[(\boldsymbol{\nabla} {\boldsymbol{u}^{M}}^{(0)}+(\boldsymbol{\nabla}{\boldsymbol{u}^{M}}^{(0)})^{\rm T}-{p^{M}}^{(0)} \boldsymbol{I})\boldsymbol{\cdot}\boldsymbol{n}]+O(De)\nonumber\\ &= \int_{\boldsymbol{r}_{p}} {\rm d}S\,\boldsymbol{r}\times {\boldsymbol{\tau}^{M}}^{(0)}+O(De)\boldsymbol{\cdot}\boldsymbol{n}=O(De). \end{align}

The $O(1)$ term in this sum is equivalent to the rotation due to the leading-order torque (in $c$) acting on the particle. This is the torque on a freely rotating particle in a Newtonian fluid and is hence zero. Thus, for a freely rotating particle for $De\ll 1$, the $O(1)$ rotation rates (in the $De$ expansion) due to the torque generated by the elastic and polymer induced solvent stress cancel each other and lead to no change in the net particle rotation. The change in the particle rotation due to the polymers arises at $O(De)$. Close to the flow–vorticity plane, the rotation rate due to the elastic torque was shown to be independent of $De$ for all $De$ in (4.41), which we have just shown balances the equivalent rotation rate due to polymer-induced solvent torque for $De\ll 1$. Thus, the $O(De)$ and higher effect of polymers on the particle's rotation rate arises entirely due to the polymer-induced solvent torque (4.37) and is given by

(4.44) \begin{equation} \lim_{De\ll 1}(\kern0.7pt\dot{\boldsymbol{p}}^{(1)}_{Elastic}+\dot{\boldsymbol{p}}^{(1)}_{PIST})=De\begin{bmatrix} -p_1p_3^2 \alpha_{De\ll 1}\\ -\dfrac{p_1^2 p_2}{2\log(2\kappa)-3}\\ p_1^2p_3\alpha_{De\ll 1} \end{bmatrix},\end{equation}

where $\alpha _{De\ll 1}$ is given in (4.38) and the errors are $O(De\epsilon _{\dot {p}})$ for all the components (where $\epsilon _{\dot {p}}$ is given in (4.27)).

For $De\gg 1$, the total effect of viscoelasticity on the particle rotation also arises mainly from the polymer-induced solvent stress. But, unlike the $De\ll 1$ regime, here the reason is that the effect of the elastic stress is $O(De)$ smaller. The net rotation rate due to the polymers is

(4.45)\begin{equation} \lim_{De\gg 1}(\kern0.7pt\dot{\boldsymbol{p}}^{(1)}_{Elastic}+\dot{\boldsymbol{p}}^{(1)}_{PIST})=De\begin{bmatrix} -p_1p_3^2 \alpha_{De\gg 1}\\ -\dfrac{p_1^2 p_2}{2\log(2\kappa)-3}\\ p_1^2p_3\alpha_{De\gg 1} \end{bmatrix},\end{equation}

where $\alpha _{De\gg 1}$ is given by (4.32). The error in $\dot {{p}}^{(1)}_{1,PIST}$ and $\dot {{p}}^{(1)}_{3,PIST}$ is of $O(De\epsilon _{\dot {p}})$ (4.27). Neglecting the elastic torque leads to an additional error of $O(\log (\kappa )/\kappa ^2)$ in the second component and $O(\log (\kappa )\max [\,p_2\kappa ^{-2},\kappa ^{-3}])$ in the first and third component of the net rotation rate due to the polymers (4.41).

4.5. Equations of motion of a freely rotating fibre in low $c$ viscoelastic fluid at all $De$

Equations (4.44) and (4.45) encompass our main result for the effect of viscoelasticity on the rotation of a particle suspended in simple shear flow at large and small $De$, respectively. In the limit of large $De$ and for a second-order fluid ($O(De)$ rotation rate in small $De$ limit), the effect of viscoelasticity is due to the first normal stress difference of the fluid as discussed in § 4.2. The governing equation for the orientation dynamics of a slender prolate spheroid that includes the viscoelastic effects near the flow–vorticity plane is

(4.46)\begin{equation} \dot{\boldsymbol{p}}=\frac{{\rm d}\boldsymbol{p}}{{\rm d}t}=\begin{bmatrix} \dfrac{p_2}{2}+\dfrac{\kappa^2-1}{\kappa^2+1}\left(\dfrac{p_2}{2}-p_1^2p_2\right)-c\,Dep_1p_3^2 \alpha\\ - \dfrac{p_1}{2}+\dfrac{\kappa^2-1}{\kappa^2+1}\left(\dfrac{p_1}{2}-p_1p_2^2\right)-c\,De \dfrac{p_1^2 p_2}{2\log(2\kappa)-3}\\ -\dfrac{\kappa^2-1}{\kappa^2+1}p_1 p_2 p_3+c\,Dep_1^2p_3\alpha \end{bmatrix},\end{equation}

where in the large $De$ limit $\alpha$ is $\alpha _{De\gg 1}$ (4.32). In the small $De$ limit we consider $\alpha _{De\ll 1}$ (4.38) up to $O(De)$ and interpolating between these two limits of $De$, we obtain the following uniformly valid approximation for $\alpha$ at all $De$:

(4.47)\begin{equation} \alpha=-\frac{4}{\kappa^2}+\frac{p_2}{4\log(2\kappa)-6}\frac{2.5p_1De}{1+2.5p_1De\left/\left(3{\rm \pi}\tan^{-1}\dfrac{p_1}{|p_3|}\right)\right.}.\end{equation}

Equation (4.46) introduces no errors in the Newtonian rotation rate of a prolate spheroidal particle. The slender body theory provides a good approximation for the Newtonian velocity field for $\kappa \gtrsim 10$. The errors in these equations due to the viscoelastic terms are of $O(c\,De\epsilon _{\dot {p}})$ for all the components ($\epsilon _{\dot {p}}$ is given in (4.27)) in the $De\ll 1$ limit. For $De\gg 1$, the errors are of $O(c\,De\epsilon _{\dot {p}})$ in the first and third component and $O(c\max [\log (\kappa )/\kappa ^2,\epsilon _{\dot {p}}])$ in the second component. Before analysing the influence of polymers on particle orientation as suggested by this equation we compare our theory at low $De$ with that of Leal (Reference Leal1975).

4.6. Comparison of second-order fluid result with Leal (Reference Leal1975)

Leal (Reference Leal1975) considered the motion of a fibre in a second-order fluid and found

(4.48)$$\begin{gather} \dot{p}_2^{Leal}=- p_1p_2 (\,p_2+Vp_1 (1-2p_2^2)) \end{gather}$$

(4.49)$$\begin{gather}\dot{p}_3^{Leal}=- p_1 p_2 p_3(1-2V p_1p_3), \end{gather}$$

where $V=-({3 \lambda \gamma }/{16 \log (\kappa )})M_1(1+2\epsilon _1)$, $\lambda =\varPhi _3 U/\mu l$ is the polymer relaxation time ($\mu$ is the zero-shear-rate viscosity, $l$ the particle half-length and $U$ a characteristic velocity scale) and $\epsilon _1=\varPhi _2/\varPhi _3$,

(4.50a,b)\begin{equation} \varPhi_2=-\lim_{\gamma\rightarrow 0}\frac{\sigma_{11}-\sigma_{22}}{2\gamma^2} \quad\text{and}\quad \varPhi_3=\lim_{\gamma\rightarrow 0}\frac{\sigma_{11}-\sigma_{33}}{\gamma^2}, \end{equation}

where $\sigma _{11}$, $\sigma _{22}$, $\sigma _{33}$ and $\gamma$ are the first, second and third normal stresses and the shear rate. and $M_1$ is a positive number depending upon the shape of the particle. In non-dimensional terms $\lambda =2De$, $V=({-3De}/{8\log (\kappa )})M_1(1+2\epsilon _1)$. Here $(1+2\epsilon _1)$, and hence $V$ in Leal's theory, is proportional to the second normal stress difference. Up to $O(De p_2)$ Leal's results are, therefore,

(4.51)$$\begin{gather} \dot{p}_2^{Leal}=-p_1p_2^2+\frac{3De}{8\log(\kappa)}M_1 (1+2\epsilon_1)p_1^2p_2 +O(De^2,p_2^2De), \end{gather}$$

(4.52)$$\begin{gather}\dot{p}_3^{Leal}=-p_1p_2p_3+O(De^2,p_2^2De). \end{gather}$$

By taking only up to $O(De)$ terms from (4.46), our results up to this order for a prolate spheroidal particle in a second-order fluid are

(4.53)$$\begin{gather} \dot{p}_2=-p_1p_2^2-\frac{c\,De}{2\log(2\kappa)-3}p_1^2p_2 +O(c\,De^2,cp_2^2De), \end{gather}$$

(4.54)$$\begin{gather}\dot{p}_3=-p_1p_2p_3+O(c\,De^2,cp_2^2De). \end{gather}$$

Hence, our theory's first viscoelastic effects at low $De$ have the same functional dependence of rotation rates on the particle orientation as that of Leal's second-order fluid theory. However, our theory predicts these effects to arise from the first normal stress difference. In contrast, Leal's theory predicts these to emerge from the second normal stress difference. According to our theory, a prolate spheroid rotating in a simple shear flow of a Boger fluid that has zero second normal stress difference (Magda et al. Reference Magda, Lou, Baek and DeVries1991), will exhibit a different rotational motion at a finite but small $c\, De$ as compared with $c=0$ (Newtonian fluid). Specifically, a particle spirals towards the stable limit cycle close to the vorticity, as discussed in the next section. However, Leal's theory predicts Jeffery rotations or no effect of viscoelasticity. Brunn (Reference Brunn1977) considered the motion of rigid particles in second-order fluid. While Brunn (Reference Brunn1977) does not calculate the rotation rates for a rod-like particle such as prolate spheroid, for a dumbbell representing two spheres joined by a thin, rigid rod, he also obtains the same functional dependence of rotations rates on the particle orientation. However, instead of the factor $M_1(1+2\epsilon _1)$, Brunn (Reference Brunn1977) obtains $(1+4\epsilon _1)$ so that the rotation rate of a dumbbell in a second-order fluid with zero second normal stress difference is finite.

5.1. Low dimensional orbits: log-rolling and tumbling

We begin our analysis with two convenient orientational states: (a) log-rolling with the particle aligned with the vorticity axis, and, (b) tumbling in the FGP ($p_3=0$). Due to symmetry, and as suggested by (4.46), viscoelasticity does not change the particle's orientation from the log-rolling state. Also, in the log-rolling state the theory predicts that the polymer-induced solvent and elastic torques due to viscoelasticity, i.e. ${{G}_3^{PIST}}^{(1)}$ from (4.30) and ${{G}_3^{Elastic}}^{(1)}$ from (4.40), are both zero. Thus, the polymers do not change the log-rolling angular velocity of the particle at $O (c)$.

In the FGP, close to the flow direction $p_1\approx 1$, the equation governing the particle orientation is

(5.2)\begin{equation} \dot{p}_2\approx-p_2^2-\frac{1}{\kappa^2}-c\,De\frac{p_2}{2\log(2\kappa)-3}.\end{equation}

Therefore, we notice that viscoelasticity ($c\, De>0$) reduces the rotation rate of the particle as compared with the Newtonian value ($c\, De=0$). This equation has an analytical solution

(5.3)\begin{equation} p_2=-b\tan({(t-t_0)b})\frac{c\,De}{4\log(2\kappa)-6},\end{equation}

where

(5.4)\begin{equation} b^2={\frac{1}{\kappa^2}-\left(\frac{c\,De}{4\log(2\kappa)-6}\right)^2}.\end{equation}

At small values of $c\, De$ when $b^2>0$, the solution is periodic with a time-period $2{\rm \pi} /b$. Upon increasing $c\, De$, $b=0$ when $c\,De_{crit}=(4\log (2\kappa )-6)/\kappa$ and an infinite period bifurcation occurs. Further increasing $c\, De$ lead to $b^2<0$ and two fixed points appear at $[\,p_{2,{flow}}^{0-},0]$ and $[\,p_{2,{flow}}^{0--},0]$ within the FGP, where

(5.5a,b)\begin{equation} p_{2,{flow}}^{0-}=-\frac{c\,De}{4\log(2\kappa)-6}+\sqrt{-b^2},\quad p_{2,{flow}}^{0--}=-\frac{c\,De}{4\log(2\kappa)-6}-\sqrt{-b^2}.\end{equation}

There are two additional fixed points at $-p_{2,{flow}}^{0-}$ and $-p_{2,{flow}}^{0--}$ that can be mathematically obtained by repeating the above analysis near $p_1=-1$. In the sense of Jeffery orbits, the fixed points are downstream of the flow–vorticity plane and very close to the flow direction when $\kappa$ is large, and $c\, De$ is small; $p_{2,{flow}}^{0--}$ is farther downstream. Within the FGP, particle trajectories approach $\pm p_{2,{flow}}^{0-}$, while they depart $\pm p_{2,{flow}}^{0--}$. Therefore, in the $b^2<0$ regime, a particle placed in the FGP approaches a steady state orientation $p_2=\pm p_{2,{flow}}^{0-}$. We will later observe in § 5.3 that in the orientation space near the flow-direction off the FGP, i.e. for a finite $0< p_3\ll 1$, the trajectories may either approach or leave these two fixed points along the vorticity direction.

5.1.1. More accurate location of fixed points on FGP in $b^2<0$ regime

Above, we analysed the equations near the flow direction under the assumption $p_1\approx 1$ and found two fixed points in the $b^2<0$ regime at orientations given by (5.5a,b). These expressions are only valid when the fixed points are near the flow direction. However, as $c\, De$ is increased at a given $\kappa$, the fixed points separate from each other in the flow gradient plane. The fixed point at $p_{2,{flow}}^-$ moves closer to the flow direction, while $p_{2,{flow}}^{--}$ moves away from the flow direction. To better estimate the latter's location we relax the assumption of $p_1\approx 1$ in the expression for ${\textrm {d}p_2}/{\textrm {d}t}$ and obtain the improved expressions

(5.6) \begin{align} \left.\begin{gathered} \tilde{p}_{2,{flow}}^{0-}=-\sqrt{\frac{-f^2/\kappa^2+{c^2De^2}/2-fc \, De\sqrt{-\tilde{b}^2}}{{c^2De^2+f^2}}},\\ \tilde{p}_{2,{flow}}^{0--}=-\sqrt{\frac{-f^2/\kappa^2+{c^2De^2}/2+fc\,De\sqrt{-\tilde{b}^2}}{{c^2De^2+f^2}}},\end{gathered}\right\} \end{align}

where $\tilde {b}^2={{(1+\kappa ^2)}/{\kappa ^4}-({c\,De}/{(2f)})^2}\approx b^2$ and $f=2\log (2\kappa )-3$. From figure 3, as $c\, De$ or $\kappa$ are varied, we observe the qualitative behaviours of $\tilde {p}_{2,{flow}}^{0-}$ and $\tilde {p}_{2,{flow}}^{0--}$ are the same as those of ${p}_{2,{flow}}^{0-}$ and ${p}_{2,{flow}}^{0--}$ mentioned above. Here $\tilde {p}_{2,{flow}}^{0-}$ is shown with dashed lines and $\tilde {p}_{2,{flow}}^{0--}$ with solid lines in the region $\widetilde {b^2}\approx {b^2}<0$. At the beginning of the dashed and solid lines, $\widetilde {b^2}\approx {b^2}=0$, for a given $\kappa$, $\tilde {p}_{2,{flow}}^{0-}=\tilde {p}_{2,{flow}}^{0--}$, because these fixed points arise out of an infinite period bifurcation (Strogatz Reference Strogatz2018), as mentioned earlier.

Figure 3. Variation of the fixed points’ location on the FGP with $c\, De$ for various $\kappa$ in the $\widetilde {b^2}\approx {b^2}<0$ regime.

Therefore, while polymers have no influence on the log-rolling motion, they slow down or stop the tumbling motion of a prolate spheroidal particle within the FGP. We now move on to consider three-dimensional orbits by first observing the particle's orientation trajectories close to the vorticity axis in § 5.2 and then the trajectories near the FGP in § 5.3.

5.2. Effect of viscoelasticity near the vorticity direction, $p_3\approx 1$

Equation (4.46) has a fixed point on the vorticity axis,

(5.7)\begin{equation} \boldsymbol{p}^0_{vort}=\left[0\quad 0\quad 1 \right]. \end{equation}

Due to the constraint $||\boldsymbol {p}||_2=1$, we consider the linear stability of the $p_1-p_2$ dynamical system using $p_3=\sqrt {1-p_1^2-p_2^2}$. In the $p_1-p_2$ coordinate system, at the fixed point $[0\ 0]$, the eigenvalues are

(5.8a,b)\begin{align} \lambda_{1,{vort}}=2c\,De/\kappa^2-1/\kappa\sqrt{4c^2De^2/\kappa^2-1},\quad \lambda_{2,{vort}}=2c\,De/\kappa^2+1/\kappa\sqrt{4c^2De^2/\kappa^2-1},\end{align}

with the corresponding eigenvectors

(5.9a,b)\begin{equation} \boldsymbol{v}_{1,{vort}}=[ \lambda_{1,{vort}}\quad 1 ],\quad \boldsymbol{v}_{2,{vort}}=[ 1\quad \lambda_{2,{vort}} ]. \end{equation}

The fixed point at the vorticity axis undergoes a Hopf bifurcation at $c\, De=0$ and is unstable for all finite $c\, De$. The vorticity axis is a centre (imaginary eigenvalues) when $c\, De=0$. In the presence of viscoelasticity, $c\, De\ne 0$, it becomes a hyperbolic fixed point (eigenvalues with non-zero real part). Therefore, the linear stability analysis provides qualitative insight into the effect of viscoelasticity on the behaviour of fibre orientation governed by the full nonlinear equation (4.46) near the vorticity axis when $c\, De\ne 0$. This follows from the Hartman–Grobman theorem (see Guckenheimer & Holmes Reference Guckenheimer and Holmes2013), which guarantees the homeomorphism between the linearized and full nonlinear system near hyperbolic fixed points, while also preserving the time parametrization. The instability of the vorticity axis arises from the polymer conformation driven by the force doublet and $O(1/\kappa ^2)$ Stokeslet discussed in § 4.1. When $0< c\,De<\kappa /2$ a small perturbation leads to a particle departing the vorticity axis in a spiral fashion (as the vorticity fixed point is an unstable spiral). However, for $c\,De>\kappa /2$ a particle departs the vorticity axis monotonically (as the vorticity fixed point is an unstable fixed point).

The linear stability analysis is confirmed by the full numerical integration of (4.46) for $\kappa =50$ at $c\, De=24$ (also for $c\, De=5$ and 6) versus $c\, De=26$ in figure 4(a) zoomed near the vorticity axis. At $\kappa =50$, $c\, De=26$ is a point in $R_{flow}^{(2)}$, $c\, De=24$ and 6 are points in $R_{flow}^{(1)}$ and $c\, De=5$ is a point in $R_{flow-vort}^{(3)}$. In $R_{flow}^{(2)}$, the particle drifts out of the vorticity axis monotonically, and in the rest, it spirals out. We will discuss the dynamical system's features of figure 4(b) in more detail later in § 5.4.1 where we will find the unstable spiral at the vorticity axes to be surrounded by a stable limit cycle up to a cut-off $c\, De$. Hence, the Hopf bifurcation occurring at $c\, De=0$ at the vorticity axis is supercritical (Strogatz Reference Strogatz2018). In $R_{flow-vort}^{(3)}$ ($c\, De=5$), the particle starting very close to the vorticity axis spirals outwards into the stable limit cycle, where it then moves periodically around the vorticity axis. The limit cycle does not exist in $R_{flow}^{(1)}$ and $R_{flow}^{(2)}$, so that, after the particle either spirals or monotonically drifts out of the vorticity axis, it ends up in the flow aligned state as it approaches the stable fixed point ($\pm \tilde {p}_{2,{flow}}^{0-}$) near the flow direction in the FGP for $c\, De=6$, 24 and 25.

Figure 4. Various orientation behaviours near the vorticity axis: trajectories of particle orientation starting very close to the vorticity axis at different $c\, De$ in $R_{flow}^{(1)}$ ($c\, De=6$, 24), $R_{flow}^{(2)}$ ($c\, De=26$) and $R_{flow-vort}^{(3)}$ ($c\, De=5$) at $\kappa =50$. All trajectories start at the same point. Panel (a) is the same as panel (b) (showing complete particle trajectory) but zoomed near the vorticity axis ($p_3=1$). The grey surface is the unit sphere, i.e. the orientation space.

The identification of the stable limit cycle near the vorticity axis is made possible by including the flow generated by force dipoles per unit length (flow from Cox (Reference Cox1971)) along the fibre in addition to the flow generated by the force per unit length (flow from Batchelor (Reference Batchelor1970)). If this dipole-generated flow is neglected from our theory the stable limit cycle will not be predicted and instead when the theory predicts particle to be repelled from the FGP it will approach the vorticity axis. A previous slender body theory for a second-order fluid by Férec et al. (Reference Férec, Bertevas, Khoo, Ausias and Phan-Thien2021), relying on only the flow generated by the force per unit length, predicts the fibre to align with the vorticity axis. When the axis of symmetry of an axisymmetric slender fibre is aligned in the flow–vorticity plane of the imposed flow no force per unit length is exerted by the fluid (Newtonian or polymeric). This is because in the flow–vorticity-aligned state the imposed flow has no variation along the fibre axis. However, according to (29) and (35) of Férec et al. (Reference Férec, Bertevas, Khoo, Ausias and Phan-Thien2021) when a fibre is in the flow–vorticity plane, the force per unit length is non-zero and proportional to the polymer relaxation time. This implies an error in the expression for the tension force in (29) of Férec et al. (Reference Férec, Bertevas, Khoo, Ausias and Phan-Thien2021).

Figure 3(a) of Wang, Yu & Lin (Reference Wang, Yu and Lin2019) clearly shows that a $\kappa =4.0$ spheroid started away from the vorticity axis in a plane Couette flow at $De=0.1$ reaches a closed orbit around the vorticity axis, instead of approaching the axis. Figures 2(a) and 2(b) of d'Avino et al. (Reference d'Avino, Hulsen, Greco and Maffettone2014) at $De=1.0$ and 2.0, respectively, for a $\kappa =4.0$ spheroid in an unbounded simple shear flow show drift towards the vorticity axis, but do not show the particle approaching the axis. Although these studies were conducted at high polymer concentration, $c$, and small $\kappa$, outside the formal range of validity of our theory, they indicate the possibility of a stable limit cycle around the vorticity axis. Furthermore, in the second-order fluid regime with $De=0.1$, careful observation of figure 8(a) of the theoretical investigation of Wang, Tai & Narsimhan (Reference Wang, Tai and Narsimhan2020) shows a $\kappa =3.0$ spheroid approaching a stable limit cycle, instead of reaching the vorticity axis of an unbounded parabolic slit flow ($u=1-y^2$). Lastly, while the boundary element formulation aided study of Phan-Thien & Fan (Reference Phan-Thien and Fan2002) in an Oldroyd-B fluid shows the $\kappa =2.0$ particle drifting towards the vorticity axis, the simulation (see figures 7 and 8 of their paper) stops before we can conclude if it will approach the axis or a stable limit cycle.

5.3. Effect of viscoelasticity near the FGP and flow direction

We described the effect of viscoelasticity on the particle's motion near the flow direction within the FGP in § 5.1. Here, we consider the particle motion near (but not exactly on) the FGP. The analysis of § 5.1 for the $p_2$ (gradient) direction is valid even outside the FGP. Near the flow direction, $p_1\approx 1$, the orientation dynamics and solution in the $p_2$ direction are given by (5.2) and (5.3). The orientation dynamics in the $p_3$ (vorticity) direction for $p_1\approx 1$ is governed by the simplified equation (from (4.46))

(5.10)\begin{equation} \dot{p}_3\approx-p_2p_3+c\,Dep_3\left(-\frac{4}{\kappa^2}+\frac{p_2}{8\log(2\kappa)-12}\frac{15{\rm \pi}^2De}{3{\rm \pi}^2+5De}\right)+O(\,p_3^2).\end{equation}

Its closed form solution is

(5.11)\begin{equation} p_3=C\exp(c\,De\zeta (t-t_0))\cos({(t-t_0)b})^{-1+c\,De\gamma/(8\log(2\kappa)-12)}\end{equation}

where $b^2$ is given in (5.4), and

(5.12a,b)\begin{equation} \gamma=\frac{15{\rm \pi}^2De}{3{\rm \pi}^2+5De},\quad \zeta=\frac{1}{4\log(2\kappa)-6}-c\,De\frac{\gamma}{2(4\log(2\kappa)-6)^2}-\frac{4}{\kappa^2}.\end{equation}

When $b$ is real-valued, i.e. $b^2>0$ (5.4), there are no fixed points on the FGP. Thus, the particle undergoes periodic motion in the FGP, as discussed in § 5.1. For a particle perturbed slightly away from the FGP, from (5.11) we note that the sign of $\zeta$ determines whether the particle will be attracted to or repelled from the FGP ($p_3=0$). Here $\zeta >0$ represents the region $R_{vort}^{(1)}$ where the particle spirals away from the FGP and $\zeta <0$ represents the region $R_{flow-vort}^{(1)}$ where the particle spirals into the FGP and continues tumbling within the plane (albeit with a larger orbit period $2{\rm \pi} /b$ than in Newtonian case). In the $b^2>0$ regime, the phase portrait of the dynamical system of (4.46) (through the simplified equations (5.2) and (5.10)) projected in the $p_2\unicode{x2013} p_3$ plane near $p_1=1$ is shown in figure 5 for $\zeta >0$ and $\zeta <0$. The phase flow in the $p_2\unicode{x2013} p_3$ plane in the case of a Newtonian fluid ($c\, De=0$) is similar to that in the region $R_{vort}^{(1)}$ (figure 5a), but it is symmetric about $p_2=0$ (not shown). Therefore, a particle in a Newtonian fluid continues in a particular (periodic) Jeffery orbit before and after $p_2=0$. At small but finite $c\, De$, i.e. in the region $R_{vort}^{(1)}$, represented figure 5(a), this symmetry about $p_2=0$ is broken and a particle comes out of the $p_2=0$ plane with a larger $\dot {p}_3$ velocity than it enters the plane. This explains the drift towards the vorticity axis (greater $p_3$). In the region $R_{flow-vort}^{(1)}$ represented by figure 5(b) we can observe that the phase flow points towards $p_3=0$, indicating migration of a particle towards the FGP. The FGP is a stable limit cycle in $R_{flow-vort}^{(1)}$ and an unstable limit cycle in $R_{vort}^{(1)}$.

Figure 5. Phase portraits of (5.2) and (5.10) in the $b^2>0$ regime in the gradient ($p_2$)–vorticity ($p_3$) plane with $\zeta >0$ (a) and $\zeta <0$ (b). When $b^2>0$ and $\zeta >0$, i.e. in (a) representing $R_{vort}^{(1)}$ close to the flow gradient plane ($p_3=0$), the $p_3$ component of the phase velocity changes sign from negative to positive along with an increase in magnitude downstream of the $p_2=0$ plane indicating a departure of a particle from the flow gradient plane at a rate higher than it approaches the plane. However, when $b^2>0$ and $\zeta <0$, i.e. in (b) representing $R_{flow-vort}^{(1)}$ the $p_3$ component of the phase velocity is negative for all $p_2$ indicating an approach towards the flow gradient plane. Since $p_2$ never approaches zero on the flow gradient plane for $b^2>0$ it is an unstable and stable limit cycle for $R_{vort}^{(1)}$ ($\zeta >0$, a) and $R_{flow-vort}^{(1)}$ ($\zeta <0$, b), respectively.

The numerically integrated trajectories for parameters chosen in $R_{vort}^{(1)}$ and $R_{flow-vort}^{(1)}$ for $c=0.005$ are shown in figures 6 and 7, respectively. The prediction of the spiral exit of the orientation trajectories from the FGP in $R_{vort}^{(1)}$ and the spiral approach towards it in $R_{flow-vort}^{(1)}$ mentioned above is confirmed from these figures. The spiralling rate in $R_{vort}^{(1)}$ increases with $c\, De$ and $\kappa$ (not shown).

Figure 6. In $R_{vort}^{(1)}$ (shown here for $c=0.005$, $\kappa =50, c\, De=0.01$), trajectories starting near the FGP (exemplified here with the blue trajectory) spiral out of the plane. Globally they approach the same stable limit cycle as the trajectories starting near the vorticity axis (exemplified here with the orange trajectory). The blue trajectory starting near the flow direction spans a larger portion of phase space in $p_2$, but we show the region close to the flow–vorticity plane to highlight the limit cycle.

Figure 7. In $R_{flow-vort}^{(1)}$ (shown here for $c=0.005$, $\kappa =10, c\, De=0.48$), trajectories of particle orientation starting near the FGP (blue) spiral into the plane. Globally they emanate from an unstable limit cycle – the boundary between blue and green trajectories (that are started very close to each other in this numerical integration). There is a stable limit cycle above this unstable limit cycle at the boundary between green and orange trajectories.

As $b^2$ is reduced by increasing $c\, De$ or $\kappa$, the time period, $2{\rm \pi} /b$, increases. The effect of viscoelasticity (increasing $c\, De$) to increase the time period is consistent with the experimental observations of Gauthier et al. (Reference Gauthier, Goldsmith and Mason1971) and Bartram et al. (Reference Bartram, Goldsmith and Mason1975). As mentioned earlier in § 5.1, when $b^2=0$, due to an infinite period bifurcation (Strogatz Reference Strogatz2018) two fixed points emerge on the FGP, and we discuss the trajectories off the FGP in the $b^2<0$ regime next.

5.3.1. Monotonic behaviour near FGP: ${b^2<0}$

Equation (5.2) has two fixed points at $p_{2,{flow}}^{0-}$ and $p_{2,{flow}}^{0--}$ from (5.5a,b) (or their more accurate values in (5.6)). For the dynamical system in the $p_2\unicode{x2013} p_3$ plane near the flow direction, defined by the system of (5.2) and (5.10), these fixed points are $[\,p_{2,{flow}}^{0-},0]$ and $[\,p_{2,{flow}}^{0--},0]$. The eigenvectors at each of the fixed points are the same, i.e.

(5.13a,b)\begin{equation} \boldsymbol{v}_{1,{flow}}=\left[1\quad 0\right],\quad \boldsymbol{v}_{2,{flow}}=\left[0\quad 1\right]. \end{equation}

The corresponding eigenvalues at $[\,p_{2,{flow}}^{0-},0]$ and $[\,p_{2,{flow}}^{0--},0]$ are $\lambda _{i,{flow}}^{0-}$ and $\lambda _{i,{flow}}^{0--}$, where $i\in [1,2]$,

(5.14)\begin{equation} \left.\begin{gathered} \lambda_{1,{flow}}^{0-}=-2\sqrt{-b^2},\quad \lambda_{2,{flow}}^{0-}=\zeta c\, De -(\zeta+4/\kappa^2) (4\log(2\kappa)-6)\sqrt{-b^2}, \\ \lambda_{1,{flow}}^{0--}=2\sqrt{-b^2},\quad \lambda_{2,{flow}}^{0--}=\zeta c\, De +(\zeta+4/\kappa^2) (4\log(2\kappa)-6)\sqrt{-b^2}. \end{gathered}\right\}\end{equation}

Both the fixed points are hyperbolic in the $b^2<0$ regime. Thus, similar to the fixed point at the vorticity axis using the Hartman–Grobman theorem, the linear stability of the fixed points allows qualitative insight into the complete nonlinear behaviour of particle orientation. Both eigenvalues of both the fixed points are real in the $b^2<0$ regime. The first eigenvalues ($\lambda _{1,{flow}}^{0-}$ and $\lambda _{1,{flow}}^{0--}$) of both the fixed points do not change sign in the $b^2<0$ regime, i.e. $\lambda _{1,{flow}}^{0-}<0$ and $\lambda _{1,{flow}}^{0--}>0$. But, the second eigenvalues ($\lambda _{1,{flow}}^{0-}$ and $\lambda _{1,{flow}}^{0--}$), that are initially positive in the $b^2<0$ regime, become negative at different parameter ($c$, $De$ and $\kappa$) values. Here $\lambda _{2,{flow}}^{0-}$ becomes negative first. $\lambda _{2,{flow}}^{0-}=0$ is the bifurcation boundary where $[\,p_{2,{flow}}^{0-},0]$ changes from a saddle point (stable manifold along gradient direction and unstable along vorticity) to a stable fixed point. Here $\lambda _{2,{flow}}^{0--}=0$ is the bifurcation boundary where $[\,p_{2,{flow}}^{0--},0]$ changes from an unstable node to a saddle point (stable manifold along vorticity direction and unstable along gradient). Therefore, three new regimes arise within the $b^2<0$ region that are labelled in figure 2 as (1) $R_{vort}^{(2)}$ ($\lambda _{1,{flow}}^{0-}>0$ and $\lambda _{1,{flow}}^{0--}>0$), (2) $R_{flow-vort}^{(2)}$ ($\lambda _{1,{flow}}^{0-}<0$ and $\lambda _{1,{flow}}^{0--}>0$) and (3) $R_{flow-vort}^{(3)}+R_{flow}$ ($\lambda _{1,{flow}}^{0-}<0$ and $\lambda _{1,{flow}}^{0--}<0$). The phase flow close to the fixed points in the gradient–vorticity plane for these three cases is shown in figure 8. We can observe the phase flow approaching $[\,p_{2,{flow}}^{0-},0]$ and leaving $[\,p_{2,{flow}}^{0--},0]$ along the gradient direction ($p_2$) for all three cases. This implies that a particle with an initial orientation close to the FGP will approach the flow direction (specifically the fixed point $[\,p_{2,{flow}}^{0-},0]$) while slowing down. For the parameters within $R_{vort}^{(2)}$, the particle will then monotonically drift away from the FGP. This is similar to the experimental observation of Bartram et al. (Reference Bartram, Goldsmith and Mason1975) discussed in § 1. For a given $\kappa$ and $c$, at larger $c\, De$ than $R_{vort}^{(2)}$, within the regions $R_{flow-vort}^{(2)}$ and $R_{flow-vort}^{(3)}$, the particle starting near the FGP (and for all starting orientations in $R_{flow}$) achieves a stable orientation near the flow direction, similar to the experiments of Iso et al. (Reference Iso, Cohen and Koch1996a).

Figure 8. Phase portraits of the system of (5.2) and (5.10) in the $b^2<0$ regime. In this regime, two fixed points exist on the FGP ($p_3=0$) close to the flow direction. Both the fixed points are downstream of the flow–vorticity plane ($p_2=0$). An unstable (red marker) and a saddle (green marker) node with its unstable manifold along the $p_3$ axis in (a) ($R_{vort}^{(2)}$) indicates a monotonic drift of the particle away from the flow gradient plane. A stable fixed point (blue marker) near the flow direction ($p_2\approx 0, p_3=0$) in (b) ($R_{flow-vort}^{(2)}$) and (c) ($R_{flow-vort}^{(3)}$) panels indicate that particles with starting orientation near the FGP may align near the flow direction. The presence of an unstable point in $R_{flow-vort}^{(2)}$ (b) instead of a saddle point with its stable manifold perpendicular to the flow gradient plane in $R_{flow-vort}^{(2)}$ (c) in addition to the stable point in these cases indicates a lower proportion of trajectories leading to flow alignment in $R_{flow-vort}^{(2)}$ than in $R_{flow-vort}^{(3)}$.

The locations of fixed points corresponding to $\tilde {p}_{2,{flow}}^{0-}$ and $\tilde {p}_{2,{flow}}^{0--}$ in the orientation space are $\left [ \pm \sqrt {1-(\tilde {p}_{2,{flow}}^{0-})^2}\ \pm \tilde {p}_{2,{flow}}^{0-} \ 0 \right ]$ and $\left [\pm \sqrt {1-(\tilde {p}_{2,{flow}}^{0--})^2} \ \pm \tilde {p}_{2,{flow}}^{0--} \ 0 \right ]$. These locations match the corresponding fixed points in the numerically integrated trajectories in $R_{vort}^{(2)}$, $R_{flow-vort}^{(2)}$, and, $R_{flow-vort}^{(3)}+R_{flow}$ shown in the plots of figures 9, 10 and 11, respectively. We mark the analytical locations of these fixed points with two different coloured markers on the FGP ($p_3=0$) in figures 9, 10 and 11 and show the nearby trajectories to approach/leave these locations in the fashion described by the linear stability theory discussed above.

Figure 9. In $R_{vort}^{(2)}$ (shown here for $c=0.005$, $\kappa =50, c\, De=0.3$), a particle's orientation trajectories approach the stable limit cycle around and near the vorticity axis. In contrast to $R_{vort}^{(1)}$ (figure 6) where the trajectories spiral away from the FGP, here they leave the FGP monotonically. We show the region close to the flow–vorticity plane as this is where the different attractors lie. The grey surface is the unit sphere, i.e. the orientation space.

Figure 10. In $R_{flow-vort}^{(2)}$ (shown here for $c=0.005$, $\kappa =100$ and $c\, De=0.9$) a particle's orientation trajectories that start close to the FGP (solid lines) approach the flow direction either on the same or the opposite side of the gradient–vorticity plane. Trajectories farther away from the FGP (dashed lines) approach the stable limit cycle near the vorticity axis. Due to a saddle and stable node close to each other on the FGP and a stable limit cycle near the vorticity axis, another saddle node emerges on the faster eigendirection of the stable node. The grey surface is the unit sphere, i.e. the orientation space.

Figure 11. In $R_{flow-vort}^{(3)}$ (shown here for $c=0.005$, $\kappa =100$ and $c\, De=1.2$) the behaviour of a particle's orientation trajectories is similar to that in $R_{flow-vort}^{(2)}$ shown in figure 10. The primary difference between $R_{flow-vort}^{(2)}$ and $R_{flow-vort}^{(3)}$ is that, in the former, the fixed point on the FGP farther from the flow direction is an unstable fixed point (figure 10), while it is a saddle node in the latter (shown here). Additionally, in $R_{flow-vort}^{(3)}$, there is an unstable node at the intersection of the stable manifold of the saddle points in FGP and flow–vorticity plane. Trajectories with solid lines end up at one of the fixed points near the flow direction, and those with dashed lines end in the stable limit cycle near the vorticity axis.

Numerical integration of the governing equation shown in figure 9 confirms the existence of the stable fixed point and unstable node predicted by linear stability analysis. In agreement with the $p_2-p_3$ phase plane analysis above, the particle starting near the flow gradient plane in $R_{vort}^{(2)}$ approaches the flow direction. It then monotonically departs the FGP along the flow–vorticity plane. The primary difference between $R_{vort}^{(1)}$ and $R_{vort}^{(2)}$ is that the particle leaves the FGP spirally in the former (figure 6) and monotonically in the latter (figure 9). As mentioned earlier, these trajectories are reminiscent of that in the experimental observations of Bartram et al. (Reference Bartram, Goldsmith and Mason1975). We can only make qualitative comparisons with Bartram et al. (Reference Bartram, Goldsmith and Mason1975) due to unknown non-Newtonian properties of their viscoelastic fluid. Furthermore, the velocity field around the $\kappa =9.1$ particles used in their experiments is likely to have a quantitative difference from the SBT approximations used in our theory.

In $R_{flow-vort}^{(2)}$, the numerically integrated trajectories shown in figure 10 reveal the presence of stable and unstable nodes in the FGP at the locations predicted by the linear stability analysis. Similarly, the existence of a stable node close to the flow direction and a saddle point farther away on the FGP is confirmed from the numerical integration shown in figure 11 for the parameters chosen in $R_{flow-vort}^{(3)}$. In $R_{flow-vort}^{(2)}$ and $R_{flow-vort}^{(3)}$ the trajectories starting close to the FGP approach a stable orientation near the flow direction (stable fixed point at $\tilde {p}_{2,{flow}}^{0-}$).

5.4. Global particle orientation dynamics

In § 5.1 we analysed the equations for the rotational motion of a particle at the vorticity axis and within the FGP. We found that viscoelasticity does not lead to a non-zero $\dot {\boldsymbol {p}}$ at the vorticity axis and it does not change the log-rolling rotational velocity. In the FGP, we observed that a small amount of viscoelasticity (small $c\, De$) leads to a larger orbit time period. Further increasing $c\, De$ leads to an infinite period bifurcation and birth of two fixed points, thereby leading to particle migration towards the flow direction when in the FGP. In § 5.2, by linearizing the governing equation at the vorticity axis, we showed that the vorticity axis is an unstable spiral for $c\, De<\kappa$. For $c\, De>\kappa$, it is an unstable node. Therefore, a small perturbation of the particle orientation from the vorticity axis leads to departure from the axis. This was corroborated by the orientation trajectories obtained by numerical integration of the system of governing equations. By analysing the system and its fixed points (when they exist) near the flow direction in the FGP and performing numerical integration of the trajectories starting near but not on the FGP in § 5.3, we determined the orientational dynamics in the region nearby the flow gradient plane. The particle either spirals away from or towards the FGP when there are no fixed points in FGP. When the fixed points in FGP arise, the one closer to the flow direction is either a saddle point (with its stable manifold in the FGP) or a stable node. Therefore, the particle first migrates towards the flow direction nearly parallel to the FGP. Then it may either settle near the flow direction or monotonically drift away from the FGP along the flow–vorticity plane. The boundaries in the $c\, De-\kappa$ space separating these different qualitative behaviours were shown in figure 2 for $c=0.005$. The qualitative nature of these boundaries is not sensitive to the exact value of $c$. As we observed in figures 6, 7, 9, 10 and 11, there are other invariant features (stable limit cycle, unstable limit cycle, saddle point, and unstable node) that arise in the orientation space off the FGP and the vorticity axis. These features may be viewed as a result of the orientation space being restricted to a unit sphere.

In $R_{vort}^{(1)}$ (figure 2), which extends to arbitrarily small (but finite) viscoelasticity (small $c\, De$), the vorticity axis is a spiral source, and the FGP is an unstable limit cycle. Therefore, at least one stable attractor must exist between these two regions on the unit sphere. In the case of a Newtonian fluid, the particle follows an initial-condition-dependent Jeffery orbit, i.e. one of a concatenation of neutrally stable periodic orbits centred at the vorticity axis. Therefore, at small (but finite) $c\, De$ the simplest change in phase-plane dynamics leads to a stable limit cycle around and near one of these limit cycles, as shown in figure 6. As shown by linear stability analysis at the vorticity direction in § 5.2 the rate of deviation of trajectories away from the vorticity axis is driven by the $O(1/\kappa ^{2})$ flow generated by the force doublet and $O(1/\kappa ^{2})$ Stokeslet. Thus, a slender particle starting at any orientation between the vorticity axis and the FGP leads to a final orientation behaviour where the particle undergoes periodic motion very close to and around the vorticity axis.

In $R_{flow-vort}^{(1)}$ (figure 2), the FGP changes from being an unstable to a stable limit cycle. Therefore, an unstable invariant object or a repeller in the form of an unstable limit cycle exists between the stable limit cycle near the vorticity direction and the FGP, as shown in figure 7. Thus, in $R_{flow-vort}^{(1)}$, a particle with an initial orientation between the unstable limit cycle and FGP spirals to the FGP, where it undergoes a tumbling motion. A particle with an initial orientation between the unstable limit cycle above the FGP and the stable limit cycle near the vorticity axis spirals towards the stable limit cycle near the vorticity axis. A particle released at an arbitrarily small (but finite) angle from the vorticity axis spirals outwards. In both these cases, eventually, the particle undergoes perpetual periodic motion on the stable limit cycle close to the vorticity axis.

In $R_{vort}^{(2)}$ (figure 2), the FGP contains two fixed points as shown in figure 9. One is an unstable node, and the other is a saddle point with its unstable manifold perpendicular to the FGP. Thus, no other invariant feature (attractor or repeller) is needed between the FGP and the stable limit cycle near the vorticity axis, which is the only stable invariant object on the orientation space. Hence a particle starting at an arbitrarily small angle from the vorticity axis or the FGP ends up in a periodic motion very close to the vorticity axis.

Similar to $R_{vort}^{(1)}$, $R_{vort}^{(2)}$ and $R_{flow-vort}^{(1)}$ discussed above, there is a stable limit cycle around the vorticity axis in $R_{flow-vort}^{(2)}$ and $R_{flow-vort}^{(3)}$. Therefore, a particle with an initial orientation arbitrarily close to the vorticity direction (but not on the vorticity axis) eventually undergoes a periodic motion around the vorticity axis. This stable limit cycle exists for all regions except for $R_{flow}^{(1)}$ and $R_{flow}^{(2)}$ and it will be further analysed in § 5.4.1.

In $R_{flow-vort}^{(2)}$ (figure 2) the fixed point in the FGP closest to the flow direction is a stable node and the other fixed point in the FGP is an unstable node as shown in figure 10. Hence, a saddle node exists near the flow–vorticity plane to repel the particle orientation trajectories towards the stable limit cycle near the vorticity axis on one side and the stable fixed point near the flow direction on the FGP on the other. This saddle point receives trajectories emanating from the unstable point on the FGP. Thus, depending upon the initial orientation of the particle, it may eventually either obtain a stable orientation near the flow direction or undergo a periodic motion close to the vorticity axis.

The dynamics in $R_{flow-vort}^{(3)}$ (figure 2) are similar to those in $R_{flow-vort}^{(2)}$, but the unstable fixed point in the FGP in $R_{flow-vort}^{(2)}$ changes to a saddle node in $R_{flow-vort}^{(3)}$ with no qualitative change to the other invariant objects present in $R_{flow-vort}^{(2)}$. Therefore, to allow smooth and continuous phase flow, an unstable node occurs at the intersection of the saddle points in the FGP and near the flow–vorticity plane, as shown in figure 11. At a given $\kappa$, the value of $c\, De$ in $R_{flow-vort}^{(3)}$ is larger than that in $R_{flow-vort}^{(2)}$ (figure 2) and the separation between two fixed points increases with increasing $c\, De$ as shown by figure 3. Due to this increased separation and the unstable node above the flow gradient plane, more of the phase flow is directed towards the FGP in $R_{flow-vort}^{(3)}$ than in $R_{flow-vort}^{(2)}$. Thus, the basin of attraction for the stable fixed point near the flow direction is larger in $R_{flow-vort}^{(3)}$ as compared with $R_{flow-vort}^{(2)}$. In other words, more of the initial particle orientations lead to a stable final orientation near the flow direction in the $R_{flow-vort}^{(3)}$ than in $R_{flow-vort}^{(2)}$ as observed by comparing the proportion of solid lines in figures 10 and 11. When the stable limit cycle near the vorticity axis exists, the saddle point on the flow–vorticity plane is important in determining the proportions of initial orientations leading to the flow alignment and towards the periodic motion around the vorticity axis. Its location will be further analysed in § 5.4.2.

The $c\, De|_{cut\textrm -off}^{Limit}$ boundary in figure 2, beyond which the stable limit cycle near the vorticity direction does not exist, may be viewed as arising due to the interaction between the saddle node near the flow–vorticity plane and the stable limit cycle near the vorticity direction as these move in the orientation space upon increasing $c\, De$. We will discuss this in § 5.4.3.

5.4.1. Stable limit cycle around the vorticity axis

The previous study of Leal (Reference Leal1975) concerning a slender particle in simple shear flow of a second-order fluid predicts that viscoelasticity (although incorrect in attributing it to the second instead of the first normal stress difference) drives the particle orientation towards the vorticity axis. Leal (Reference Leal1975) only considered the Stokeslet fluid flow (of order $O(1/\log (\kappa ))$ and $O(1/\log (\kappa )^2)$) from the slender body theory. This velocity field is proportional to $p_2$ and hence has no effect when the particle is in the flow–vorticity plane. However, our study shows that a stable limit cycle exists around the vorticity axis due to the competition between polymer-driven torque arising from the doublet flow causing the particle to spiral away from the vorticity axis and that from the Stokeslet flow at $O(1/\log (\kappa ))$ causing it move towards the region around vorticity.

The real part of the eigenvalues of the fixed point at the vorticity axis in $R_{vort}^{(1)}$ is $2c\,De/\kappa ^2$ (5.8a,b), and the rate at which trajectories leave the FGP is $c\,De\zeta$ (5.11). Figure 12 shows the variation of $2c\,De/\kappa ^2$ and $c\,De\zeta$ with $\kappa$ for a few values of $c\, De$ at $c=0.005$. From this figure, we conclude that for $\kappa \gtrapprox 15$, trajectories leave the flow gradient plane much faster than the outward spiralling rate from the vorticity axis. Hence, the stable limit cycle shifts closer to the vorticity axis as $\kappa$ increases. Furthermore, as $\kappa$ increases beyond $\kappa \approx 20$, $c\,De\zeta$ remains almost constant while the rate of spiralling from the vorticity axis, $\kappa ^{-2}$, decreases. Therefore, the stable limit cycle approaches the vorticity axes at larger $\kappa$.

Figure 12. The variation with $\kappa$ of the rate of deviation of orientation trajectories away from the vorticity axis, $2c\,De/\kappa ^2$ (orange), and the flow gradient plane, $c\,De\zeta$ (black), in $R_{vort}^{(1)}$ ($b^2>0$) for a few values of $c\, De$ at $c=0.005$. On each curve, $b^2>0$ is to the left of the solid markers. An increasing gap between orange and black curves, with almost horizontal black curves, with $\kappa$ indicates a larger increase in the departure rate of the orientation trajectories from the flow gradient plane than that from the vorticity axis, implying a shift of the stable limit cycle closer to the vorticity axis.

We quantify the location of the limit cycle observed in the numerical integration of (4.46). For this purpose we transform this equation into the $C-\tau$ coordinate system (Leal & Hinch Reference Leal and Hinch1971), using

(5.15a,b)\begin{equation} C=\frac{\sqrt{p_2^2+p_1^2/\kappa^2}}{p_3}, \quad\tau=\frac{p_1}{p_2}, \end{equation}

to obtain

(5.16)$$\begin{gather} \frac{{\rm d} C}{{\rm d}t}=c\,De\tau^2\left(\frac{4(1+C^2)}{C\kappa^4}-\frac{C}{(4\log(2\kappa)-3)(\kappa^2+\tau^2)}\right)p_2^2+c\,DeO(\,p_2^4), \end{gather}$$

(5.17)$$\begin{gather}\frac{{\rm d}\tau}{{\rm d}t}=1+c\,De\tau p_2^2\left(\frac{\tau^2}{4\log(2\kappa)-3}+4\frac{\kappa^2+\tau^2}{C^2\kappa^4}\right)+O(c\,De). \end{gather}$$

In the Newtonian limit, ${\textrm {d}C}/{\textrm {d}t}=0$ and the particle trajectory is determined by its initial condition $C(t)=C(0)=C_N$. Each $C_N=[0,\infty ]$ represents a different Jeffery orbit. Here $C=0$ represents the log-rolling motion, i.e. the particle rotating about its major axis which is aligned with the vorticity axis, and $C=\infty$ represents the major axis rotating within the FGP. Here $\tau =t$ represents the phase on the Jeffery orbit. A periodic orbit can thus be ascertained using a single parameter $C$ that repeats periodically. In the range, $10\le \kappa \le 200$, for $c\, De\ge 0.01$ we evolve the system of equations for an initial condition starting near the vorticity direction, at $C=10^{-7}$. We evolve the equations for each choice of $c\, De$ and $\kappa$ until the change in the period-averaged $C$ across 100 successive Jeffery periods is less than $10^{-5}$. Finally, we define $\bar {C}_{Limit}^{numerical}$ as the average $C$ over the last 100 Jeffery periods. The $\bar {C}_{Limit}^{numerical}$ versus $c\, De$ for a few values of $\kappa \in [10,200]$ is shown in figure 13(a).

Figure 13. (a) The location of the stable limit cycle at different $c\, De$ and $\kappa$ at $c=0.005$ (nearly identical curves are obtained for other polymer concentrations in the range $0.01\le c\le 0.2$) indicated by the average $C=\sqrt {p_2^2+p_1^2/\kappa ^2}/{p_3}$, $\bar {C}_{Limit}^{numerical}$ on the limit cycle. Here $C=0$ is the vorticity axis and $C=\infty$ is the FGP. (b) The $\bar {C}_{Limit}^{numerical}$ on the limit cycle versus $\kappa$ at low $c\, De$.

Now, $\bar {C}_{Limit}^{numerical}$ is smaller than $0.1$ for $\kappa >10$, decreases with $\kappa$ and varies slightly with $c\, De$ at each $\kappa$. For $\kappa >50$ $\bar {C}_{Limit}^{numerical}$ is less than $0.01$ for a range of $c\, De$. Therefore, the limit cycle is close to the vorticity axis. Beyond a certain value of $c\, De=c\, De|_{cut\textrm -off}^{Limit}$, the limit cycle does not exist. This marks the boundary between $R_{flow-vort}^{(2)}$ and $R_{flow}^{(1)}$ for a given $\kappa$, as shown in figure 2. Now, $c\, De|_{cut\textrm -off}^{Limit}$ increases with $\kappa$ and for $10\le \kappa \le 200$ we find

(5.18)\begin{equation} c\, De|_{cut\text-off}^{Limit}\approx 0.1056\kappa-0.5183. \end{equation}

At larger $\kappa$, the stable limit cycle is closer to the vorticity axes and exists for a larger range of $c\, De$. Hence, marked as $R_{flow-vort}^{(1)}$ and $R_{flow-vort}^{(2)}$ in figure 2, there is a significant range of $c\, De$ at large $\kappa$ where two possible types of orientation dynamics are possible: a steady state orientation close to the flow direction and a periodic orbit near the vorticity direction. For $\kappa \lessapprox 20$, increasing $c\, De$ from very small values does not affect the position of the limit cycle much. However, as $c\, De$ approaches $c\, De|_{cut\textrm -off}^{Limit}$, the limit cycle moves away from the vorticity axis. For $\kappa \gtrapprox 20$, the limit cycle moves towards the vorticity axis upon increasing $c\, De$ from very small values. However, beyond a given $c\, De$, up to $c\, De|_{cut\textrm -off}^{Limit}$ the limit cycle's angular position is not influenced by $c\, De$. These features are reflected by the trends in $\bar {C}_{Limit}^{numerical}$ in figure 13(a). At low $c\, De$, the limit cycle is very close to one of the degenerate Jeffery orbits of the Newtonian case, i.e. $C$ does not change much on the limit cycle. The variation of $\bar {C}_{Limit}^{numerical}$ with $\kappa$ at low $c\, De=0.01, 0.02$ and 0.03 is shown in figure 13(b), where we can observe $\bar {C}_{Limit}^{numerical}$ to scale as $\kappa ^{-1.5}\log (2\kappa -3)$ for $10\le \kappa \le 100$ and approximately $\kappa ^{-1}$ for larger $\kappa$.

We find the particle orientation trajectories approaching this stable limit cycle close to the vorticity axis instead of spiralling towards the vorticity axis as concluded from the experiments of Gauthier et al. (Reference Gauthier, Goldsmith and Mason1971) and Iso et al. (Reference Iso, Koch and Cohen1996b). Gauthier et al. (Reference Gauthier, Goldsmith and Mason1971) claim the approach towards the vorticity axis by extrapolating the observed data. In some of the experimental observations of Iso et al. (Reference Iso, Koch and Cohen1996b) at small elasticity after undergoing initial spiralling, the particle oscillates around the vorticity axis without settling at a steady orientation. This may indicate the limit cycle discussed above. The theoretical rate of spiralling away from the vorticity axis ($2c\,De/\kappa ^2$) is small (figure 12), and small imperfections from the Couette cell might have led to secondary flows in the experiments. Thus, the existence and size of the stable limit cycle around the vorticity axis must be further tested in future experiments and numerical simulations.

5.4.2. Saddle node near the flow–vorticity plane in $R_{flow-vort}^{(2)}$, $R_{flow-vort}^{(3)}$, $R_{flow}^{(1)}$ and $R_{flow}^{(2)}$

From the trajectories of particle orientation in $R_{flow-vort}^{(2)}$ and $R_{flow-vort}^{(3)}$ shown in figures 10 and 11 we find the stable limit cycle around the vorticity axis and the stable node on the flow gradient plane near the flow axis to be separated by a saddle point. The saddle point arises on the boundary between $R_{vort}^{(2)}$ and $R_{flow-vort}^{(2)}$ when the saddle point near the flow direction on the FGP changes into a stable node. It also arises on the boundary between $R_{flow-vort}^{(3)}$ and $R_{flow-vort}^{(1)}$ when the unstable limit cycle in the FGP vanishes. The saddle point persists in the region $R_{flow-vort}^{(3)}$ when approaching it by increasing $c\, De$ from $R_{flow-vort}^{(2)}$. In the regions $R_{flow-vort}^{(2)}$ and $R_{flow-vort}^{(3)}$, the unstable direction of the saddle point leads to the limit cycle on one side, and the stable node on the other. Based on numerical evidence from figures 10 and 11 we assume the saddle point, $\boldsymbol {p}_{saddle}=[{p}_{1,saddle}\ {p}_{2,saddle}\ {p}_{3,saddle}]$, to approximately lie at the same $p_2$ coordinate as the stable node in the flow gradient plane:

(5.19a,b)\begin{equation} {p}_{2,saddle}={p}_{2,{flow}}^{0-}+{p}_{2,saddle}',\quad {p}_{2,saddle}' \ll {p}_{2,{flow}}^{0-} \end{equation}

where ${p}_{2,{flow}}^{0-}$ is given by (5.5a,b). From $\dot {p}_{2,saddle}=0$, we obtain

(5.20)\begin{equation} {p}_{2,saddle}=\frac{\sqrt{-b^2}c\, De+(4\log(2\kappa)-6)b^2}{c\, De(1-p_{1,saddle})-\sqrt{-b^2}(4\log(2\kappa)-6)}.\end{equation}

From (4.46), $\dot {p}_{1,saddle}=\dot {p}_{3,saddle}=0$ is equivalent to

(5.21)\begin{equation} p_{1,saddle} \alpha|_{p_1=p_{1,saddle},p_3=p_{3,saddle}}= p_{2,saddle}/(c\, De). \end{equation}

A general solution to the above equation is intractable, but we can solve it in the limit of small $p_{1,saddle}$ or small $p_{3,saddle}$. We denote the value of $p_{1,saddle}$ in the limit of small $p_{1,saddle}$ as $\tilde {p}_{1,saddle}$ and the value of $p_{3,saddle}$ in the limit of small $p_{3,saddle}$ as $\tilde {p}_{3,saddle}$. We obtain

(5.22) \begin{equation} \tilde{p}_{1,saddle}=\frac{\begin{array}{l}-8 (c\, De)^2 f^2+16 c \,De f^3 \sqrt{-b^2}\\ -\sqrt{\begin{array}{l} 4c \,De f \left(-(c \,De)^2 \kappa^2+2 (c \,De) f \kappa^2 \sqrt{-b^2}+4 f^2\right) \\ 4\left(-3 {\rm \pi}(c \,De)^2 \kappa^2+2 c \,De f \left(8 f+3 {\rm \pi}\kappa^2 \sqrt{-b^2}\right)+12 {\rm \pi}f^2\right)\\ +\,128 c \,De f^3 \left((c \,De)-2 f \sqrt{-b^2}\right)^2 \end{array}}\end{array}} {(c \,De) \left(3 {\rm \pi} (c \,De)^2 \kappa^2-2 (c \,De) f \left(8 f+3 {\rm \pi} \kappa^2 \sqrt{-b^2}\right)-12 {\rm \pi} f^2\right)},\end{equation}

and

(5.23) \begin{align} \tilde{p}_{3,saddle}=-\frac{\begin{array}{l} 3 {\rm \pi} (c \,De) \left((c \,De)^2 \kappa^2-2 c \,De f \kappa^2 \sqrt{-b^2}-4 f^2\right)\\ +\sqrt{\begin{array}{l} 9 {\rm \pi}^2 (c \,De)^2 \left(-(c \,De)^2 \kappa^2+2 c \,De f \kappa^2 \sqrt{-b^2}+4 f^2\right)^2\\ +\dfrac{1}{2}( c \,De)^3 \kappa^2 \left(3 {\rm \pi}^2 c \,De-2 f \left(3 {\rm \pi}^2 \sqrt{-b^2}+2\right)\right)\\ +\dfrac{1}{2} c \,Def^2\left(\,f \left(16-64 c \,De \sqrt{-b^2}\right)+c \,De \left(8 \kappa^2 \sqrt{-b^2}-12 {\rm \pi}^2\right)\right) \\ \left(3 {\rm \pi}^2 (c \,De)^2 \kappa^2-2 c \,De f \left(16 f+3 {\rm \pi}^2 \kappa^2 \sqrt{-b^2}\right)-4 f^2 \left(16 f \sqrt{-b^2}+3 {\rm \pi} ^2\right)\right)\end{array}}\end{array}}{2 \left(\dfrac{3}{4} {\rm \pi}^2 c \,De \left((c \,De)^2 \kappa^2-2 c \,De f \kappa^2 \sqrt{-b^2}-4 f^2\right)-8 c \,De f^2 \left((c \,De)+2 f \sqrt{-b^2}\right)\right)},\end{align}

where $f=2\log (2\kappa )-3$ and $\sqrt {-b^2}$ is in (5.12a,b). The variation of $\tilde {p}_{3,saddle}$ and $\sqrt {1-\tilde {p}_{1,saddle}^2-{p}_{2,saddle}^2}$ with $c\, De$ for a few $\kappa$ is shown in figure 14(a). From either of these approximations, we find that the saddle point moves closer to the vorticity axis as $c\, De$ is increased or $\kappa$ is reduced. Both approximations yield a similar location for the saddle point. The best approximation for the location of the saddle point for any $c\, De$ and $\kappa$, is $\boldsymbol {p}_{saddle}=[{p}_{1,saddle}\ {p}_{2,saddle}\ {p}_{3,saddle}]$, where

(5.24)\begin{gather} \left.\begin{gathered} {p}_{1,saddle}=\tilde{p}_{1,saddle},\quad {p}_{3,saddle}=\sqrt{1-\tilde{p}_{1,saddle}^2-{p}_{2,saddle}^2}, \quad\text{if }|\tilde{p}_{1,saddle}|\le|\tilde{p}_{3,saddle}|\\ {p}_{1,saddle}=\sqrt{1-{p}_{2,saddle}^2-\tilde{p}_{3,saddle}^2},\quad {p}_{3,saddle}=\tilde{p}_{3,saddle},\quad \text{if } |\tilde{p}_{1,saddle}|>|\tilde{p}_{3,saddle}| \end{gathered},\right\}\end{gather}

and $p_2$ is from (5.20). We show this saddle point's more accurate location of $p_3$ in figure 14(b). The location of this analytically predicted saddle point, shown with a green marker in figures 10, 11 and 15, agrees well with that inferred from the of numerically integrated trajectories.

Figure 14. The $p_3$ coordinate of the saddle point in the flow–vorticity plane in $R_{flow-vort}^{(2)}$ and $R_{flow-vort}^{(3)}$. (a) Solid lines (the graph of $\sqrt {1-\tilde {p}_{1,saddle}^2-p_{2,saddle}^2}$ from (5.22) and (5.20)) provide the more accurate $p_3$ location of the saddle point when it is closer to the vorticity axis ($p_{3,saddle}\approx 1$) and dashed lines (the graph of (5.23)) represent the more accurate $p_3$ location of the saddle point when it is closer to the flow axis ($p_{3,saddle}\approx 0$). (b) Most accurate location of the saddle point constructed from (a) using the crossover point between the solid and dashed lines for each $\kappa$, i.e. using (5.24).

Figure 15. Trajectories at large $c\, De$ for $\kappa =20$ in $R_{flow-vort}^{(3)}$, $R_{flow}^{(1)}$ and $R_{flow}^{(2)}$ zoomed near the flow–vorticity plane. At very large $c\, De$ (corresponding to regions $R_{flow}^{(1)}$ and $R_{flow}^{(2)}$), the stable limit cycle does not exist around the vorticity axis and a particle starting anywhere apart from vorticity axis ends up being nearly flow aligned.

In $R_{flow-vort}^{(2)}$ and $R_{flow-vort}^{(3)}$, the unstable manifold of this saddle point leads to two stable attractors, i.e. either a limit cycle near the vorticity axis or a stable node near the flow direction. Therefore, the location of the saddle point partially dictates the relative sizes of the basins of attraction of these stable attractors. As $c\, De$ increases, due to the changing location of the saddle point, $\boldsymbol {p}_{saddle}$, the particle is more likely to attain a final orientation towards the flow axis as compared with approaching the limit cycle around the vorticity axis (see figure 14 and compare the proportion of the solid lines in 10 and 11). The saddle point also persists in the large $c\, De$ regions $R_{flow}^{(1)}$ and $R_{flow}^{(2)}$ as shown in the different plots in figure 15.