# Rotation of a fiber in simple shear flow of a dilute polymer solution

Arjun Sharma1<EMAIL_ADDRESS>Donald L. Koch2  1Sibley School of
Mechanical and Aerospace Engineering, Cornell University, Ithaca, NY, 14853,
USA 2Robert Frederick Smith School of Chemical and Biomolecular Engineering,
Cornell University, Ithaca, NY, 14853, USA

###### Abstract

The motion of a freely rotating prolate spheroid in a simple shear flow of a
dilute polymeric solution is examined in the limit of large particle aspect
ratio, $\kappa$. A regular perturbation expansion in the polymer
concentration, $c$, a generalized reciprocal theorem, and slender body theory
to represent the velocity field of a Newtonian fluid around the spheroid are
used to obtain the $\mathcal{O}(c)$ correction to the particle’s orientational
dynamics. The resulting dynamical system predicts a range of orientational
behaviors qualitatively dependent upon $c\cdot De$ ($De$ is the imposed shear
rate times the polymer relaxation time) and $\kappa$ and quantitatively on
$c$. At a small but finite $c\cdot De$, the particle spirals towards a limit
cycle near the vorticity axis for all initial conditions. Upon increasing
$\kappa$, the limit cycle becomes smaller. Thus, ultimately the particle
undergoes a periodic motion around and at a small angle from the vorticity
axis. At moderate $c\cdot De$, a particle starting near the flow-gradient
plane departs it monotonically instead of spirally, as this plane (a limit
cycle at smaller $c\cdot De$) obtains a saddle and an unstable node. The
former is close to the flow direction. Upon further increasing $c\cdot De$,
the saddle-node changes to a stable node. Therefore, depending upon the
initial condition, a particle may either approach a periodic orbit near the
vorticity axis or obtain a stable orientation near the flow direction. Upon
further increasing $c\cdot De$, the limit cycle near the vorticity axis
vanishes, and the particle aligns with the flow direction for all starting
orientations.

###### keywords:

##  1 Introduction

A particle-filled viscoelastic polymeric fluid undergoes simple shear flow in
many industrial applications such as fiber spinning (Breitenbach, 2002; Huang
et al., 2003; Nakajima et al., 1994; Chae & Kumar, 2008) and roll-to-roll
manufacturing of high aspect ratio, low resistance films for flexible and
transparent electronics (Mutiso et al., 2013; Yin et al., 2010). Simple shear
flow occurs in the spinneret during fiber spinning and in the patterning
channel during roll-to-roll manufacturing. The suspending viscoelastic fluid
may include large aspect ratio particles to impart strength to the final
product in fiber spinning or provide a desired anisotropy to the low
resistance film. In the simple shear flow of an inertia-less Newtonian fluid,
a fiber/ slender particle undergoes an initial condition-dependent periodic
motion in orientational trajectories termed Jeffery (1922) orbits as shown in
figure 1 for particles with aspect ratio, $\kappa=10$ and 50. However, the
interaction between the polymers in a viscoelastic fluid and the fiber breaks
this degenerate periodic behavior. Previous experiments (Gauthier et al.,
1971; Bartram et al., 1975; Stover & Cohen, 1990; Johnson et al., 1990; Iso et
al., 1996b; Gunes et al., 2008) indicate that depending on the $\kappa$,
imposed shear rate and the properties of the viscoelastic fluid such as
polymer concentration and relaxation time, a particle may exhibit various
orientation dynamics. A slender particle, i.e., one with a large $\kappa$, may
either spiral or monotonically drift towards the vorticity axis, align near
the flow direction or settle somewhere within the flow-vorticity plane.
Therefore, careful design and choice of flow parameters during the simple
shear regime are essential in obtaining a final product with desired particle
orientation and material strength. Theoretical studies are useful due to the
many parameters required for characterizing a viscoelastic fluid. Polymers
lead to new features in a viscoelastic fluid flow such as shear thinning or a
finite first and second normal stress difference as compared to a Newtonian
fluid flow. Leal (1975) predicts that a slender particle in a slow flow will
spiral towards the vorticity axis due to the second normal stress difference
in the fluid. Whereas, operating in a double limit of small polymer
concentration and large Deborah number, $De$, (the product of the imposed
shear rate and the polymer relaxation time) Harlen & Koch (1993) also predict
the spiraling of the particle towards the vorticity axis, but due to first
normal stress difference in the fluid. Neither of these theories captures any
other orientation behavior observed experimentally. In this paper, using a
regular perturbation expansion in polymer concentration, $c$, we develop a
slender body theory that spans a range of $De$. It encapsulates the
$\mathcal{O}(c)$ effect of particle-polymer interaction and qualitatively
describes the rich orientation dynamics seen in previous experiments.

(a)  (b)

Figure 1: Jeffery (1922) orbits or orientation trajectories in a simple shear
flow of Newtonian fluid of a prolate spheroidal particle with aspect ratio,
$\kappa=$ (a) 10, and (b) 50.

Simple shear flow is ubiquitous in industrial applications as the near wall
flow can always be considered locally simple shear. Many scenarios include
laminar flow between two parallel walls; in these cases, the local flow is
simple shear everywhere. In many industrial scenarios, particle concentration
in the suspension is dilute, and due to negligible particle-particle
interaction, each particle can be considered to be suspended in an unbounded
fluid. In their experiments with highly elastic fluids, Iso et al. (1996a)
observed an isolated fiber to obtain a stable orientation close to the
direction of the imposed flow. They found that fibers in a moderate particle
concentration suspension also obtained a highly peaked orientation
distribution near the flow direction for the same parameters. Therefore, the
particle-particle interaction may sometimes be ignored even with higher
particle concentrations, making the studies of an isolated particle/ fiber
freely rotating in a simple shear flow of a viscoelastic fluid invaluable.

Based on their qualitative nature, the Jeffery (1922) orbits, illustrated for
particles with aspect ratio, $\kappa=10$ and 50 in figure 1, may be classified
into log-rolling, wobbling, flipping, and tumbling. Log-rolling occurs when
the particle is aligned with the axis of rotation of the imposed flow or the
vorticity axis. Here the particle rotates about its major axis. Except in the
log-rolling motion, the particle’s angular velocity changes throughout its
orbit. When initially placed in the flow-gradient plane, the particle remains
in the plane. It rotates about its minor axis while tumbling from one side of
the flow axis to the other. In the tumbling and flipping orbits, a large
aspect ratio particle or fiber spends only an $\mathcal{O}(1/\kappa)$
proportion (non-dimensionalized with shear rate) of the Jefferey time-period
of $2\pi\kappa$ away from the flow-vorticity plane indicated by dashed red
lines in the plots of figure 1. In the flipping orbits (which are three-
dimensional extensions of the tumbling orbit when the particle is not confined
to the flow-gradient plane), the particle comes within $\mathcal{O}(1/\kappa)$
of the flow direction. In these orbits, the particle’s orientation rapidly
flips from being aligned with the positive to the negative flow axis. During
its rapid flipping phase, a particle in a flipping orbit spans a large portion
of the orientation space in the gradient direction. In wobbling orbits, which
are smaller in circumference than flipping orbits, the particle does not come
very close to the flow direction. In these orbits, the particle gradually
wobbles in its orbit around the vorticity axis.

Gauthier et al. (1971) conducted experiments with $\kappa=16.1$ nylon rods in
a viscoelastic fluid made with 0.03 wt. % polyacrylamide solution in water.
They conducted experiments at a shear rate of $0.53s^{-1}$ and found that a
particle starting close to the flow-gradient plane spirals towards the
vorticity axis as it is exposed to the Couette flow. Using a similar
viscoelastic fluid and a polyethylene rod of $\kappa=9.0$ Bartram et al.
(1975) also found a similar behavior with shear rates up to $5s^{-1}$. Upon
further increase in shear rate, they found that a $\kappa=9.0$ rod released
near the gradient direction initially moves within the flow-gradient plane
towards an orientation near the flow direction. From here, introducing a
disturbance made the particle move monotonically along the flow-vorticity
plane away from the flow direction. When the particle was sufficiently close
to the vorticity direction, it started to spiral toward it. Bartram et al.
(1975) observed similar behavior with a $\kappa=5.6$ rod. However, unlike the
$\kappa=9.0$ rod, no disturbance was required when the particle came close to
the flow direction after being placed near the gradient direction. The time
period of particle rotation about the vorticity axis, that is already very
large at large $\kappa$ in Newtonian fluid ($2\pi\kappa$), is further
increased in experiments with viscoelastic fluids (Gauthier et al., 1971;
Bartram et al., 1975). In the experiments where the orientation of the
particle centerline was found to be spiraling towards the vorticity axis,
complete alignment with the vorticity axis was not shown.

Iso et al. (1996b) observed rotations of different high $\kappa$ isolated
fibers in a Boger fluid consisting of 1000 ppm polyisobutylene (PIB) in
polybutene (PB) in a simple shear flow with different shear rates. The polymer
relaxation time and concentration, $c$, defined as the polymer-to-solvent
viscosity ratio, were 3$s$ and 0.39, respectively. With a $\kappa=19$ fiber,
they found the particle to be orientated very close to the vorticity axis when
$De=1.5$. Increasing $\kappa$ to 34.4, or $De$ to 3.0 or both, they found
that, after initial spiraling away from the flow-gradient plane, the particle
obtained a steady orientation between 5∘ and 60 ∘ from the vorticity axis in
the flow-vorticity plane. With $\kappa=34.4$ and $De=3.0$, they report two
additional observations with no initial spiraling in contrast to other
experiments at identical parameters. The authors attributed different initial
orientations and fluid rheology due to slight changes in room temperature as
the causes of the lack of initial spiraling.

Across their two studies Iso et al. (1996a, b) also conducted experiments in a
viscoelastic liquid obtained by adding a certain amount of high molecular
weight polymer polyacrylamide (PAA) to a Newtonian solvent. The shear rate in
these experiments was $0.5s^{-1}$, and the fluid was slightly shear-thinning.
They observed various behaviors as the amount of PAA was increased from 100ppm
to 2000ppm (although the exact value of $c$ is not available, it increases
with PAA amount). For 100ppm ($\kappa=14$) and 500 ppm ($\kappa=24$) solutions
of PAA at a shear rate of 0.5$s^{-1}$ (Iso et al., 1996a, b), fibers either
end up in a trajectory where they oscillate in a small periodic orbit close to
the vorticity axis or obtain a stable orientation in the flow-vorticity plane
at a particular angle from the vorticity axis similar to the Boger fluid
(Boger, 1977) experiments at a higher shear rate of $1.0s^{-1}$ (Iso et al.,
1996b). With 1000 ppm PAA, the $\kappa=24$ fiber obtains a stable orientation
at the flow direction or 20∘ from the flow direction in the flow-vorticity
plane. Irrespective of the initial condition, a $\kappa=24$ fiber in 2000 ppm
PAA solution stably aligns with the flow direction. Therefore, fibers become
more flow aligned with increasing elasticity or polymer concentration.

The latest available experimental results measuring the effect of
viscoelasticity on the rotation of an anisotropic particle in simple shear
flow are by Gunes et al. (2008). They considered hematite spheroidal particles
with a much smaller aspect ratio, $\kappa$, between 2 and 7.5, than the
previous experimental studies. In a 20% hydroxypropylcellulose solution in
water, they found $\kappa=3.8$ particles to be oriented close to vorticity and
flow directions at low and large shear rates or $De$, respectively. At
moderate shear rates, particles exhibited a bimodal orientation distribution.
In a 2% poly-(ethyleneoxide), flow alignment was obtained at a lower shear
rate for a higher $\kappa$ or $c$. Most of the fluids they used were shear-
thinning. They reported one experiment of a non-shear-thinning Boger fluid
(Boger, 1977), consisting of a 0.1% polyisobutylene in polybutadiene solution,
with $\kappa=3.8$. Here the particles were close to the vorticity axis at all
shear rates reported.

Due to the variety of non-Newtonian parameters needed to fully characterize a
viscoelastic fluid and several physical phenomena that polymers may undergo
simultaneously, it is difficult to quantitatively compare the previous
experiments and identify the source of various behaviors of particle
orientation. For example, temperature changes during an experiment may change
the rheology of the fluid, subsequently changing the polymer’s relaxation
time, or multiple relaxation times may be needed to represent the fluid fully,
or adding more polymers to a solution may not only increase the polymer
concentration, but it may also change the relaxation time of polymers as they
entangle with one another. Thus numerical and theoretical modeling of the
relevant system is an important tool in isolating and understanding different
physical mechanisms that can complement or inspire future experiments.

Recently, d’Avino et al. (2014) reported numerical simulations of a
$\kappa=4.0$ prolate spheroidal particle in a simple shear flow of a Giesekus
fluid that models polymer melts (Bird et al., 1987), with large polymer
concentration. They used $c=10$ and found spiraling towards the vorticity axis
for $De\lesssim 1$ and alignment close to the flow direction for $De\gtrsim
3$. For intermediate $De$, depending on $De$, the particle obtained either one
or two stable orientations between the flow and vorticity axis. They mentioned
similar observations in unreported simulations with $\kappa=8$ spheroids with
the transition from vorticity- to flow-aligned particle orientation occurring
at a smaller $De$. Therefore, the trend of the particle’s final orientation
with the shear rate or $De$ and $\kappa$ between the experiments of Gunes et
al. (2008) and the numerical study of d’Avino et al. (2014), both conducted at
small $\kappa$, is similar. However, it is unclear from the aforementioned
numerical findings if the novel orientation dynamics are due to first or
second normal stress difference, shear thinning, or synergistic effects of
various non-Newtonian behaviors exhibited by the Giesekus fluid. Also, at
intermediate shear rates or $De$, while the orientation behavior is bimodal,
i.e., either along vorticity or flow directions, in the experiments of Gunes
et al. (2008), the final orientations in numerical results of d’Avino et al.
(2014) lie between the flow and vorticity axes in the flow-vorticity plane.
The latter is instead similar to some of the experimental observations of Iso
et al. (1996a) at larger $\kappa$.

A richer orientation behavior is illustrated in the previous experiments of
Gauthier et al. (1971); Bartram et al. (1975); Iso et al. (1996a, b) at larger
$\kappa$ as compared to the more recent studies of Gunes et al. (2008);
d’Avino et al. (2014) at smaller $\kappa$. Numerical studies with large
$\kappa$ in a viscoelastic fluid are expensive due to the large velocity and
polymer conformation gradients near the particle surface. Resolving these
large gradients and accurately modeling the shape of a slender particle
requires very fine spatial resolution and hence smaller time steps to ensure
numerical stability. Therefore, theoretical studies at large $\kappa$ are
invaluable, and Leal (1975); Harlen & Koch (1993); Abtahi & Elfring (2019) are
such pre-existing examples.

Using the slender body theory of Batchelor (1970), Leal (1975) predicts that a
fiber rotating in a slow, simple shear flow of a second order fluid, will
spiral towards the vorticity axis due to the second normal stress difference,
$\psi_{2}$, of the fluid. $\psi_{2}$ is usually smaller than the first normal
stress difference, $\psi_{1}$, and it is zero for Boger fluids (Boger, 1977).
Hence, Leal’s theory predicts that a slender particle rotating in a Boger
fluid undergoing a simple shear flow with a small shear rate rotates in the
same manner as in a Newtonian fluid. However, the low shear rate experiments
of Iso et al. (1996b) with a Boger fluid show spiraling towards vorticity. For
a large Deborah number, $De\gg 1$, small polymer concentration, $c\ll 1$ and
also using the slender body theory of Batchelor (1970), Harlen & Koch (1993)
predict the fiber in an Oldroyd-B fluid to spiral towards the vorticity axis,
but, due to $\psi_{1}$ (an Oldroyd-B fluid has $\psi_{2}=0$). Shear-thinning,
another property exhibited by polymeric fluids, does not play a role in either
of these theories. Abtahi & Elfring (2019) conducted a theoretical study on an
asymptotically weakly shear thinning liquid and found that a prolate spheroid
rotates in a closed periodic orbit but with a longer time period compared with
the Jeffery orbit in a Newtonian fluid. In other words, shear thinning makes a
prolate spheroid’s trajectory equivalent to that of a larger aspect ratio
particle in a Newtonian fluid but does not qualitatively alter the topology of
the trajectories.

Spiraling towards vorticity, as indicated by the two previous theories (Leal,
1975; Harlen & Koch, 1993) predicting a qualitative change in particle
orientation trajectory, is only one of the many qualitative influences of
viscoelasticity observed in the previous simple shear experiments of Gauthier
et al. (1971); Bartram et al. (1975); Iso et al. (1996a, b) at large particle
aspect ratio, $\kappa$. In this paper, assuming a small polymer concentration,
we theoretically revisit the effect of viscoelasticity on a slender fiber
rotating in a simple shear flow to explain the richer qualitative behavior of
a large $\kappa$ particle’s orientation in a polymeric fluid observed in the
experiments. We use the Oldroyd-B model to isolate the effect of elasticity
from shear thinning. Also, any predicted viscoelastic effects will exclude the
impact of the second normal stress difference and re-examine the fiber
rotation in Boger fluid undergoing simple shear flow.

In the absence of fluid inertia, fluid stress at any point in the viscoelastic
fluid surrounding a suspended particle is decomposed into three components: a)
Particle motion induced solvent stress, i.e., the stress around the particle
rotating in a Newtonian fluid, b) Elastic stress or the polymer stress, and,
c) Polymer-induced solvent stress, i.e., the stress created by the
perturbations in fluid velocity and pressure due to the forcing by polymer
stress. Therefore, the torque acting on the particle is the sum of particle
motion-induced solvent torque (MIST), elastic torque, and polymer-induced
solvent torque (PIST). We consider a freely rotating or torque-free particle
where the sum of the three torque components is zero.

The rest of the paper is organized as follows. In section 2 we discuss
different torque generating mechanisms along with the mathematical formulation
relevant to any freely rotating (torque-free) particle in a viscoelastic
fluid. For any polymer concentration, $c$, using a generalized reciprocal
theorem, we derive the formulation where PIST on any such particle can be
expressed in terms of the polymer stress instead of the polymer-induced
solvent stress. After section 2 we are concerned with viscoelastic fluid with
a small polymer concentration, $c$. Therefore, in section 3 we briefly
describe the regular perturbation expansion of the relevant flow variables in
$c$ and the procedure to obtain the $\mathcal{O}(c)$ correction to the
rotation of a torque-free particle in a low $c$ viscoelastic fluid. The
formulation in section 2 that expresses PIST in terms of the polymer stress
allows us to circumvent the numerical discretization of the partial
differential equations to calculate the $\mathcal{O}(c)$ polymer-induced
solvent stress. Therefore, the $\mathcal{O}(c)$ PIST (and also the elastic
torque) can be evaluated using the leading order or Newtonian velocity field
around the particle. In section 4, we use the matched asymptotic/ slender body
theory (SBT) solution for the Newtonian flow field around a slender prolate
spheroid to calculate the torques and obtain the $\mathcal{O}(c)$ correction
to the Jeffery (1922) rotation rates for large aspect ratio prolate spheroids
due to viscoelasticity. In SBT, in the inner region close to the particle,
owing to a large $\kappa$, the velocity field is obtained by assuming the flow
to vary slowly along the length of the particle compared to its variation
perpendicular to the particle surface. This solution is taken from Cox (1970).
Further from the particle surface, in the outer region, the particle center
line is replaced by a line of Stokeslets and doublets. The velocity
disturbance created by these are taken from Batchelor (1970), and Cox (1971)
respectively. In the SBT (Cox, 1970, 1971; Batchelor, 1970), the inner and the
outer solution approach one another in the matching reaching, i.e., in the
dual limit of the radial distance from the particle centerline written in
inner and outer variables approaching infinity and zero, respectively. In
section 5 we study the dynamical system defined by these equations for
different $c$ and $De$ and illustrate various orientation dynamics of the
spheroid predicted by this system. We also provide a qualitative comparison of
the theoretical orientation trajectories with the previous experimental
observations. Finally, we conclude our findings in section 6.

##  2 Mathematical formulation and different torque generating mechanisms in
viscoelastic fluids

In the absence of inertia, equations governing the conservation of mass and
momentum in the viscoelastic fluid surrounding a particle are,

$\nabla\cdot\mathbf{u}=0,\hskip 14.45377pt\nabla\cdot\boldsymbol{\sigma}=0,$
(1)

where $\mathbf{u}$ and $\boldsymbol{\sigma}$ are the fluid velocity vector and
stress tensor field. We consider a particle with its center of mass at the
origin of a simple shear flow such that it freely rotates with an angular
velocity $\boldsymbol{\omega}_{p}$, but does not translate. Therefore, the no-
slip boundary condition on particle surface and the far-field (as
$|\mathbf{r}|\rightarrow\infty$) imposed flow boundary condition are,

$\displaystyle{\mathbf{u}=\boldsymbol{\omega}_{p}\times\mathbf{r},\hskip
14.45377pt\text{ on particle surface },\hskip 14.45377pt\text{and,
}\mathbf{u}=\mathbf{r}\cdot(\nabla\mathbf{u})_{\infty},\hskip 14.45377pt\text{
as }|\mathbf{r}|\rightarrow\infty,}$ (2)

where $(\nabla\mathbf{u})_{\infty}$ is the velocity gradient tensor of the
imposed flow. The equations are non-dimensionalized with the particle’s
maximum length and the inverse of the imposed shear rate as length and time
scales. The fluid stress,

$\boldsymbol{\sigma}=\boldsymbol{\tau}+\mathbf{\Pi},$ (3)

is the sum of the solvent stress,
$\boldsymbol{\tau}=-p\mathbf{I}+\nabla\mathbf{u}+(\nabla\mathbf{u})^{T}$ and
the polymeric stress, $\boldsymbol{\Pi}$. The solvent viscosity is the
viscosity scale in our non-dimensionalization. In the solvent stress, $p$ is
the fluid pressure. We use the Oldroyd-B model (Bird et al., 1987) where the
polymers are modeled as dumbells consisting of two beads attached to a
linearly elastic spring. The polymeric stress is,

$\boldsymbol{\Pi}=\frac{c}{De}(\boldsymbol{\Lambda}-\mathbf{I}),$ (4)

where $c$ is the polymer concentration, $De$ is a non-dimensional product of
the polymer relaxation time and imposed shear rate, $\boldsymbol{\Lambda}$ is
the polymer conformation (defined as the average over possible polymer
conformations of the outer product of the polymer end to end vector with
itself) and $\mathbf{I}$ is the identity tensor. $\boldsymbol{\Lambda}$ is
affected by convection due to the fluid velocity, stretching and rotation by
the velocity gradients, and relaxation of the polymer to its equilibrium
orientation, $\boldsymbol{\Lambda}=\mathbf{I}$. It is governed by,

$\frac{\partial\boldsymbol{\Lambda}}{\partial
t}+\mathbf{u}\cdot\nabla\boldsymbol{\Lambda}=(\nabla\mathbf{u})^{T}\cdot\boldsymbol{\Lambda}+\boldsymbol{\Lambda}\cdot\nabla\mathbf{u}+\frac{1}{De}(\mathbf{I}-\boldsymbol{\Lambda}).$
(5)

Due to the absence of any non-linear term explicitly involving the velocity
and pressure variables in the mass and momentum conservation equations, we can
split equation (1) and its boundary conditions into two parts. The first part
is the same as the flow around a particle rotating with an angular velocity
$\boldsymbol{\omega}_{p}$ in an imposed simple shear flow of a Newtonian
fluid. It is governed by,

${\nabla\cdot\mathbf{u}^{\text{M}}=0,\hskip
14.45377pt\nabla\cdot\boldsymbol{\tau}^{\text{M}}=0,}$ (6)

where the particle motion induced solvent stress,
$\boldsymbol{\tau}^{\text{M}}=-p^{\text{M}}\mathbf{I}+\nabla\mathbf{u}^{\text{M}}+(\nabla\mathbf{u}^{\text{M}})^{T}$,
is the sum of the pressure and rate of strain tensor in the solvent due to the
particle motion in an inertia-less Newtonian fluid. These equations are
subject to boundary conditions,

$\displaystyle\mathbf{u}^{\text{M}}=\boldsymbol{\omega}_{p}\times\mathbf{r},\hskip
14.45377pt\text{ on particle surface},\hskip 14.45377pt\text{and,
}\mathbf{u}^{\text{M}}=\mathbf{r}\cdot(\nabla\mathbf{u})_{\infty},\hskip
14.45377pt\text{ as }|\mathbf{r}|\rightarrow\infty.$ (7)

The second part is the balance of the divergence of the sum of the polymeric
stress and $\boldsymbol{\tau}^{\text{P}}$,

$\nabla\cdot\mathbf{u}^{\text{P}}=0,\hskip
14.45377pt\nabla\cdot\boldsymbol{\tau}^{\text{P}}+\nabla\cdot\boldsymbol{\Pi}=0,$
(8)

where the solvent stress generated due to the forcing by the polymeric stress
is
$\boldsymbol{\tau}^{\text{P}}=-p^{\text{P}}\mathbf{I}+\nabla\mathbf{u}^{\text{P}}+(\nabla\mathbf{u}^{\text{P}})^{T}$.
$p^{\text{P}}$ and
$(\nabla\mathbf{u}^{\text{P}}+(\nabla\mathbf{u}^{\text{P}})^{T})/2$ are the
modification of pressure and rate of strain tensor by the polymers. Boundary
conditions for equation (8) are,

$\displaystyle\mathbf{u}^{\text{P}}=0,\hskip 14.45377pt\text{ on particle
surface},\hskip 14.45377pt\text{and, }\mathbf{u}^{\text{P}}=0,\hskip
14.45377pt\text{ as }|\mathbf{r}|\rightarrow\infty.$ (9)

The two parts are coupled via equation (5), i.e., the polymer constitutive
equation, where

$\mathbf{u}=\mathbf{u}^{\text{M}}+\mathbf{u}^{\text{P}},$ (10)

and the torque-free condition (equation (17)). In this framework, the total
fluid stress,

$\boldsymbol{\sigma}=\boldsymbol{\tau}^{\text{M}}+\boldsymbol{\tau}^{\text{P}}+\boldsymbol{\Pi},$
(11)

and the total torque acting on the particle surface, $\boldsymbol{\sigma}$,

$\mathbf{G}(\boldsymbol{\sigma})=\int_{\mathbf{r}_{\text{p}}}dS\hskip
7.22743pt\mathbf{r}\times(\boldsymbol{\sigma}\cdot\mathbf{n}),$ (12)
$\mathbf{G}(\boldsymbol{\sigma})=\mathbf{G}(\boldsymbol{\tau}^{\text{M}})+\mathbf{G}(\boldsymbol{\tau}^{\text{P}})+\mathbf{G}(\boldsymbol{\Pi})=\mathbf{G}^{\text{MIST}}+\mathbf{G}^{\text{PIST}}+\mathbf{G}^{\text{Elastic}},$
(13)

are decomposed into three underlying mechanisms, where,

$\displaystyle\mathbf{G}^{\text{MIST}}$
$\displaystyle=\mathbf{G}({\boldsymbol{\tau}^{\text{M}}})=\int_{\mathbf{r}_{\text{p}}}dS\hskip
7.22743pt\mathbf{r}\times\boldsymbol{\tau}^{\text{M}}\cdot\mathbf{n},$ (14)
$\displaystyle\mathbf{G}^{\text{PIST}}$
$\displaystyle=\mathbf{G}({\boldsymbol{\tau}^{\text{P}}})=\int_{\mathbf{r}_{\text{p}}}dS\hskip
7.22743pt\mathbf{r}\times\boldsymbol{\tau}^{\text{P}}\cdot\mathbf{n},\text{
and},$ (15) $\displaystyle\mathbf{G}^{\text{Elastic}}$
$\displaystyle=\mathbf{G}(\boldsymbol{\Pi})=\int_{\mathbf{r}_{\text{p}}}dS\hskip
7.22743pt\mathbf{r}\times\boldsymbol{\Pi}\cdot\mathbf{n},$ (16)

are the particle motion induced solvent, the polymer induced solvent, and the
elastic or polymeric torques, respectively. The angular velocity,
$\boldsymbol{\omega}_{p}$, of a freely rotating particle, is determined by the
torque-free condition,

$\boldsymbol{\omega}_{p}=\\{\boldsymbol{\omega}_{p}:\mathbf{G}=\mathbf{G}^{\text{MIST}}+\mathbf{G}^{\text{PIST}}+\mathbf{G}^{\text{Elastic}}=0\\}.$
(17)

We consider the motion of a freely rotating particle in a viscoelastic fluid.
Our main motivation is to study the behavior of a prolate spheroid in simple
shear flow. Due to symmetry, this particle has zero net hydrodynamic force if
it moves with the local velocity of the imposed flow. Its physically motivated
decomposed components (particle motion induced solvent, polymer induced
solvent, and elastic forces) are individually zero by symmetry. However, the
force balance can be considered similar to the torque balance discussed above
for determining the motion of a freely translating particle (of any shape) in
a general linear flow or a particle sedimenting under gravity (where the net
hydrodynamic force must balance the force due to gravity).

###  2.1 Using a generalized reciprocal theorem to obtain
${\mathbf{G}^{\text{PIST}}}$ in terms of $\boldsymbol{\Pi}$

In its original form,
$\mathbf{G}^{\text{PIST}}=\int_{\mathbf{r}_{\text{p}}}dS\hskip
7.22743pt\mathbf{r}\times\boldsymbol{\tau}^{\text{P}}\cdot\mathbf{n}$ is a
function of $\boldsymbol{\tau}^{\text{P}}$ which in turn is driven by
$\boldsymbol{\Pi}$ through equation (8). In this section, using the balance of
the divergence of the polymeric stress ($\boldsymbol{\Pi}$) and
$\boldsymbol{\tau}^{\text{P}}$ from equation (8) and using a generalized
reciprocal theorem we derive an expression to obtain
$\mathbf{G}^{\text{PIST}}$ directly from $\boldsymbol{\Pi}$ without the need
to compute $\boldsymbol{\tau}^{\text{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-2 pseudo-tensor $\mathbf{b}$, ‘pressure’ field that is a pseudovector,
$\mathbf{q}$, and ‘fluid stress’ field that is a rank-3 pseudo-tensor
$\mathbf{B}$,

$\displaystyle\frac{\partial B_{ijk}}{\partial r_{i}}=0,\hskip
14.45377pt\frac{\partial b_{ij}}{\partial r_{i}}=0,\hskip
14.45377ptB_{ijk}=-\delta_{ij}q_{k}+\frac{\partial b_{jk}}{\partial
r_{i}}+\frac{\partial b_{ik}}{\partial r_{j}},$ (18)

with boundary condition

$b_{ij}=\epsilon_{ijk}r_{k},\hskip 14.45377pt\text{on particle surface},\hskip
14.45377pt\text{and},\hskip 14.45377ptb_{ij}=0,\hskip 14.45377pt\text{as
}|\mathbf{r}|\rightarrow\infty,$ (19)

where $\epsilon_{ijk}$ is the permutation tensor.
$\mathbf{b}\cdot\boldsymbol{\omega}_{\text{auxillary}}$ is the velocity field
around a particle rotating with an angular velocity
$\boldsymbol{\omega}_{\text{auxillary}}$ in a quiescent Newtonian fluid. Using
the definitions of $\boldsymbol{\tau}^{\text{P}}$ and $\mathbf{B}$ in terms of
the respective velocities, $\mathbf{u}^{\text{P}}$ and $\mathbf{b}$,
incompressibility of the velocities in equations (8) and (18),
$\nabla\cdot\mathbf{B}=0$ from equation (18) and the symmetry of $B_{lki}$ and
${\tau}^{\text{P}}_{lk}$ about the $l$ and $k$ indices we obtain,

${{\tau}^{\text{P}}_{lk}\frac{\partial b_{ki}}{\partial r_{l}}=\frac{\partial
B_{lki}{u}^{\text{P}}_{k}}{\partial r_{l}}.}$ (20)

Using the balance of the divergence of $\boldsymbol{\Pi}$ and
$\boldsymbol{\tau}^{\text{P}}$ from equation (8), the volume integral of
equation (20) in a region bounded by the particle surface,
$\mathbf{r}_{\text{p}}$, and a far-field spherical surface at
$|\mathbf{r}|\rightarrow\infty$ is

${\int_{\text{Fluid}}dV\hskip 7.22743pt\frac{\partial}{\partial
r_{l}}[{\tau}^{\text{P}}_{lk}b_{ki}-B_{lki}{u}^{\text{P}}_{k}]=-\int_{\text{Fluid}}dV\hskip
7.22743ptb_{ki}\frac{\partial{\Pi}_{lk}}{\partial r_{l}}.}$ (21)

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

$\displaystyle\begin{split}&\int_{\text{Fluid}}dV\hskip
7.22743pt\frac{\partial}{\partial
r_{l}}[{\tau}^{\text{P}}_{lk}b_{ki}-B_{lki}{u}^{\text{P}}_{k}]=\\\
&\int_{|\mathbf{r}|\rightarrow\infty}dS\hskip
7.22743ptn_{l}[{\tau}^{\text{P}}_{lk}b_{ki}-B_{lki}{u}^{\text{P}}_{k}]-\int_{\mathbf{r}_{\text{p}}}dS\hskip
7.22743ptn_{l}[{\tau}^{\text{P}}_{lk}b_{ki}-B_{lki}{u}^{\text{P}}_{k}],\end{split}$
(22)

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 $\mathbf{u}^{\text{P}}$ and $\mathbf{b}$
scale as $1/r$ and $1/{r}^{2}$, respectively, and the stresses
$\boldsymbol{\tau}^{\text{P}}$ and $\mathbf{B}$ scale as $1/{r}^{2}$ and
$1/{r}^{3}$. Therefore, the first surface integral in equation (22) vanishes.
Using the boundary conditions for $\mathbf{u}^{\text{P}}$ and $\mathbf{b}$
from equations (9) and (19) in the second surface integral in equation (22) we
obtain,

$\int_{\text{Fluid}}dV\hskip 7.22743pt\frac{\partial}{\partial
r_{l}}({\tau}^{\text{P}}_{lk}b_{ki}-B_{lki}{u}^{\text{P}}_{k})=\int_{\mathbf{r}_{\text{p}}}dS\hskip
7.22743ptn_{l}{\tau}^{\text{P}}_{lk}\epsilon_{kim}r_{m}=-{G}^{\text{PIST}}_{i}.$
(23)

Therefore, from equation (21) and (23), for a particle of any shape,

$\boldsymbol{G}^{\text{PIST}}=\int_{\text{Fluid}}dV\hskip
7.22743pt(\nabla\cdot{\boldsymbol{\Pi}})\cdot\mathbf{b}.$ (24)

This completes the derivation expressing $\boldsymbol{G}^{\text{PIST}}$ in
terms of the polymer stress $\boldsymbol{\Pi}$ and a quasi-steady 2-tensor
field $\mathbf{b}$ that is dependent on the particle shape and is a solution
of the auxiliary Stokes problem defined by equations (18) and (19).

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

$\boldsymbol{G}^{\text{PIST}}=\int_{\text{Fluid}}dV\hskip
7.22743pt(\nabla\cdot{(\boldsymbol{\Pi}-\boldsymbol{\Pi}^{U})})\cdot\mathbf{b}.$
(25)

Using the chain rule and divergence theorem,

${G}^{\text{PIST}}_{i}=\int_{|\mathbf{r}|\rightarrow\infty}dS\hskip
0.72229pt({\Pi}_{lk}-{\Pi}^{U}_{lk})b_{ki}n_{l}-\int_{\mathbf{r}_{\text{p}}}dS\hskip
0.72229pt({\Pi}_{lk}-{\Pi}^{U}_{lk})b_{ki}n_{l}-\int_{\text{Fluid}}dV\hskip
0.72229pt({\Pi}_{lk}-{\Pi}^{U}_{lk})\frac{\partial b_{ki}}{\partial r_{l}}.$
(26)

The disturbance of the polymer stress created by the particle,
${\Pi}_{lk}-{\Pi}^{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., 1987), this far-field scaling of ${\Pi}_{lk}-{\Pi}^{U}_{lk}$ is
ascertained by linearizing the polymer constitutive equation (for example (4)
and (5) for the Oldroyd-B model) about the far-field velocity and polymer
conformation to obtain a governing equation for ${\Pi}_{lk}-{\Pi}^{U}_{lk}$.
Therefore, the first surface integral in the above equation vanishes, and
using the surface boundary condition of equation (19) leads to,

${G}^{\text{PIST}}_{i}=-\int_{\mathbf{r}_{\text{p}}}dS\hskip
7.22743pt({\Pi}_{lk}-{\Pi}^{U}_{lk})\epsilon_{kif}r_{f}n_{l}-\int_{\text{Fluid}}dV\hskip
7.22743pt({\Pi}_{lk}-{\Pi}^{U}_{lk})\frac{\partial b_{ki}}{\partial r_{l}},$
(27)

and using ${G}^{\text{Elastic}}_{i}=\int_{\mathbf{r}_{\text{p}}}dS\hskip
7.22743pt{\Pi}_{lk}\epsilon_{kif}r_{f}n_{l}$,

${G}^{\text{PIST}}_{i}+{G}^{\text{Elastic}}_{i}=\int_{\mathbf{r}_{\text{p}}}dS\hskip
7.22743pt{\Pi}^{U}_{lk}\epsilon_{kif}r_{f}n_{l}-\int_{\text{Fluid}}dV\hskip
7.22743pt({\Pi}_{lk}-{\Pi}^{U}_{lk})\frac{\partial b_{ki}}{\partial r_{l}}.$
(28)

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

$\displaystyle{\boldsymbol{G}^{\text{PIST}}}+{\boldsymbol{G}^{\text{Elastic}}}=-\int_{\text{Fluid}}dV\hskip
7.22743pt(\boldsymbol{\Pi}-\boldsymbol{\Pi}^{U}):\nabla\mathbf{b},\text{
and,}$ (29)
$\displaystyle{\boldsymbol{G}^{\text{Elastic}}}=-\int_{\text{Fluid}}dV\hskip
7.22743pt\nabla\cdot((\boldsymbol{\Pi}-\boldsymbol{\Pi}^{U})\cdot\mathbf{b}).$
(30)

Equation (28) (or (29) for a fore-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{\Pi}$, and its
undisturbed value, $\boldsymbol{\Pi}^{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,
${\Pi}^{U}_{ij}=c(2De\delta_{i1}\delta_{j1}+(\delta_{i2}\delta_{j1}+\delta_{i1}\delta_{j2}))$.

##  3 Regular perturbation expansion for small polymer concentration

The leading order solution in a regular perturbation expansion for $c\ll 1$
corresponds to a freely rotating particle in simple shear flow of a Newtonian
fluid. At this order the stresses and hence the respective torques arising due
to the polymers, i.e. ${\boldsymbol{G}^{\text{PIST}}}$ and
${\boldsymbol{G}^{\text{Elastic}}}$ are zero and the particle rotates with an
angular velocity that satisfies ${\boldsymbol{G}^{\text{MIST}}}=0$ (equation
(17)). This is simply the Jeffery (1922) rotation. At the leading order the
polymer configuration is driven by the leading order velocity field (equation
(5)) and this leads to an $\mathcal{O}(c)$ polymer stress (equation (4)).
Therefore, the torques ${\boldsymbol{G}^{\text{PIST}}}$ and
${\boldsymbol{G}^{\text{Elastic}}}$ are $\mathcal{O}(c)$ (equations (25),
(14), (15), (16)). Hence, the particle rotation must be modified at
$\mathcal{O}(c)$ such that the sum all three torques
${\boldsymbol{G}^{\text{MIST}}}$, ${\boldsymbol{G}^{\text{PIST}}}$ and
${\boldsymbol{G}^{\text{Elastic}}}$ is zero at $\mathcal{O}(c)$. The regular
perturbation expansion of the relevant flow variables is,

$\displaystyle\mathbf{u}^{\text{M}}=$
$\displaystyle{\mathbf{u}^{\text{M}}}^{(0)}+c{\mathbf{u}^{\text{M}}}^{(1)}+\mathcal{O}(c^{2}),$
(31) $\displaystyle p^{\text{M}}=$
$\displaystyle{p^{\text{M}}}^{(0)}+c{p^{\text{M}}}^{(1)}+\mathcal{O}(c^{2}),$
(32) $\displaystyle\boldsymbol{\tau}^{\text{M}}=$
$\displaystyle{\boldsymbol{\tau}^{\text{M}}}^{(0)}+c{\boldsymbol{\tau}^{\text{M}}}^{(1)}+\mathcal{O}(c^{2}),$
(33) $\displaystyle\mathbf{u}^{\text{P}}=$ $\displaystyle
c{\mathbf{u}^{\text{P}}}^{(1)}+\mathcal{O}(c^{2}),$ (34) $\displaystyle
p^{\text{P}}=$ $\displaystyle c{p^{\text{P}}}^{(1)}+\mathcal{O}(c^{2}),$ (35)
$\displaystyle\boldsymbol{\tau}^{\text{P}}=$ $\displaystyle
c{\boldsymbol{\tau}^{\text{P}}}^{(1)}+\mathcal{O}(c^{2}),$ (36)
$\displaystyle\boldsymbol{\Lambda}=$
$\displaystyle\boldsymbol{\Lambda}^{(0)}+c\boldsymbol{\Lambda}^{(1)}+\mathcal{O}(c^{2}),$
(37) $\displaystyle\boldsymbol{\Pi}=$ $\displaystyle
c\boldsymbol{\Pi}^{(1)}+\mathcal{O}(c^{2}),$ (38)
$\displaystyle\boldsymbol{\omega}_{p}=$
$\displaystyle{\boldsymbol{\omega}^{(0)}_{p}}+c{\boldsymbol{\omega}^{(1)}_{p}}+\mathcal{O}(c^{2}),$
(39) $\displaystyle\mathbf{G}^{\text{MIST}}=$ $\displaystyle
c{\mathbf{G}^{\text{MIST}}}^{(1)}+\mathcal{O}(c^{2}),$ (40)
$\displaystyle\mathbf{G}^{\text{PIST}}=$ $\displaystyle
c{\mathbf{G}^{\text{PIST}}}^{(1)}+\mathcal{O}(c^{2}),$ (41)
$\displaystyle\mathbf{G}^{\text{Elastic}}=$ $\displaystyle
c{\mathbf{G}^{\text{Elastic}}}^{(1)}+\mathcal{O}(c^{2}).$ (42)

In an inertia-less Newtonian fluid undergoing simple shear a particle rotating
at an angular velocity $\boldsymbol{\omega}^{(0)}_{p}$ generates the velocity
field ${\mathbf{u}^{\text{M}}}^{(0)}$. As mentioned earlier, the leading order
angular velocity, ${\boldsymbol{\omega}^{(0)}_{p}}$, is the Jeffery (1922)
rotation. It allows the leading order torque free condition,
$\int_{\mathbf{r}_{\text{p}}}dS\hskip
7.22743pt\mathbf{r}\times{\boldsymbol{\tau}^{\text{M}}}^{(0)}\cdot\mathbf{n}=0$
to be satisfied. The leading order polymer constitutive equation is

$\frac{\partial\boldsymbol{\Lambda}^{(0)}}{\partial
t}+{\mathbf{u}^{\text{M}}}^{(0)}\cdot\nabla\boldsymbol{\Lambda}^{(0)}=(\nabla{\mathbf{u}^{\text{M}}}^{(0)})^{T}\cdot\boldsymbol{\Lambda}^{(0)}+\boldsymbol{\Lambda}^{(0)}\cdot\nabla{\mathbf{u}^{\text{M}}}^{(0)}+\frac{1}{De}(\mathbf{I}-\boldsymbol{\Lambda}^{(0)}).$
(43)

Solving this equation, one obtains the $\mathcal{O}(c)$ polymer stress,

$\displaystyle\boldsymbol{\Pi}^{(1)}=\frac{1}{De}(\boldsymbol{\Lambda}^{(0)}-\mathbf{I}),$
(44)

and hence the $\mathcal{O}(c)$ polymer induced solvent and elastic torques,

$\displaystyle{\mathbf{G}^{\text{PIST}}}^{(1)}=$
$\displaystyle\int_{\text{Fluid}}dV\hskip
14.45377pt\nabla\cdot(\boldsymbol{\Pi}^{(1)}-\boldsymbol{\Pi}^{U}/c)\cdot\mathbf{b},$
(45) $\displaystyle{\mathbf{G}^{\text{Elastic}}}^{(1)}=$
$\displaystyle\int_{\mathbf{r}_{\text{p}}}dS\hskip
7.22743pt\mathbf{r}\times\boldsymbol{\Pi}^{(1)}\cdot\mathbf{n}.$ (46)

The $\mathcal{O}(c)$ angular velocity, ${\boldsymbol{\omega}^{(1)}_{p}}$ is
the one that allows the $\mathcal{O}(c)$ torque-free condition to be
satisfied,

${\mathbf{G}^{\text{MIST}}}^{(1)}=\int_{\mathbf{r}_{\text{p}}}dS\hskip
7.22743pt\mathbf{r}\times{\boldsymbol{\tau}^{\text{M}}}^{(1)}\cdot\mathbf{n}=-{\mathbf{G}^{\text{PIST}}}^{(1)}-{\mathbf{G}^{\text{Elastic}}}^{(1)},$
(47)

where
${\boldsymbol{\tau}^{\text{M}}}^{(1)}=-{p^{\text{M}}}^{(1)}\mathbf{I}+\nabla{\mathbf{u}^{\text{M}}}^{(1)}+(\nabla{\mathbf{u}^{\text{M}}}^{(1)})^{T}$
is obtained from the $\mathcal{O}(c)$ equations,

$\nabla\cdot{\mathbf{u}^{\text{M}}}^{(1)}=0,\hskip
14.45377pt\nabla\cdot{\boldsymbol{\tau}^{\text{M}}}^{(1)}=0,$ (48)

subject to the velocity boundary conditions,

$\displaystyle{\mathbf{u}^{\text{M}}}^{(1)}={\boldsymbol{\omega}^{(1)}_{p}}\times\mathbf{r},\hskip
14.45377pt\text{ on particle surface},\hskip 14.45377pt\text{and,
}{\mathbf{u}^{\text{M}}}^{(1)}=0,\hskip 14.45377pt\text{ as
}|\mathbf{r}|\rightarrow\infty.$ (49)

Equations (48) and (49) represent the rotation of the particle in a quiescent
Newtonian fluid i.e. a Stokes flow. Due to the linearity of the Stokes flow,
${\mathbf{G}^{\text{MIST}}}^{(1)}$ is a linear function of
${\boldsymbol{\omega}^{(1)}_{p}}$. The latter can be viewed as the additional
angular velocity of the particle that generates a large enough viscous torque,
${\mathbf{G}^{\text{MIST}}}^{(1)}$, to balance the sum of polymer induced
solvent and elastic torque.

The formulation in this section is valid for any polymer constitutive model
and particle shape. By using the formulation of $\mathbf{G}^{\text{PIST}}$ in
equation (25) to express its $\mathcal{O}(c)$ value in equation (46), we have
avoided dealing with the $\mathcal{O}(c)$ equation for the balance of the
divergence of the polymer stress and the polymer induced solvent stress
(obtained by regularly expanding equation (8) in $c$). Otherwise, obtaining
$\mathbf{G}^{\text{PIST}}$ would have required a numerical solution via
discretization of the governing partial differential equations.

##  4 Rotation of a fiber due to simple shear flow in viscoelastic fluid with
small polymer concentration, $c$

In this section, we calculate torques from the various physical mechanisms
discussed in section 2 on a large aspect ratio prolate spheroid (considered a
slender fiber) freely rotating in a simple shear flow of a viscoelastic fluid
with a small polymer concentration, $c$, using the procedure indicated in
section 3. These torques are then balanced to obtain the $\mathcal{O}(c)$
correction to the particle’s rotation rate due to the presence of the
polymers. The Jeffery orbit period of a slender fiber with aspect ratio
$\kappa$ is $2\pi\kappa$ and the proportion of time spent outside
$|p_{2}|>\mathcal{O}(1/\kappa)$ is only $\mathcal{O}(1/\kappa)$, where
$p_{2}=0$ defines the flow-vorticity plane. Therefore, a slender fiber
suspended in a Newtonian fluid spends most of its Jeffery (1922) orbit close
to the flow-vorticity plane, where the particle rotation rate is very small.
Hence, most of the elasticity influence arises when the particle is in this
orientation, and polymer conformation when the fiber is close to the flow-
vorticity plane is quasi-steady. Thus,the polymer constitutive equation (43)
is simplified to,

${\frac{\partial\boldsymbol{\Lambda}^{(0)}}{\partial t}\approx
0\rightarrow{\mathbf{u}^{\text{M}}}^{(0)}\cdot\nabla\boldsymbol{\Lambda}^{(0)}\approx(\nabla{\mathbf{u}^{\text{M}}}^{(0)})^{T}\cdot\boldsymbol{\Lambda}^{(0)}+\boldsymbol{\Lambda}^{(0)}\cdot\nabla{\mathbf{u}^{\text{M}}}^{(0)}+\frac{1}{De}(\mathbf{I}-\boldsymbol{\Lambda}^{(0)}).}$
(50)

The Stokes flow solution of a Newtonian fluid around a large aspect ratio or a
slender prolate spheroid is analytically solvable via a matched asymptotic
expansion in $1/\kappa$ called slender body theory (SBT). In SBT, the fluid
velocity close to the particle or in the inner region is considered quasi-two-
dimensional. It is solved by ignoring the end effects and treating the slender
particle as an infinite cylinder. In terms of the radial distance from the
centerline (non-dimensionalized with the major radius of the particle) $\rho$,
the inner region is defined by $\rho\ll 1$. Away from the particle in the
outer region, defined by $\rho\gg 1/\kappa$ ($1/\kappa$ is the minor radius),
the particle is assumed to be a line of point forces (Batchelor, 1970; Cox,
1970) and force doublets (Cox, 1971). These singularity solutions are used to
represent the velocity and pressure fields. For a large enough $\kappa$ (which
is required for SBT to be valid), an intermediate or matching region exists
which is defined by $1/\kappa\ll\rho\ll 1$. In this region, the
$\rho\rightarrow\infty$ asymptote of the flow in the inner region and
$\rho\rightarrow 0$ asymptote of the flow in the outer region are matched. See
previous SBT calculations in (Cox, 1970, 1971; Batchelor, 1970) for more
details. In the rest of this section, we will use these SBT results to obtain
the effect of viscoelasticity on the rotation of a slender prolate spheroid.
In contrast to a fiber with blunt ends, such as a cylinder, a slender prolate
spheroid is more convenient for analysis due to the absence of localized
forces at the ends.

###  4.1 Flow of Newtonian fluid around a slender fiber and particle rotation
rates

Due to the relative velocity between the particle’s centerline and the imposed
flow, the flow disturbance at $\mathcal{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 centerline. This flow is
considered by the general slender body theory (SBT) of Cox (1970) up to
$\mathcal{O}(1/\log(\kappa)^{2})$. An equivalent theory by Batchelor (1970) 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 (1970)
instead of Cox (1970). For a prolate spheroidal particle fixed in the flow-
vorticity plane of a simple shear flow or rotating about its centerline in a
quiescent fluid, the force Stokeslet from these two theories is zero. The flow
driven by the local velocity gradients, which acts at
$\mathcal{O}(1/\kappa^{2})$, is dominant in these cases. Considering the
velocity gradients in the far-field relative to the particle, Cox (1971)
provides the solution for this flow. It is a combination of flows driven by
$\mathcal{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 fiber
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 $\mathbf{p}$,

$\mathbf{p}=\begin{bmatrix}p_{1}&p_{2}&p_{3}\end{bmatrix}^{T},$ (51)

close to the flow-vorticity plane is obtained by superposition of flows from
Batchelor (1970), and Cox (1971) with the assumption $p_{2}\ll 1$ and
$\kappa\gg 1$. The outer flow in the presence of the fiber is

$\mathbf{u}_{\text{out}}^{\text{M}{(0)}}(\mathbf{r};\boldsymbol{\omega}^{(0)}_{p},\mathbf{p})=\mathbf{r}\cdot(\nabla\mathbf{u})_{\infty}+{\mathbf{u}_{\text{out}}^{\text{M}{(0)}}}^{\prime}(\mathbf{r};\boldsymbol{\omega}^{(0)}_{p},\mathbf{p}).$
(52)

The disturbance velocity field,

$\displaystyle{\mathbf{u}_{\text{out}}^{\text{M}{(0)}}}^{\prime}(\mathbf{r};\boldsymbol{\omega}^{(0)}_{p},\mathbf{p})=\frac{1}{8\pi}\int_{-1}^{1}d\lambda\hskip
7.22743pt\Big{\\{}\frac{\boldsymbol{h}^{(0)}(\lambda;\boldsymbol{\omega}^{(0)}_{p},\mathbf{p})}{\log(2\kappa)-1.5}+\frac{\mathbf{f}^{(0)}(\lambda;\mathbf{p})}{\kappa^{2}}\Big{\\}}\cdot\Big{\\{}\frac{\mathbf{I}}{|\mathbf{r}-\lambda\mathbf{p}|}+\frac{(\mathbf{r}-\lambda\mathbf{p})(\mathbf{r}-\lambda\mathbf{p})}{|\mathbf{r}-\lambda\mathbf{p}|^{3}}\Big{\\}}+$
$\displaystyle\frac{1}{8\pi}\int_{-1}^{1}d\lambda\hskip
7.22743pt\frac{1}{\kappa^{2}}\boldsymbol{\mathcal{G}}^{(0)}(\lambda;\boldsymbol{\omega}^{(0)}_{p},\mathbf{p}):\nabla\Big{\\{}\frac{\mathbf{I}}{|\mathbf{r}-\lambda\mathbf{p}|}+\frac{(\mathbf{r}-\lambda\mathbf{p})(\mathbf{r}-\lambda\mathbf{p})}{|\mathbf{r}-\lambda\mathbf{p}|^{3}}\Big{\\}},$
(53)

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

$\displaystyle\boldsymbol{h^{(0)}}(\lambda;\boldsymbol{\omega}^{(0)}_{p},\mathbf{p})=-4\pi\lambda[\mathbf{p}\cdot(\nabla\mathbf{u})_{\infty}-\boldsymbol{\omega}^{(0)}_{p}\times\mathbf{p}]\cdot\Big{(}\frac{1}{2}\mathbf{pp}+(\mathbf{I}-\mathbf{p}\mathbf{p})\frac{\log(2\kappa)-1.5}{\log(2\kappa)-0.5}\Big{)},$
$\displaystyle\mathbf{f}^{(0)}(\lambda;\mathbf{p})=-4\pi\lambda\begin{bmatrix}0&p_{1}&0\end{bmatrix}^{T}\Big{(}1-\frac{1}{\log(2\kappa)-0.5}\Big{)}+\mathcal{O}(p_{2}),$
(54)
$\displaystyle\boldsymbol{\mathcal{G}}^{(0)}(\lambda;\boldsymbol{\omega}^{(0)}_{p},\mathbf{p})=2\pi(1-\lambda^{2})\Bigg{(}\begin{bmatrix}0&1+\frac{p_{3}^{2}}{2}&0\\\
\frac{p_{3}^{2}}{2}&0&-\frac{1}{2}p_{1}p_{3}\\\
0&-\frac{1}{2}p_{1}p_{3}&0\end{bmatrix}+(\boldsymbol{\omega}^{(0)}_{p}\cdot\mathbf{p})\boldsymbol{\epsilon}\cdot\mathbf{p}\Bigg{)}+\mathcal{O}(p_{2}).$

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, 1922; Kim & Karrila, 2013),

$\dot{\mathbf{p}}^{(0)}=\boldsymbol{\omega}^{(0)}_{p}\times\mathbf{p}=\boldsymbol{\omega}_{\infty}\times\mathbf{p}+\frac{\kappa^{2}+1}{\kappa^{2}-1}(\mathbf{E}_{\infty}\cdot\mathbf{p})\cdot(\mathbf{I}-\mathbf{pp}),$
(55)

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

$\boldsymbol{\omega}_{\infty}=-\frac{1}{2}\begin{bmatrix}0&0&1\end{bmatrix}^{T},\hskip
14.45377pt\mathbf{E}_{\infty}=\frac{1}{2}\begin{bmatrix}0&1&0\\\ 1&0&0\\\
0&0&0\end{bmatrix}.$ (56)

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

$\displaystyle\boldsymbol{\omega}^{(0)}_{p}\cdot\mathbf{p}=-\frac{p_{3}}{2}+\mathcal{O}(p_{2}).$
(57)

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}\cdot\mathbf{p}$ from equation (57) and
$\boldsymbol{\omega}^{(0)}_{p}\times\mathbf{p}$ from equation (55) are used to
obtain $\boldsymbol{\mathcal{G}}^{(0)}$ and $\boldsymbol{h^{(0)}}$,
respectively, in equation (54).

The flow due to $\boldsymbol{h}^{(0)}$ is taken from Batchelor (1970). 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 equation (53)
and (4.1). The flow generated at the next order in $\kappa$ arises at
$\mathcal{O}(\kappa^{-2})$ and is generated by $\mathbf{f}^{(0)}$ and
$\boldsymbol{\mathcal{G}}^{(0)}$. It is taken from Cox (1971) where only the
$\mathcal{O}(\kappa^{-2})$ flow is available. Thus, the disturbance velocity
in the outer region given by equation (53) has an overall error of
$\mathcal{O}(\kappa^{-3})$. If only $\boldsymbol{h}^{(0)}$ is considered, the
error is of $\mathcal{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 $\mathbf{f}^{(0)}$ and
$\boldsymbol{\mathcal{G}}^{(0)}$ capture the (highest)
$\mathcal{O}(\kappa^{-2})$ disturbance created by the particle. Accounting for
the flow generated by $\mathbf{f}^{(0)}$ and $\boldsymbol{\mathcal{G}}^{(0)}$
allows us to consider the influence of elasticity within the flow-vorticity
plane.

As mentioned just before equation (50) 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}|\leq\mathcal{O}(1/\kappa)$.
Therefore, in $\boldsymbol{h}^{(0)}$ we only consider the flow at
$\mathcal{O}(p_{2})$ i.e. a flow of $\mathcal{O}(p_{2}/\log(\kappa))$ with an
error of $\mathcal{O}(p_{2}^{2}/\log(\kappa))$. Since the primary purpose of
using the $\mathbf{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 $\mathbf{f}^{(0)}$ and
$\boldsymbol{\mathcal{G}}^{(0)}$ with $p_{2}=0$ (as expressed in equation
(54)) which leads to a flow of $\mathcal{O}(1/\kappa^{2})$ with an error of
$\mathcal{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}|\leq\mathcal{O}(1/\kappa)$ regime considered, the errors in
$\boldsymbol{h}^{(0)}$ driven flow, i.e.,
$\mathcal{O}(p_{2}^{2}/\log(\kappa))$, are always smaller than the flow
generated by $\mathbf{f}^{(0)}$ and $\boldsymbol{\mathcal{G}}^{(0)}$, i.e.,
$\mathcal{O}(1/\kappa^{2})$. As $p_{2}$ approaches zero, the
$\boldsymbol{h}^{(0)}$ driven flow and the associated errors fall rapidly to
zero making $\mathbf{f}^{(0)}$ and $\boldsymbol{\mathcal{G}}^{(0)}$ driven
flow at $\mathcal{O}(1/\kappa^{2})$ the dominating one. Hence, the error terms
arising from either of the $\boldsymbol{h}^{(0)}$ or $\mathbf{f}^{(0)}$ and
$\boldsymbol{\mathcal{G}}^{(0)}$ driven flow are always lower than the actual
flow when $|p_{2}|\leq\mathcal{O}(1/\kappa)$. For
$|p_{2}|>\mathcal{O}(1/\kappa)$, the Newtonian rotation rate of the fiber
dominates over the changes due to elasticity, and we consider the exact
Jeffery (1922) rotation to account for it.

We have considered spheroidal particles in our theory. However, in experiments
with slender particles (Gauthier et al., 1971; Bartram et al., 1975; Iso et
al., 1996b, a) it is convenient to fabricate cylindrical particles. The forces
generated by the blunt ends of a slender cylinder lead to an additional torque
of $\mathcal{O}(1/\kappa^{2})$ (Cox, 1971) rendering the
$\mathcal{O}(1/\kappa^{2})$ flow generated by $\mathbf{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, 1970). Taking the
$\mathcal{O}(p_{2})$ terms with this velocity field can still allow us to
consider the flow accurately up to $\mathcal{O}(p_{2}/\log(\kappa))$. Thus the
upcoming theoretical development for the orientation dynamics of a slender
prolate spheroid leading up to equation (95) can still be used for a slender
cylinder while ignoring the $\mathbf{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 section 5 while
analyzing 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 $\mathcal{O}(c)$ is due to a fiber rotating with the
perturbed angular velocity, $\boldsymbol{\omega}^{(1)}_{p}$, in a quiescent
fluid, i.e., equations (48)-(49). In the outer region, the fluid velocity is,

$\displaystyle{{\mathbf{u}_{\text{out}}^{\text{M}}}^{(1)}}$
$\displaystyle(\mathbf{r};\boldsymbol{\omega}^{(1)}_{p},\mathbf{p})=\frac{1}{8\pi}\int_{-1}^{1}d\lambda\hskip
7.22743pt\frac{\boldsymbol{h}^{(1)}(\lambda;\boldsymbol{\omega}^{(1)}_{p},\mathbf{p})}{\log(2\kappa)-0.5}\cdot\Big{\\{}\frac{\mathbf{I}}{|\mathbf{r}-\lambda\mathbf{p}|}+\frac{(\mathbf{r}-\lambda\mathbf{p})(\mathbf{r}-\lambda\mathbf{p})}{|\mathbf{r}-\lambda\mathbf{p}|^{3}}\Big{\\}}+$
$\displaystyle\frac{1}{8\pi}\int_{-1}^{1}d\lambda\hskip
7.22743pt\frac{1}{\kappa^{2}}\boldsymbol{\mathcal{G}}^{(1)}(\lambda;\boldsymbol{\omega}^{(1)}_{p},\mathbf{p}):\nabla\Big{\\{}\frac{\mathbf{I}}{|\mathbf{r}-\lambda\mathbf{p}|}+\frac{(\mathbf{r}-\lambda\mathbf{p})(\mathbf{r}-\lambda\mathbf{p})}{|\mathbf{r}-\lambda\mathbf{p}|^{3}}\Big{\\}},$
(58)

where,

$\boldsymbol{h^{(1)}}(\lambda;\boldsymbol{\omega}^{(1)}_{p},\mathbf{p})=4\pi\lambda\boldsymbol{\omega}^{(1)}_{p}\times\mathbf{p},\hskip
14.45377pt\boldsymbol{\mathcal{G}}^{(1)}(\lambda;\boldsymbol{\omega}^{(1)}_{p},\mathbf{p})=-2\pi(1-\lambda^{2})(\boldsymbol{\omega}^{(1)}_{p}\cdot\mathbf{p})\boldsymbol{\epsilon}\cdot\mathbf{p}.$
(59)

The torque generated by this flow is,

$\displaystyle\begin{split}&{\mathbf{G}^{\text{MIST}}}^{(1)}(\boldsymbol{\omega}^{(1)}_{p},\mathbf{p})=\\\
&\int_{-1}^{1}d\lambda\hskip
7.22743pt\Big{(}\frac{1}{\kappa^{2}}\boldsymbol{\mathbf{\epsilon}}:\boldsymbol{\mathcal{G}}^{(1)}-\frac{\lambda\mathbf{p}\times\boldsymbol{h}^{(1)}}{\log(2\kappa)-0.5}\Big{)}=\frac{8\pi}{3}\Big{(}\frac{(\boldsymbol{\omega}^{(1)}_{p}\cdot\mathbf{p})\mathbf{p}}{\kappa^{2}}-\frac{\boldsymbol{\omega}^{(1)}_{p}\cdot(\mathbf{I}-\mathbf{pp})}{\log(2\kappa)-0.5}\Big{)}\end{split}$
(60)

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

$\dot{\mathbf{p}}^{(1)}=\boldsymbol{\omega}^{(1)}_{p}\times\mathbf{p}=-\frac{3}{8\pi}({\log(2\kappa)-0.5})({\mathbf{G}^{\text{MIST}}}^{(1)}(\boldsymbol{\omega}^{(1)}_{p},\mathbf{p})\times\mathbf{p}),$
(61)

Using the torque balance at $\mathcal{O}(c)$ from equation (47) we obtain,

$\dot{\mathbf{p}}^{(1)}=\dot{\mathbf{p}}^{(1)}_{\text{PIST}}+\dot{\mathbf{p}}^{(1)}_{\text{Elastic}},$
(62)

where

$\dot{\mathbf{p}}^{(1)}_{\text{PIST}}=\frac{3}{8\pi}({\log(2\kappa)-0.5})({\mathbf{G}^{\text{PIST}}}^{(1)}\times\mathbf{p}),$
(63)

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

$\dot{\mathbf{p}}^{(1)}_{\text{Elastic}}=\frac{3}{8\pi}({\log(2\kappa)-0.5})({\mathbf{G}^{\text{Elastic}}}^{(1)}\times\mathbf{p}),$
(64)

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

###  4.2 Rotation due to polymer induced solvent stress

From equations (4), (44), (46) and (63), the rotation rate due to the polymer
induced solvent stress is,

$\dot{\mathbf{p}}^{(1)}_{\text{PIST}}=-\frac{3}{8\pi
De}({\log(2\kappa)-0.5})\mathbf{p}\times\int_{\text{Fluid}}dV\hskip
7.22743pt(\nabla\cdot({\boldsymbol{\Lambda}^{(0)}}-{\boldsymbol{\Lambda}^{U}}))\cdot\mathbf{b}.$
(65)

We remind the reader that here $\mathbf{b}$ is the auxillary ‘velocity’ field
used in the reciprocal theorem for deriving the polymer induced solvent torque
(equation (24) or (46)). $\mathbf{b}\cdot\boldsymbol{\omega}$ corresponds to
the fluid velocity around a fiber rotating with an angular velocity
$\boldsymbol{\omega}$ in a quiescent fluid. Therefore,

$\mathbf{b}=\nabla_{\boldsymbol{\omega}^{(1)}_{p}}{\mathbf{u}^{\text{M}}}^{(1)},$
(66)

where ${\mathbf{u}^{\text{M}}}^{(1)}$ is the solution to equations (48) and
(49). The volume integral in equation (65) 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,

$\dot{\mathbf{p}}^{(1)}_{\text{PIST}}\approx-\frac{3}{8\pi
De}({\log(2\kappa)-0.5})\mathbf{p}\times\int d^{3}\mathbf{r}\hskip
7.22743pt\nabla\cdot({\boldsymbol{\Lambda}_{\text{out}}^{(0)}}-{\boldsymbol{\Lambda}^{U}})\cdot\mathbf{b}_{\text{out}},$
(67)

where the volume integral is taken over the entire space,
${\boldsymbol{\Lambda}_{\text{out}}^{(0)}}$ is the polymer conformation in the
outer region and from equation (58), (59) and (66),

$\displaystyle\mathbf{b}_{\text{out}}=$
$\displaystyle\frac{1}{8\pi}\int_{-1}^{1}d\lambda\hskip
7.22743pt\frac{4\pi\lambda\mathbf{p}\cdot\boldsymbol{\epsilon}}{\log(2\kappa)-0.5}\cdot\Big{(}\frac{\mathbf{I}}{|\mathbf{r}-\lambda\mathbf{p}|}+\frac{(\mathbf{r}-\lambda\mathbf{p})(\mathbf{r}-\lambda\mathbf{p})}{|\mathbf{r}-\lambda\mathbf{p}|^{3}}\Big{)}+\mathcal{O}(\kappa^{-2}).$
(68)

It is straightforward to show,

$\dot{\mathbf{p}}^{(1)}_{\text{PIST}}\approx\frac{3}{2}(\mathbf{I}-\mathbf{p}\mathbf{p})\cdot\int_{-1}^{1}d\lambda\hskip
7.22743pt\lambda\mathbf{u}_{\text{out}}^{\text{polymer}}(\lambda\mathbf{p}),$
(69)

where,

$\mathbf{u}_{\text{out}}^{\text{polymer}}(\lambda\mathbf{p})=\frac{1}{8\pi}\int
d^{3}\mathbf{r}\hskip
7.22743pt\frac{1}{De}\nabla\cdot({\boldsymbol{\Lambda}_{\text{out}}^{(0)}}-{\boldsymbol{\Lambda}^{U}})\cdot\Big{(}\frac{\mathbf{I}}{|\mathbf{r}-\lambda\mathbf{p}|}+\frac{(\mathbf{r}-\lambda\mathbf{p})(\mathbf{r}-\lambda\mathbf{p})}{|\mathbf{r}-\lambda\mathbf{p}|^{3}}\Big{)},$
(70)

is the velocity field evaluated on the particle center-line that is produced
by polymeric force,
$\frac{1}{De}\nabla\cdot({\boldsymbol{\Lambda}_{\text{out}}^{(0)}}-{\boldsymbol{\Lambda}^{U}})$,
acting in the outer region. This form of the rotation rate, equations (69) and
(70), was considered by (Harlen & Koch, 1993). However, they did not account
for all the relevant terms in calculating
${\boldsymbol{\Lambda}_{\text{out}}^{(0)}}$ as we show below.

As discussed in section 3 and shown by equations (43) and (50), the polymer
conformation, ${\boldsymbol{\Lambda}^{(0)}}$ depends on the leading order
Newtonian velocity. In the outer region, the velocity disturbance created by
the particle is $\mathcal{O}(\max[p_{2}/\log(\kappa),1/\kappa^{2}])$ smaller
than the velocity of the imposed simple shear (equations (52), (53) and (54)).
Thus, we linearize equation (50) in the outer region about the imposed flow
field, $\mathbf{r}\cdot(\nabla\mathbf{u})_{\infty}$, and obtain the governing
equation for the disturbance in the polymer conformation from its undisturbed
value,

${\boldsymbol{\Lambda}_{\text{out}}^{(0)}}^{\prime}=\boldsymbol{\Lambda}_{\text{out}}^{(0)}-\boldsymbol{\Lambda}_{\text{out}}^{U},$
(71)

to be

$\displaystyle\begin{split}\Big{(}\mathbf{r}\cdot(\nabla\mathbf{u})_{\infty}\cdot\nabla+\frac{1}{De}\Big{)}{\boldsymbol{\Lambda}_{\text{out}}^{(0)}}^{\prime}&-(\nabla\mathbf{u})_{\infty}^{T}\cdot{\boldsymbol{\Lambda}_{\text{out}}^{(0)}}^{\prime}-{\boldsymbol{\Lambda}_{\text{out}}^{(0)}}^{\prime}\cdot(\nabla\mathbf{u})_{\infty}=\\\
&(\nabla{\mathbf{u}_{\text{out}}^{\text{M}{(0)}}}^{\prime})^{T}\cdot\boldsymbol{\Lambda}^{U}+\boldsymbol{\Lambda}^{U}\cdot\nabla{\mathbf{u}_{\text{out}}^{\text{M}{(0)}}}^{\prime}.\end{split}$
(72)

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

$\displaystyle\begin{split}\Big{(}-k_{1}\frac{\partial}{\partial
k_{2}}+\frac{1}{De}\Big{)}{\hat{\boldsymbol{\Lambda}}_{\text{out}}^{(0)}}^{\prime}&-(\nabla\mathbf{u})_{\infty}^{T}\cdot{\hat{\boldsymbol{\Lambda}}_{\text{out}}^{(0)}}^{\prime}-{\hat{\boldsymbol{\Lambda}}_{\text{out}}^{(0)}}^{\prime}\cdot(\nabla\mathbf{u})_{\infty}\\\
&=i({{\hat{\mathbf{u}}_{\text{out}}^{\text{M}}}^{(0)}}^{\prime}\mathbf{k}\cdot\boldsymbol{\Lambda}^{U}+\boldsymbol{\Lambda}^{U}\cdot\mathbf{k}{{\hat{\mathbf{u}}_{\text{out}}^{\text{M}}}^{(0)}}^{\prime}),\end{split}$
(73)

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

$\dot{\mathbf{p}}^{(1)}_{\text{PIST}}\approx\frac{3}{2}(\mathbf{I}-\mathbf{p}\mathbf{p})\cdot\int_{-1}^{1}d\lambda\int
d^{3}\mathbf{k}\hskip 7.22743pt\lambda\exp(2\pi
i\lambda\mathbf{k}\cdot\mathbf{p})\hat{\mathbf{u}}_{\text{out}}^{\text{polymer}},$
(74)

where from the convolution theorem,

$\hat{\mathbf{u}}_{\text{out}}^{\text{polymer}}=i\frac{1}{De}\mathbf{k}\cdot{\hat{\boldsymbol{\Lambda}}_{\text{out}}^{(0)}}^{\prime}(\mathbf{k};\mathbf{p})\cdot\hat{\mathbf{J}}(\mathbf{k}),\text{
with
}\hat{\mathbf{J}}(\mathbf{k})=\mathcal{F}\Big{(}\frac{1}{8\pi}\Big{(}\frac{\mathbf{I}}{r}-\frac{\mathbf{rr}}{r^{3}}\Big{)}\Big{)}=\frac{1}{k^{2}}\Big{(}\mathbf{I}-\frac{\mathbf{kk}}{k^{2}}\Big{)}.$
(75)

We solve the polymer constitutive equation (73) in Fourier space using the
method of characteristics to obtain
${\hat{\boldsymbol{\Lambda}}_{\text{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\text{sgn}(k_{1})$ and $k_{2}$ and the
boundary condition is
${\hat{\boldsymbol{\Lambda}}_{\text{out}}^{(0)}}^{\prime}(\infty\text{sgn}(k_{1}))=0$.
We use computer algebra for this calculation.
${\hat{\boldsymbol{\Lambda}}_{\text{out}}^{(0)}}^{\prime}$ hence obtained
allows us to evaluate the integral in equation (74). The expression used for
the disturbance in the Newtonian velocity field,
${\mathbf{u}_{\text{out}}^{\text{M}{(0)}}}^{\prime}$, has errors of
$\mathcal{O}(\text{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 $\mathcal{O}(\epsilon_{\dot{p}})$, where,

$\epsilon_{\dot{p}}=\text{max}\Bigg{[}\frac{p_{2}^{2}}{\log(\kappa)},\frac{p_{2}}{\kappa^{2}},\frac{1}{\kappa^{3}}\Bigg{]}.$
(76)

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

$\dot{{p}}^{(1)}_{\text{2,PIST}}=-\frac{Dep_{1}^{2}p_{2}}{2\log(2\kappa)-3}+\mathcal{O}(\epsilon_{\dot{p}}),$
(77)

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

${\mathbf{G}^{\text{PIST}}}^{(1)}\approx\frac{1}{De}\int d^{3}\mathbf{r}\hskip
7.22743pt\nabla\cdot({\boldsymbol{\Lambda}_{\text{out}}^{(0)}}-{\boldsymbol{\Lambda}^{U}})\cdot\mathbf{b}_{\text{out}},$
(78)

we obtain the 1st and 3rd components of the polymer induced solvent torque,

$\frac{{{G}_{1}^{\text{PIST}}}^{(1)}}{p_{3}}=-\frac{{{G}_{3}^{\text{PIST}}}^{(1)}}{p_{1}}=\frac{4De^{2}p_{1}^{2}p_{2}}{3(\log(2\kappa)-0.5)(\log(2\kappa)-1.5)}+\mathcal{O}\Bigg{(}\frac{\epsilon_{\dot{p}}}{\log(\kappa)}\Bigg{)}.$
(79)

Since,
$\dot{\mathbf{p}}^{(1)}_{\text{PIST}}={3}/{(8\pi)}({\log(2\kappa)-0.5})({\mathbf{G}^{\text{PIST}}}^{(1)}\times\mathbf{p})$
(equation (63)) we may obtain
$\dot{{p}}^{(1)}_{\text{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}^{\text{PIST}}}^{(1)}$ or rotation rates in flow and
vorticity direction, $\dot{{p}}^{(1)}_{\text{1,PIST}}$ and
$\dot{{p}}^{(3)}_{\text{1,PIST}}$ are not tractable for a general $De$.
Therefore, to obtain these components we consider small and large $De$ limits
separately in sections 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)}_{\text{2, PIST}}$ to be dependent on the
first normal stress difference of the Oldroyd-B fluid. This is because
$\dot{\mathbf{p}}^{(1)}_{\text{PIST}}$ in equation (74) is directly
proportional to ${\hat{\boldsymbol{\Lambda}}_{\text{out}}^{(0)}}^{\prime}$
that is in turn dependent on $\boldsymbol{\Lambda}^{U}$ via equation (73).
This first normal stress difference dependence is further confirmed by
repeating the above calculation by artificially setting
${\Lambda}^{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 section
4.2.1 to arise from the first normal stress difference of the Oldroyd-B fluid.
The rotation rate for a second order fluid ($\mathcal{O}(De)$ rotation rate in
the small $De$ calculation of section 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 modeling 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., 1987) of simplicity and hence better
analytical tractability. As mentioned in section 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$, equation
(73) simplifies to,

$-k_{1}\frac{\partial}{\partial
k_{2}}{\hat{\boldsymbol{\Lambda}}_{\text{out}}^{(0)}}^{\prime}-(\nabla\mathbf{u})_{\infty}^{T}\cdot{\hat{\boldsymbol{\Lambda}}_{\text{out}}^{(0)}}^{\prime}-{\hat{\boldsymbol{\Lambda}}_{\text{out}}^{(0)}}^{\prime}\cdot(\nabla\mathbf{u})_{\infty}=i({{\hat{\mathbf{u}}_{\text{out}}^{\text{M}}}^{(0)}}^{\prime}\mathbf{k}\cdot\boldsymbol{\Lambda}^{U}+\boldsymbol{\Lambda}^{U}\cdot\mathbf{k}{{\hat{\mathbf{u}}_{\text{out}}^{\text{M}}}^{(0)}}^{\prime}).$
(80)

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

${\lim\limits_{De\gg
1}\dot{\mathbf{p}}^{(1)}_{\text{PIST}}=De\begin{bmatrix}-p_{1}p_{3}^{2}\alpha_{De\gg
1}\\\ -\frac{p_{1}^{2}p_{2}}{2\log(2\kappa)-3}\\\ p_{1}^{2}p_{3}\alpha_{De\gg
1}\end{bmatrix},\text{ with }\alpha_{De\gg
1}=\frac{-4}{\kappa^{2}}+\frac{3\pi}{4\log(2\kappa)-6}p_{2}\tan^{-1}\frac{p_{1}}{|p_{3}|}.}$
(81)

The large $De$ approximation can be viewed as the leading $\mathcal{O}(De)$
term in an expansion in powers of $1/De$. Thus the error due to the expansion
in $1/De$ is $\mathcal{O}(1)$ in the 1 and 3 components of $\lim\limits_{De\gg
1}\dot{\mathbf{p}}^{(1)}_{\text{PIST}}$. There is no error due to expansion in
$1/De$ for $\dot{{p}}^{(1)}_{\text{2,PIST}}$ as same value is obtained via
$De\gg 1$ approximation as that for a general $De$ in equation (77). As
mentioned earlier, $\dot{{p}}^{(1)}_{\text{2,PIST}}$ has an error of
$\mathcal{O}(\epsilon_{\dot{p}})$ (equation (76)).
$\dot{{p}}^{(1)}_{\text{1,PIST}}$ and $\dot{{p}}^{(1)}_{\text{3,PIST}}$ have
an error of $\mathcal{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
(1993). However they did not account for the stretching of the conformation
disturbance, ${\hat{\boldsymbol{\Lambda}}_{\text{out}}^{(0)}}^{\prime}$, by
the mean velocity gradients i.e. they omitted
$(\nabla\mathbf{u})_{\infty}^{T}\cdot{\hat{\boldsymbol{\Lambda}}_{\text{out}}^{(0)}}^{\prime}+{\hat{\boldsymbol{\Lambda}}_{\text{out}}^{(0)}}^{\prime}\cdot(\nabla\mathbf{u})_{\infty}$
term in equation (80).

####  4.2.2 Small $De$

When $De\ll 1$, a solution of equation (72) or (73) is obtained via a regular
perturbation expansion of
${\hat{\boldsymbol{\Lambda}}_{\text{out}}^{(0)}}^{\prime}$ in the powers of
$De$

${{\hat{\boldsymbol{\Lambda}}_{\text{out}}^{(0)}}^{\prime}}=\Sigma_{n=0}^{\infty}De^{n}{{\hat{\boldsymbol{\Lambda}}_{\text{out}}^{(0)}}^{\prime}}^{(n)}.$
(82)

The first two terms in this expansion are

$\displaystyle\begin{split}&{{\hat{\boldsymbol{\Lambda}}_{\text{out}}^{(0)}}^{\prime}}^{(0)}=\mathbf{0},\\\
&{{\hat{\boldsymbol{\Lambda}}_{\text{out}}^{(0)}}^{\prime}}^{(1)}=i(\mathbf{k}{\hat{\mathbf{u}}^{\text{M}}}^{(0)}+{\hat{\mathbf{u}}^{\text{M}}}^{(0)}\mathbf{k}).\end{split}$
(83)

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

$\displaystyle\begin{split}{\hat{\boldsymbol{\Lambda}}_{\text{out}}^{(0^{\prime})}}^{(2)}=&(\nabla\mathbf{u})_{\infty}^{T}\cdot{\hat{\boldsymbol{\Lambda}}_{\text{out}}^{(0^{\prime})}}^{(1)}+{\hat{\boldsymbol{\Lambda}}_{\text{out}}^{(0^{\prime})}}^{(1)}\cdot(\nabla\mathbf{u})_{\infty}\\\
&+i({\hat{\mathbf{u}}^{\text{M}}}^{(0)}\mathbf{k}\cdot{\boldsymbol{\Lambda}^{U}}^{(1)}+{\boldsymbol{\Lambda}^{U}}^{(1)}\cdot\mathbf{k}{\hat{\mathbf{u}}^{\text{M}}}^{(0)})+k_{1}\frac{\partial}{\partial
k_{2}}{\hat{\boldsymbol{\Lambda}}_{\text{out}}^{(0^{\prime})}}^{(1)},\\\
{\hat{\boldsymbol{\Lambda}}_{\text{out}}^{(0^{\prime})}}^{(3)}=&(\nabla\mathbf{u})_{\infty}^{T}\cdot{\hat{\boldsymbol{\Lambda}}_{\text{out}}^{(0^{\prime})}}^{(2)}+{\hat{\boldsymbol{\Lambda}}_{\text{out}}^{(0^{\prime})}}^{(2)}\cdot(\nabla\mathbf{u})_{\infty}\\\
&+i({\hat{\mathbf{u}}^{\text{M}}}^{(0)}\mathbf{k}\cdot{\boldsymbol{\Lambda}^{U}}^{(2)}+{\boldsymbol{\Lambda}^{U}}^{(2)}\cdot\mathbf{k}{\hat{\mathbf{u}}^{\text{M}}}^{(0)})+k_{1}\frac{\partial}{\partial
k_{2}}{\hat{\boldsymbol{\Lambda}}_{\text{out}}^{(0^{\prime})}}^{(2)},\\\
{\hat{\boldsymbol{\Lambda}}_{\text{out}}^{(0^{\prime})}}^{(n)}=&(\nabla\mathbf{u})_{\infty}^{T}\cdot{\hat{\boldsymbol{\Lambda}}_{\text{out}}^{(0^{\prime})}}^{(n-1)}+{\hat{\boldsymbol{\Lambda}}_{\text{out}}^{(0^{\prime})}}^{(n-1)}\cdot(\nabla\mathbf{u})_{\infty}+k_{1}\frac{\partial}{\partial
k_{2}}{\hat{\boldsymbol{\Lambda}}_{\text{out}}^{(0^{\prime})}}^{(n-1)},n\geq
4.\end{split}$ (84)

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

$\displaystyle\begin{split}\lim\limits_{De\ll
1}\dot{\mathbf{p}}^{(1)}_{\text{PIST}}=&-\frac{3}{8\pi}({\log(2\kappa)-0.5})\mathbf{p}\times\Sigma_{n=1}^{\infty}De^{n-1}\int_{\text{Fluid}}dV\hskip
7.22743pt(\nabla\cdot({{{\boldsymbol{\Lambda}}_{\text{out}}^{(0)}}^{\prime}}^{(n)}))\cdot\mathbf{b}.\end{split}$
(85)

We will discuss the contribution of the $n=1$ term, i.e. $\mathcal{O}(1)$ term
in $De$ expansion of $\dot{\mathbf{p}}^{(1)}_{\text{PIST}}$ in the above
equation in section 4.4, along with the similar term in the expansion of
$\dot{\mathbf{p}}^{(1)}_{\text{Elastic}}$. Higher order contributions to
$\dot{\mathbf{p}}^{(1)}_{\text{PIST}}$ can be evaluated in Fourier space. From
equations (82)-(84), using expansion of
${{\boldsymbol{\Lambda}}_{\text{out}}^{(0)}}^{\prime}$ up to
$\mathcal{O}(De^{8})$ (even higher orders can be easily evaluated), we obtain,

$\displaystyle\begin{split}\lim\limits_{De\ll
1}\dot{\mathbf{p}}^{(1)}_{\text{PIST}}=&-\frac{3}{8\pi}({\log(2\kappa)-0.5})\mathbf{p}\times\int_{\text{Fluid}}dV\hskip
7.22743pt(\nabla\cdot({{{\boldsymbol{\Lambda}}_{\text{out}}^{(0)}}^{\prime}}^{(1)}))\cdot\mathbf{b}\\\
&+De\begin{bmatrix}-p_{1}p_{3}^{2}\alpha_{De\ll 1}\\\
-\frac{p_{1}^{2}p_{2}}{2\log(2\kappa)-3}\\\ p_{1}^{2}p_{3}\alpha_{De\ll
1}\end{bmatrix},\end{split}$ (86)

where,

$\alpha_{De\ll
1}=-\frac{4}{\kappa^{2}}+\frac{p_{1}p_{2}}{8\log(2\kappa)-12}\Big{(}5De-\frac{3}{2}(2+p_{3}^{2})De^{3}+\frac{13}{10}(8+p_{3}^{2}(4+3p_{3}^{2}))De^{5}\Big{)}.$
(87)

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

###  4.3 Rotation due to elastic stress

The elastic torque, ${\mathbf{G}^{\text{Elastic}}}^{(1)}$, is due to the
polymer stress on the particle surface (equations (44) and (46)). 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 quasi-steady constitutive equation (50) simplifies to a set of
coupled algebraic equations,

${(\nabla{\mathbf{u}_{\text{in}}^{\text{M}}}^{(0)}(\mathbf{r}_{p};\boldsymbol{\omega}^{(0)}_{p},\mathbf{p}))^{T}\cdot\boldsymbol{\Lambda}^{(0)}_{p}+\boldsymbol{\Lambda}^{(0)}_{p}\cdot\nabla{\mathbf{u}_{\text{in}}^{\text{M}}}^{(0)}(\mathbf{r}_{p};\boldsymbol{\omega}^{(0)}_{p},\mathbf{p})+\frac{1}{De}(\mathbf{I}-\boldsymbol{\Lambda}^{(0)}_{p})=0,}$
(88)

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

$\displaystyle\begin{split}{\mathbf{G}}^{(1)}_{\text{Elastic}}&=\frac{8\pi
c}{3\kappa^{2}}\begin{bmatrix}p_{1}p_{3}\\\ 0\\\
-p_{1}^{2}\end{bmatrix}+\mathcal{O}(\max[p_{2}\kappa^{-2},\kappa^{-3}]).\end{split}$
(89)

and the corresponding rotation rate,

$\displaystyle\begin{split}\dot{\mathbf{p}}^{(1)}_{\text{Elastic}}&=-\frac{3}{8\pi}({\log(2\kappa)-0.5})\mathbf{p}\times\int
dS\hskip
7.22743pt\Big{(}\mathbf{r}\times\frac{1}{De}(\boldsymbol{\Lambda}^{(0)}_{p}-\mathbf{I})\cdot\mathbf{n}\Big{)}\\\
&=\frac{{\log(2\kappa)-0.5}}{\kappa^{2}}p_{1}\begin{bmatrix}0\\\ 1\\\
0\end{bmatrix}+\mathcal{O}(\log(\kappa)\max[p_{2}\kappa^{-2},\kappa^{-3}]).\end{split}$
(90)

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

###  4.4 Net rotation due to the polymers

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

For $De\ll 1$, we can obtain the $\mathcal{O}(1)$ term in the $De$ expansion
of $\dot{\mathbf{p}}^{(1)}_{\text{PIST}}$ and
$\dot{\mathbf{p}}^{(1)}_{\text{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 equation (83), i.e.,

${{{{\boldsymbol{\Lambda}}^{(0)}}}=\mathbf{I}+De(\nabla{\mathbf{u}^{\text{M}}}^{(0)}+(\nabla{\mathbf{u}^{\text{M}}}^{(0)})^{T})+\mathcal{O}(De^{2}).}$
(91)

From the combined effect of polymers from equations (46), (62), (63) and (64),
along with using the value of ${{{\boldsymbol{\Lambda}}^{(0)}}}$ from equation
(91), divergence theorem, the leading order momentum equation (i.e.
$\nabla{p^{\text{P}}}^{(0)}=\nabla\cdot[\nabla{\mathbf{u}^{\text{P}}}^{(1)}+(\nabla{\mathbf{u}^{\text{P}}}^{(1)})^{T}]$),
and the auxillary Stokes problem (equations (18) and (19)), we obtain,

$\displaystyle\begin{split}&\lim\limits_{De\ll
1}{\dot{\mathbf{p}}^{(1)}}_{\text{Elastic}}+{\dot{\mathbf{p}}^{(1)}}_{\text{PIST}}\\\
&=\int_{\mathbf{r}_{\text{p}}}dS\hskip
7.22743pt\mathbf{r}\times[(\nabla{\mathbf{u}^{\text{M}}}^{(0)}+(\nabla{\mathbf{u}^{\text{M}}}^{(0)})^{T}-{p^{\text{M}}}^{(0)}\mathbf{I})\cdot\mathbf{n}]+\mathcal{O}(De)\\\
&=\int_{\mathbf{r}_{\text{p}}}dS\hskip
7.22743pt\mathbf{r}\times{\boldsymbol{\tau}^{\text{M}}}^{(0)}+\mathcal{O}(De)\cdot\mathbf{n}=\mathcal{O}(De).\end{split}$
(92)

The $\mathcal{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 $\mathcal{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 $\mathcal{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 equation (90), which we have just shown balances the equivalent
rotation rate due to polymer induced solvent torque for $De\ll 1$. Thus the
$\mathcal{O}(De)$ and higher effect of polymers on the particle’s rotation
rate arises entirely due to the polymer induced solvent torque (equation (86))
and is given by,

$\lim\limits_{De\ll
1}(\dot{\mathbf{p}}^{(1)}_{\text{Elastic}}+\dot{\mathbf{p}}^{(1)}_{\text{PIST}})=De\begin{bmatrix}-p_{1}p_{3}^{2}\alpha_{De\ll
1}\\\ -\frac{p_{1}^{2}p_{2}}{2\log(2\kappa)-3}\\\ p_{1}^{2}p_{3}\alpha_{De\ll
1}\end{bmatrix},$ (93)

where, $\alpha_{De\ll 1}$ is given in equation (87) and the errors are
$\mathcal{O}(De\epsilon_{\dot{p}})$ for all the components (where
$\epsilon_{\dot{p}}$ is given in equation (76)).

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
$\mathcal{O}(De)$ smaller. The net rotation rate due to the polymers is,

$\lim\limits_{De\gg
1}(\dot{\mathbf{p}}^{(1)}_{\text{Elastic}}+\dot{\mathbf{p}}^{(1)}_{\text{PIST}})=De\begin{bmatrix}-p_{1}p_{3}^{2}\alpha_{De\gg
1}\\\ -\frac{p_{1}^{2}p_{2}}{2\log(2\kappa)-3}\\\ p_{1}^{2}p_{3}\alpha_{De\gg
1}\end{bmatrix}.$ (94)

where $\alpha_{De\gg 1}$ is given by equation (81). The error in
$\dot{{p}}^{(1)}_{\text{1,PIST}}$ and $\dot{{p}}^{(1)}_{\text{3,PIST}}$ is of
$\mathcal{O}(De\epsilon_{\dot{p}})$ (equation (76)). Neglecting the elastic
torque leads to an additional error of $\mathcal{O}(\log(\kappa)/\kappa^{2})$
in the second component and
$\mathcal{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 (equation (90)).

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

Equations (93) and (94) 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 ($\mathcal{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 section 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,

$\displaystyle\dot{\mathbf{p}}=\frac{d\mathbf{p}}{dt}=\begin{bmatrix}\frac{p_{2}}{2}+\frac{\kappa^{2}-1}{\kappa^{2}+1}\Big{(}\frac{p_{2}}{2}-p_{1}^{2}p_{2}\Big{)}-cDep_{1}p_{3}^{2}\alpha\\\
-\frac{p_{1}}{2}+\frac{\kappa^{2}-1}{\kappa^{2}+1}\Big{(}\frac{p_{1}}{2}-p_{1}p_{2}^{2}\Big{)}-cDe\frac{p_{1}^{2}p_{2}}{2\log(2\kappa)-3}\\\
-\frac{\kappa^{2}-1}{\kappa^{2}+1}p_{1}p_{2}p_{3}+cDep_{1}^{2}p_{3}\alpha\end{bmatrix},$
(95)

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

$\alpha=-\frac{4}{\kappa^{2}}+\frac{p_{2}}{4\log(2\kappa)-6}\frac{2.5p_{1}De}{1+2.5p_{1}De/(3\pi\tan^{-1}\frac{p_{1}}{|p_{3}|})},$
(96)

Equation (95) 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
$\mathcal{O}(cDe\epsilon_{\dot{p}})$ for all the components
($\epsilon_{\dot{p}}$ is given in equation (76)) in the $De\ll 1$ limit. For
$De\gg 1$, the errors are of $\mathcal{O}(cDe\epsilon_{\dot{p}})$ in the 1 and
3 component and
$\mathcal{O}(c\max[\log(\kappa)/\kappa^{2},\epsilon_{\dot{p}}])$ in the 2
component. Before analyzing the influence of polymers on particle orientation
as suggested by this equation we compare our theory at low $De$ with that of
Leal (1975).

###  4.6 Comparison of second order fluid result with Leal (1975)

Leal (1975) considered the motion of a fiber in a second order fluid and
found,

$\displaystyle\dot{p}_{2}^{\text{Leal}}=$ $\displaystyle-
p_{1}p_{2}(p_{2}+Vp_{1}(1-2p_{2}^{2}))$ (97)
$\displaystyle\dot{p}_{3}^{\text{Leal}}=$ $\displaystyle-
p_{1}p_{2}p_{3}(1-2Vp_{1}p_{3}),$ (98)

where, $V=-\frac{3\lambda\gamma}{16\log(\kappa)}M_{1}(1+2\epsilon_{1})$,
$\lambda=\Phi_{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}=\Phi_{2}/\Phi_{3}$,

$\displaystyle\Phi_{2}=-\lim\limits_{\gamma\rightarrow
0}\frac{\sigma_{11}-\sigma_{22}}{2\gamma^{2}},\hskip
14.45377pt\text{and}\hskip 14.45377pt\Phi_{3}=\lim\limits_{\gamma\rightarrow
0}\frac{\sigma_{11}-\sigma_{33}}{\gamma^{2}},$ (99)

$\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=\frac{-3De}{8\log(\kappa)}M_{1}(1+2\epsilon_{1})$.
$(1+2\epsilon_{1})$ and hence $V$ in Leal’s theory is proportional to the
second normal stress difference. Up to $\mathcal{O}(Dep_{2})$ Leal’s results
are therefore,

$\displaystyle\dot{p}_{2}^{\text{Leal}}=-p_{1}p_{2}^{2}+\frac{3De}{8\log(\kappa)}M_{1}(1+2\epsilon_{1})p_{1}^{2}p_{2}+\mathcal{O}(De^{2},p_{2}^{2}De)$
(100)
$\displaystyle\dot{p}_{3}^{\text{Leal}}=-p_{1}p_{2}p_{3}+\mathcal{O}(De^{2},p_{2}^{2}De).$
(101)

By taking only up to $\mathcal{O}(De)$ terms from equation (95), our results
up to this order for a prolate spheroidal particle in a second order fluid
are,

$\displaystyle\dot{p}_{2}=-p_{1}p_{2}^{2}-\frac{cDe}{2\log(2\kappa)-3}p_{1}^{2}p_{2}+\mathcal{O}(cDe^{2},cp_{2}^{2}De)$
(102)
$\displaystyle\dot{p}_{3}=-p_{1}p_{2}p_{3}+\mathcal{O}(cDe^{2},cp_{2}^{2}De).$
(103)

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., 1991), will
exhibit a different rotational motion at a finite but small $c\cdot De$ as
compared to $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 (1977) considered the motion of rigid particles in
second order fluid. While Brunn (1977) 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 (1977) 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 Analysis of particle orientation dynamics

Before a more detailed analytical and numerical treatment of equation (95), we
qualitatively describe the changes in orientation dynamics introduced by
viscoelasticity as given by this equation. We obtain three primary regions in
$c\cdot De-\kappa$ space where viscoelasticity leads to different final fiber
orientation behaviors. These regions, $R_{vort}=R_{vort}^{(1)}\cup
R_{vort}^{(2)}$, $R_{flow}=R_{flow}^{(1)}\cup R_{flow}^{(2)}$, and, $R_{flow-
vort}=R_{flow-vort}^{(1)}\cup R_{flow-vort}^{(2)}\cup R_{flow-vort}^{(3)}$,
separated by the curves $\lambda_{2,flow}^{-}=0$ and $c\cdot De=c\cdot
De|_{cut-off}^{Limit}$ are shown in figure 2 for $c=0.005$. In $R_{vort}$,
i.e. the region to the left of the curve $\lambda_{2,flow}^{-}=0$ in figure 2,
irrespective of the initial orientation, the particle is attracted towards a
stable limit cycle near the vorticity axis where the particle revolves in a
small periodic orbit around the vorticity axis. In $R_{flow}$, i.e. the region
to the right of the curve $c\cdot De=c\cdot De|_{cut-off}^{Limit}$ in figure
2, a particle aligns close to the flow direction for all initial orientations
as a stable fixed point at this location is the only stable attractor in the
orientation space. This flow alignment behavior is similar to the particle
trajectories observed by Iso et al. (1996a) at large elasticity (large $c$).
In $R_{flow-vort}$, defined as the region between $\lambda_{2,flow}^{-}=0$ and
$c\cdot De=c\cdot De|_{cut-off}^{Limit}$ in figure 2, depending upon the
initial orientation, the particle can either obtain a final orientation within
the flow-gradient plane (it may either obtain a stable orientation or rotate
within the plane) or rotate in a small periodic orbit around the vorticity
axis.

Figure 2: $\kappa-c\cdot De$ parameter space, for $c=0.005$, divided into
regions with different qualitative behavior of particle’s orientation
dynamics. $b^{2}$ and $\zeta$ are given in equations (107) and (115)

respectively. $\lambda_{2,\text{flow}}^{0-}$ and
$\lambda_{2,\text{flow}}^{0--}$ are in equation (117). Procedure for
numerically obtaining $c\cdot De|_{\text{cut-off}}^{\text{Limit}}$ is
described in section 5.4.1.

There are subdivisions of the regions $R_{vort}$, $R_{flow-vort}$ and
$R_{flow}$ based on the behavior of the particle’s orientation trajectory near
the flow-gradient plane and the vorticity axis. In $R_{vort}^{(1)}$ (figure
2), a particle starting near the flow-gradient plane spirals away from the
plane and towards the limit cycle near the vorticity axis. In contrast, in
$R_{vort}^{(2)}$ (figure 2), once a particle comes close to the flow
direction, it drifts along the flow-vorticity plane in a monotonic fashion.
$R_{vort}^{(1)}$ trajectories are reminiscent of the observations of Gauthier
et al. (1971) at low shear rates (low $De$) and also some of the particle
trajectories of Iso et al. (1996b, a) at low to medium elasticity (small to
medium $c$). The trajectories in $R_{vort}^{(2)}$, on the other hand, are
similar to that observed by Bartram et al. (1975) at larger shear rates or
$De$ (than the earlier low shear rate experiments with the same fluid reported
by the same laboratory in Gauthier et al. (1971)). Within the region
$R_{flow}$ shown in figure 2, trajectories in $R_{flow}^{(1)}$ and
$R_{flow}^{(2)}$ have different behavior near the vorticity axis. In the
former, to the left of the curve $\kappa=c\cdot De$ in figure 2, the vorticity
axis is an unstable spiral. Thus the particle leaves the region close to the
vorticity axis in a spiral motion. In the latter, the vorticity axis is an
unstable node, and the particle leaves the vicinity of the vorticity axis in a
monotonic fashion. The trajectories near the vorticity axis in also in
$R_{vort}$ and $R_{flow-vort}$ are similar to those in $R_{flow}^{(1)}$. As
mentioned above, in the region of multiple final particle orientations,
$R_{flow-vort}$, in addition to a stable limit cycle/ periodic orbit close to
the vorticity axis, a stable attractor exists in the flow-gradient plane.
$R_{flow-vort}$ can be further subdivided into three sub-regions. In $R_{flow-
vort}^{(1)}$ (figure 2), the particle spirals toward the flow-gradient plane
and then tumbles within the plane. In this case, the particle slows down close
to the flow direction before speeding up again as it departs from this
orientation. In $R_{flow-vort}^{(2)}$ and $R_{flow-vort}^{(3)}$ (figure 2),
there is a stable fixed point very close to the flow direction. Hence a
particle starting close to the flow-gradient plane ends up being flow-aligned.
In $R_{flow-vort}^{(2)}$, there is an unstable node within the flow-gradient
plane (further away from the flow direction than the stable fixed point). In
contrast, in $R_{flow-vort}^{(3)}$, the unstable node is replaced with a
saddle point with its stable manifold perpendicular to the flow-gradient
plane. Therefore, the basin of attraction of the stable fixed point (with the
stable limit cycle being the other stable attractor) is larger for $R_{flow-
vort}^{(3)}$ than for $R_{flow-vort}^{(2)}$.

The above is a qualitative description of all the different behaviors
exhibited by the dynamical system of equation (95). However, the predictions
for regions of large $c\cdot De$ (such as $R_{flow}^{(2)}$) are speculative.
This is because the first normal stress difference in the polymer stress is
proportional to $c\cdot De$ and our regular perturbation theory is based on
the polymer stress being smaller than the Newtonian stress. We nevertheless
include these regions here to provide a complete description of the dynamical
system and later numerical studies of the original governing equations
described in section 2 may be useful in determining the quantitative and
qualitative validation of the theory considered here.

In the rest of this section we will derive the boundaries that determine the
aforementioned divisions of $c\cdot De-\kappa$ space and provide a more
detailed analytical and numerical treatment of equation (95). $b^{2}=0$ and
$\kappa=2c\cdot De$ boundaries shown in figure 2 will be determined in
sections 5.1 and 5.2 respectively. $\lambda_{2,\text{flow}}^{0-}=0$,
$\lambda_{2,\text{flow}}^{0--}=0$ and $\zeta=0$ boundaries are derived in
section 5.3. $c\cdot De|_{\text{cut-off}}^{\text{Limit}}$, is numerically
obtained (the other boundaries of figure 2 are analytical) in section 5.4.1.
To numerically integrate the orientation trajectory we transform equation (95)
into $\theta-\phi$ coordinates, where,

$\displaystyle\mathbf{p}=\begin{bmatrix}p_{1}\\\ p_{2}\\\
p_{3}\end{bmatrix}=\begin{bmatrix}\sin(\theta)\sin(\phi)\\\
\sin(\theta)\cos(\phi)\\\ \cos(\theta)\end{bmatrix}\rightarrow$
$\displaystyle\begin{split}&\frac{d\theta}{dt}=-\frac{1}{\sin(\theta)}\frac{\partial
p_{3}}{\partial t},\\\
&\frac{d\phi}{dt}=\frac{1}{\sin^{2}(\theta)}\Big{(}\frac{\partial
p_{2}}{\partial t}\cos(\phi)-\frac{\partial p_{1}}{\partial
t}\sin(\phi)\Big{)}\end{split}.$ (104)

Integrating this $\theta-\phi$ system instead of the equivalent equation (95)
directly for $\mathbf{p}$ numerically preserves $||\mathbf{p}||_{2}=1$
constraint.

The boundaries in the $c\cdot De-\kappa$ space shown in figure 2 are obtained
for $c=0.005$. The primary $c$ dependence of the various boundaries is in the
form of $c\cdot De$ and similar partitioning of the $c\cdot De-\kappa$ space
is found with different values of $c$, albeit with small quantitative changes.
The boundaries associated with $b^{2}=0$ (equation (107)), and $\kappa=2cDe$
only depend on $c\cdot De$ for a given $\kappa$. We find (section 5.4.1) the
$c\cdot De|_{\text{cut-off}}^{\text{Limit}}$ boundary to be insensitive to
changes in $c$ at constant $c\cdot De$ for $c\lessapprox 0.1$. However,
increasing $c$ moves the $\lambda_{2,\text{flow}}^{0-}=0$,
$\lambda_{2,\text{flow}}^{0--}=0$ and $\zeta=0$ boundaries in the $c\cdot
De-\kappa$ kappa space to the right or larger $c\cdot De$ (not shown) even at
small $c$. The $\lambda_{2,\text{flow}}^{0--}=0$ curve moves more than
$\lambda_{2,\text{flow}}^{0-}=0$ upon increasing $c$. Thus, in the $c\cdot
De-\kappa$ space, upon increasing $c$, $R_{vort}^{(1)}$, $R_{flow}^{(1)}$ and
$R_{flow}^{(2)}$ are unchanged, $R_{vort}^{(2)}$ and $R_{flow-vort}^{(2)}$ are
enlarged, and, $R_{flow-vort}^{(1)}$ and $R_{flow-vort}^{(3)}$ are reduced in
size.

###  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
flow-gradient plane ($p_{3}=0$). Due to symmetry, and as suggested by
equations (95), 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}^{\text{PIST}}}^{(1)}$ from equation (79) and
${{G}_{3}^{\text{Elastic}}}^{(1)}$ from equation (89), are both zero. Thus,
the polymers do not change the log-rolling angular velocity of the particle at
$\mathcal{O}(c)$.

In the flow-gradient plane, close to the flow direction $p_{1}\approx 1$, the
equation governing the particle orientation is,

$\dot{p}_{2}\approx-
p_{2}^{2}-\frac{1}{\kappa^{2}}-cDe\frac{p_{2}}{2\log(2\kappa)-3}.$ (105)

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

$p_{2}=-b\tan({(t-t_{0})b})\frac{cDe}{4\log(2\kappa)-6},$ (106)

where,

$b^{2}={\frac{1}{\kappa^{2}}-\Big{(}\frac{cDe}{4\log(2\kappa)-6}\Big{)}^{2}}.$
(107)

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

$p_{2,\text{flow}}^{0-}=-\frac{cDe}{4\log(2\kappa)-6}+\sqrt{-b^{2}},\hskip
14.45377ptp_{2,\text{flow}}^{0--}=-\frac{cDe}{4\log(2\kappa)-6}-\sqrt{-b^{2}}.$
(108)

There are two additional fixed points at $-p_{2,\text{flow}}^{0-}$ and
$-p_{2,\text{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\cdot De$ is small.
$p_{2,\text{flow}}^{0--}$ is further downstream. Within the flow-gradient
plane, particle trajectories approach $\pm p_{2,\text{flow}}^{0-}$, while they
depart $\pm p_{2,\text{flow}}^{0--}$. Therefore, in the $b^{2}<0$ regime, a
particle placed in the flow-gradient plane approaches a steady state
orientation $p_{2}=\pm p_{2,\text{flow}}^{0-}$. We will later observe in
section 5.3 that in the orientation space near the flow-direction off the
flow-gradient plane, 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 flow-gradient plane in
$b^{2}<0$ regime

Above, we analyzed 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 equation (108). These expressions are only valid when
the fixed points are near the flow direction. However, as $c\cdot 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,\text{flow}}^{-}$ moves
closer to the flow direction, while $p_{2,\text{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 ${dp_{2}}/{dt}$ and
obtain the improved expressions,

$\displaystyle\begin{split}\tilde{p}_{2,\text{flow}}^{0-}&=-\sqrt{\frac{-f^{2}/\kappa^{2}+{c^{2}De^{2}}/2-fc\cdot
De\sqrt{-\tilde{b}^{2}}}{{c^{2}De^{2}+f^{2}}}},\\\
\tilde{p}_{2,\text{flow}}^{0--}&=-\sqrt{\frac{-f^{2}/\kappa^{2}+{c^{2}De^{2}}/2+fc\cdot
De\sqrt{-\tilde{b}^{2}}}{{c^{2}De^{2}+f^{2}}}},\end{split}$ (109)

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

Figure 3: Variation of the fixed points’ location on the flow-gradient plane
with $c\cdot De$ for various $\kappa$ in the $\tilde{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 flow-gradient plane. We now move on to consider three-dimensional orbits
by first observing the particle’s orientation trajectories close to the
vorticity axis in section 5.2 and then the trajectories near the flow-gradient
plane in section 5.3.

###  5.2 Effect of viscoelasticity near the vorticity direction, $p_{3}\approx
1$

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

$\mathbf{p}^{0}_{\text{vort}}=\begin{bmatrix}0&0&1\end{bmatrix}.$ (110)

Due to the constraint $||\mathbf{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 $\begin{bmatrix}0&0\end{bmatrix}$, the eigenvalues are,

$\lambda_{1,\text{vort}}=2cDe/\kappa^{2}-1/\kappa\sqrt{4c^{2}De^{2}/\kappa^{2}-1},\lambda_{2,\text{vort}}=2cDe/\kappa^{2}+1/\kappa\sqrt{4c^{2}De^{2}/\kappa^{2}-1},$
(111)

with the corresponding eigenvectors

$\mathbf{v}_{1,\text{vort}}=\begin{bmatrix}\lambda_{1,\text{vort}}&1\end{bmatrix},\mathbf{v}_{2,\text{vort}}=\begin{bmatrix}1&\lambda_{2,\text{vort}}\end{bmatrix}.$
(112)

The fixed point at the vorticity axis undergoes a Hopf bifurcation at $c\cdot
De=0$ and is unstable for all finite $c\cdot De$. The vorticity axis is a
center (imaginary eigenvalues) when $c\cdot De=0$. In the presence of
viscoelasticity, $c\cdot De\neq 0$, it becomes a hyperbolic fixed point
(eigenvalues with non-zero real part). Therefore, the linear stability
analysis provides qualitative insight into effect of viscoelasticity on the
behavior of fiber orientation governed by the full nonlinear equation (95)
near the vorticity axis when $c\cdot De\neq 0$. This follows from the Hartman-
Grobman theorem (see Guckenheimer & Holmes (2013)), which guarantees the
homeomorphism between the linearized and full non-linear 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 $\mathcal{O}(1/\kappa^{2})$ Stokeslet discussed in
section 4.1. When $0<cDe<\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 $cDe>\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 equation (95) for $\kappa=50$ at $c\cdot De=24$ (also for $c\cdot De=5$ and
6) vs. $c\cdot De=26$ in the left panel of figure 4 zoomed near the vorticity
axis. At $\kappa=50$, $c\cdot De=26$ is a point in $R_{flow}^{(2)}$, $c\cdot
De=24$ and 6 are points in $R_{flow}^{(1)}$ and $c\cdot 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 the right panel of figure 4 in more detail
later in section 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\cdot De$.
Hence, the Hopf bifurcation occurring at $c\cdot De=0$ at the vorticity axis
is supercritical (Strogatz, 2018). In $R_{flow-vort}^{(3)}$ ($c\cdot 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,\text{flow}}^{0-}$) near the flow
direction in the flow-gradient plane for $c\cdot De=6$, 24 and 25.

(a)  (b)

Figure 4: Various orientation behaviors near the vorticity axis: trajectories
of particle orientation starting very close to the vorticity axis at different
$c\cdot De$ in $R_{flow}^{(1)}$ ($cDe=6$, 24), $R_{flow}^{(2)}$ ($cDe=26$),
and $R_{flow-vort}^{(3)}$ ($cDe=5$) at $\kappa=50$. All trajectories start at
the same point. Left panel is same as right panel (showing complete particle
trajectory) but zoomed near the vorticity axis ($p_{3}=1$). The gray 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 (1971)) along the fiber in addition to the flow generated by
the force per unit length (flow from Batchelor (1970)). 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 flow-gradient plane it will approach the vorticity axis. A previous
slender body theory for a second order fluid by Férec et al. (2021), relying
on only the flow generated by the force per unit length, predicts the fiber to
align with the vorticity axis. When the axis of symmetry of an axisymmetric
slender fiber 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 fiber axis. However, according to equations (29) and (35)
of Férec et al. (2021) when a fiber is in the flow-vorticity plane, the force
per unit length is nonzero and proportional to the polymer relaxation time.
This implies an error in the expression for the tension force in equation (29)
of Férec et al. (2021).

Figure 3a of Wang et al. (2019) 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 2a and 2b of d’Avino et al. (2014) 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 8a of the theoretical investigation of
Wang et al. (2020) 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 (2002) in an Oldroyd-B fuid 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 flow-gradient plane and flow
direction

We described the effect of viscoelasticity on the particle’s motion near the
flow direction within the flow-gradient plane (FGP) in section 5.1. Here, we
consider the particle motion near (but not exactly on) the FGP. The analysis
of section 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 equation (105) and (106). The
orientation dynamics in the $p_{3}$ (vorticity) direction for $p_{1}\approx 1$
is governed by the simplified equation (from equation (95)),

$\dot{p}_{3}\approx-
p_{2}p_{3}+cDep_{3}\Big{(}-\frac{4}{\kappa^{2}}+\frac{p_{2}}{8\log(2\kappa)-12}\frac{15\pi^{2}De}{3\pi^{2}+5De}\Big{)}+\mathcal{O}(p_{3}^{2}).$
(113)

Its closed form solution is,

$p_{3}=C\exp(cDe\zeta(t-t_{0}))\cos({(t-t_{0})b})^{-1+cDe\gamma/(8\log(2\kappa)-12)}$
(114)

where $b^{2}$ is given in equation (107) and,

$\gamma=\frac{15\pi^{2}De}{3\pi^{2}+5De},\hskip
7.22743pt\zeta=\frac{1}{4\log(2\kappa)-6}-cDe\frac{\gamma}{2(4\log(2\kappa)-6)^{2}}-\frac{4}{\kappa^{2}}.$
(115)

When $b$ is real-valued, i.e. $b^{2}>0$ (equation (107)), there are no fixed
points on the FGP. Thus, the particle undergoes periodic motion in the FGP, as
discussed in section 5.1. For a particle perturbed slightly away from the FGP,
from equation (114) we note that the sign of $\zeta$ determines whether the
particle will be attracted to or repelled from the FGP ($p_{3}=0$). $\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\pi/b$ than in Newtonian case). In the $b^{2}>0$
regime, the phase portrait of the dynamical system of equation (95) (through
the simplified equations (105) and (113)) projected in $p_{2}-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}-p_{3}$ plane in the case of a Newtonian fluid ($c\cdot
De=0$) is similar to that in the region $R_{vort}^{(1)}$ (left panel of figure
5), 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\cdot De$, i.e. in the region
$R_{vort}^{(1)}$, represented by the left panel of figure 5, 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 the right panel of figure 5 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)}$.

(a)  (b)

Figure 5: Phase portraits of equations (105) and (113) in the $b^{2}>0$ regime
in the gradient ($p_{2}$)- vorticity ($p_{3}$) plane with $\zeta>0$ (left) and
$\zeta<0$ (right). When $b^{2}>0$ and $\zeta>0$, i.e., in the left panel
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 the right panel 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$, left panel) and $R_{flow-vort}^{(1)}$
($\zeta<0$, right panel) 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 spiraling rate in $R_{vort}^{(1)}$ increases with $c\cdot De$ and $\kappa$
(not shown).

(a)  Figure 6: In $R_{vort}^{(1)}$ (shown here for $c=0.005$,
$\kappa=50,c\cdot De=0.01$), trajectories starting near the flow-gradient
plane (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.  (a)  Figure 7: In
$R_{flow-vort}^{(1)}$ (shown here for $c=0.005$, $\kappa=10,c\cdot De=0.48$),
trajectories of particle orientation starting near the flow-gradient plane
(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\cdot De$ or $\kappa$, the time period,
$2\pi/b$, increases. The effect of viscoelasticity (increasing $c\cdot De$) to
increase the time period is consistent with the experimental observations of
Gauthier et al. (1971) and Bartram et al. (1975). As mentioned earlier in
section 5.1, when $b^{2}=0$, due to an infinite period bifurcation (Strogatz,
2018) two fixed points emerge on the FGP, and we discuss the trajectories off
the FGP in the $b^{2}<0$ regime next.

####  ${b^{2}<0}$: Monotonic behavior near flow-gradient plane

Equation (105) has two fixed points at $p_{2,\text{flow}}^{0-}$ and
$p_{2,\text{flow}}^{0--}$ from equation (108) (or their more accurate values
in equation (109)). For the dynamical system in the $p_{2}-p_{3}$ plane near
the flow direction, defined by the system of equations (105) and (113), these
fixed points are $[p_{2,\text{flow}}^{0-},0]$ and
$[p_{2,\text{flow}}^{0--},0]$. The eigenvectors at each of the fixed points
are the same, i.e.

$\mathbf{v}_{1,\text{flow}}=\begin{bmatrix}1&0\end{bmatrix},\mathbf{v}_{2,\text{flow}}=\begin{bmatrix}0&1\end{bmatrix}.$
(116)

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

$\displaystyle\begin{split}&\lambda_{1,\text{flow}}^{0-}=-2\sqrt{-b^{2}},\hskip
14.45377pt\lambda_{2,\text{flow}}^{0-}=\zeta
cDe-(\zeta+4/\kappa^{2})(4\log(2\kappa)-6)\sqrt{-b^{2}},\\\
&\lambda_{1,\text{flow}}^{0--}=2\sqrt{-b^{2}},\hskip
14.45377pt\lambda_{2,\text{flow}}^{0--}=\zeta
cDe+(\zeta+4/\kappa^{2})(4\log(2\kappa)-6)\sqrt{-b^{2}}.\end{split}$ (117)

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 non-linear behavior of particle orientation. Both eigenvalues of both
the fixed points are real in the $b^{2}<0$ regime. The first eigenvalues
($\lambda_{1,\text{flow}}^{0-}$ and $\lambda_{1,\text{flow}}^{0--}$) of both
the fixed points do not change sign in $b^{2}<0$ regime i.e.
$\lambda_{1,\text{flow}}^{0-}<0$ and $\lambda_{1,\text{flow}}^{0--}>0$. But,
the second eigenvalues ($\lambda_{1,\text{flow}}^{0-}$ and
$\lambda_{1,\text{flow}}^{0--}$), that are initially positive in the $b^{2}<0$
regime, become negative at different parameter ($c$, $De$ and $\kappa$)
values. $\lambda_{2,\text{flow}}^{0-}$ becomes negative first.
$\lambda_{2,\text{flow}}^{0-}=0$ is the bifurcation boundary where
$[p_{2,\text{flow}}^{0-},0]$ changes from a saddle point (stable manifold
along gradient direction and unstable along vorticity) to a stable fixed
point. $\lambda_{2,\text{flow}}^{0--}=0$ is the bifurcation boundary where
$[p_{2,\text{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
labeled in figure 2 as (1) $R_{vort}^{(2)}$ ($\lambda_{1,\text{flow}}^{0-}>0$
and $\lambda_{1,\text{flow}}^{0--}>0$), (2) $R_{flow-vort}^{(2)}$
($\lambda_{1,\text{flow}}^{0-}<0$ and $\lambda_{1,\text{flow}}^{0--}>0$), and
(3) $R_{flow-vort}^{(3)}+R_{flow}$ ($\lambda_{1,\text{flow}}^{0-}<0$ and
$\lambda_{1,\text{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,\text{flow}}^{0-},0]$ and
leaving $[p_{2,\text{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 flow-gradient plane (FGP) will approach the flow direction
(specifically the fixed point $[p_{2,\text{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. (1975) discussed in section 1. For a given
$\kappa$ and $c$, at larger $c\cdot De$ than $R_{vort}^{(2)}$, within the
regions $R_{flow-vort}^{(2)}$ and $R_{flow-vort}^{(3)}$, the particle starting
near the flow-gradient plane (and for all starting orientations in $R_{flow}$)
achieves a stable orientation near the flow direction, similar to the
experiments of Iso et al. (1996a).

(a)

(b)  (c)

Figure 8: Phase portraits of the system of equations (105) and (113) in the
$b^{2}<0$ regime. In this regime, two fixed points exist on the flow-gradient
plane ($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 the top panel ($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 the bottom left
($R_{flow-vort}^{(2)}$) and right ($R_{flow-vort}^{(3)}$) panels indicate that
particles with starting orientation near the flow gradient plane may align
near the flow direction. The presence of an unstable point in $R_{flow-
vort}^{(2)}$ (bottom left) instead of a saddle point with its stable manifold
perpendicular to the flow gradient plane in $R_{flow-vort}^{(2)}$ (bottom
right) 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,\text{flow}}^{0-}$ and $\tilde{p}_{2,\text{flow}}^{0--}$ in the
orientation space are
$\begin{bmatrix}\pm\sqrt{1-(\tilde{p}_{2,\text{flow}}^{0-})^{2}}&\pm\tilde{p}_{2,\text{flow}}^{0-}&0\end{bmatrix}$
and
$\begin{bmatrix}\pm\sqrt{1-(\tilde{p}_{2,\text{flow}}^{0--})^{2}}&\pm\tilde{p}_{2,\text{flow}}^{0--}&0\end{bmatrix}$.
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 plots of figures 9, 10, and, 11
respectively. We mark the analytical locations of these fixed points with two
different colored markers on the flow-gradient plane ($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.

(a)  Figure 9: In $R_{vort}^{(2)}$ (shown here for $c=0.005$,
$\kappa=50,c\cdot 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 flow-
gradient plane, here they leave the flow-gradient plane monotonically. We show
the region close to the flow-vorticity plane as this is where the different
attractors lie. The gray surface is the unit sphere, i.e., the orientation
space.

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 flow-gradient
plane (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. (1975). We can only make
qualitative comparisons with Bartram et al. (1975) 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 flow-
gradient plane 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 further away on the flow-gradient plane 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 flow-gradient plane approach a stable
orientation near the flow direction (stable fixed point at
$\tilde{p}_{2,\text{flow}}^{0-}$).

(a)  Figure 10: In $R_{flow-vort}^{(2)}$ (shown here for $c=0.005$,
$\kappa=100$ and $c\cdot De=0.9$) a particle’s orientation trajectories that
start close to the flow-gradient plane (solid lines) approach the flow
direction either on the same or the opposite side of the gradient-vorticity
plane. Trajectories further away from the flow-gradient plane (dashed lines)
approach the stable limit cycle near the vorticity axis. Due to a saddle and
stable node close to each other on the flow-gradient plane and a stable limit
cycle near the vorticity axis, another saddle-node emerges on the faster
eigen-direction of the stable node. The gray surface is the unit sphere, i.e.,
the orientation space.  (a)  Figure 11: In $R_{flow-vort}^{(3)}$ (shown here
for $c=0.005$, $\kappa=100$ and $c\cdot De=1.2$) the behavior 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 flow-gradient
plane further 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 flow-gradient 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.

###  5.4 Global particle orientation dynamics

In section 5.1 we analyzed the equations for the rotational motion of a
particle at the vorticity axis and within the flow-gradient plane (FGP). We
found that viscoelasticity does not lead to a nonzero $\dot{\mathbf{p}}$ at
the vorticity axis and it does not change the log-rolling rotational velocity.
In the flow-gradient plane, we observed that a small amount of viscoelasticity
(small $c\cdot De$) leads to a larger orbit time period. Further increasing
$c\cdot 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 section 5.2, by linearizing the governing equation at the
vorticity axis, we showed that the vorticity axis is an unstable spiral for
$c\cdot De<\kappa$. For $c\cdot 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