%20-%20https://sketch.com%20--%3e%3ctitle%3e3ECBBF50-FD08-412A-943D-2CFC5E025A16%3c/title%3e%3cdesc%3eCreated%20with%20sketchtool.%3c/desc%3e%3cg%20id='Stylesheet'%20stroke='none'%20stroke-width='1'%20fill='none'%20fill-rule='evenodd'%3e%3cg%20transform='translate(-926.000000,%20-668.000000)'%20id='BUTTONS'%3e%3cg%20transform='translate(100.000000,%20597.000000)'%3e%3cg%20id='BTN-/-RIGHTS'%20transform='translate(816.000000,%2065.000000)'%3e%3cg%20id='ICON-/-Rights'%20transform='translate(10.000000,%206.000000)'%3e%3cpath%20d='M13.5752841,0%20L1.34943182,0%20C0.59375,0%200,0.586666667%200,1.33333333%20L0,22.6666667%20C0,23.4133333%200.620738636,24%201.34943182,24%20L17.6505682,24%20C18.40625,24%2019,23.3866667%2019,22.6666667%20L19,5.36%20L13.5752841,0%20Z%20M13.5752841,1.89333333%20L17.0838068,5.36%20L13.5752841,5.36%20L13.5752841,1.89333333%20Z%20M1.34943182,22.6666667%20L1.34943182,1.33333333%20L12.1988636,1.33333333%20L12.1988636,5.36%20C12.1988636,6.10666667%2012.7926136,6.69333333%2013.5482955,6.69333333%20L17.6235795,6.69333333%20L17.6505682,22.6666667%20L1.34943182,22.6666667%20Z'%20id='Shape'%20fill='%23595959'%20fill-rule='nonzero'%3e%3c/path%3e%3ccircle%20id='Oval'%20stroke='%23595959'%20stroke-width='1.5'%20fill='%23595959'%20cx='9.5'%20cy='13.5355339'%20r='5.5'%3e%3c/circle%3e%3cpolyline%20id='Path-2'%20stroke='%23FFFFFF'%20stroke-width='1.5'%20transform='translate(9.500000,%2013.035534)%20rotate(45.000000)%20translate(-9.500000,%20-13.035534)%20'%20points='8%2016.0355339%2011%2016.0355339%2011%2010.0355339'%3e%3c/polyline%3e%3c/g%3e%3c/g%3e%3c/g%3e%3c/g%3e%3c/g%3e%3c/svg%3e)
1. Introduction
Characterizing the rheological behaviour of a dilute Stokesian suspension of hydrodynamically interacting particles requires the solution of the two-body problem at the microscale. Einstein (Reference Einstein1906) first obtained the effective viscosity (
$\mu _e$) of a dilute suspension of non-interacting rigid spherical particles (the one-body problem) in terms of the viscosity of the suspending fluid 
$(\mu )$ as 
$\mu _e = \mu (1+\frac {5}{2}\phi )$, 
$\phi$ being the volume fraction of the suspended spheres; the factor 
$5/2$ is referred to as the Einstein coefficient. Much later, Batchelor & Green (Reference Batchelor and Green1972a,Reference Batchelor and Greenb) investigated the effect of pair interactions between spherical particles, in an ambient linear flow, in an attempt to calculate the 
$O(\phi ^{2})$ correction to the effective viscosity. In obtaining this correction, they examined the pair-sphere pathlines in simple shear flow, and showed that these pathlines had the same character as the streamlines around a single sphere in a simple shear flow (Cox, Zia & Mason Reference Cox, Zia and Mason1968). In both instances, the trajectories may be classified into two groups. The first group consists of fore–aft symmetric open trajectories that extend to upstream and downstream infinity. The second group consists of closed orbits around the test sphere. The two groups are demarcated by a separatrix surface consisting of open trajectories that asymptote to the flow axis infinitely far away in both the upstream and downstream directions. The surface itself is axisymmetric, the axis of symmetry being the gradient direction of the ambient simple shear. Both the fore–aft symmetry of open trajectories, and the existence of closed ones, arise from the reversibility of the Stokes equations. Closed streamlines in the single-sphere problem are known to profoundly affect heat and mass transfer, leading to diffusion-limited transport at large Péclet numbers (Acrivos Reference Acrivos1971; Subramanian & Koch Reference Subramanian and Koch2006a,Reference Subramanian and Kochc; Krishnamurthy & Subramanian Reference Krishnamurthy and Subramanian2018). Likewise, for the pair-sphere problem, Batchelor & Green (Reference Batchelor and Green1972a) showed that the pair-distribution function is rendered indeterminate, in the region of closed pathlines, in the pure hydrodynamic limit. This in turn leads to an indeterminate 
$O(\phi ^{2})$ coefficient for the suspension stress in all linear flows with regions of closed pair pathlines, simple shear flow being a special case (Kao, Cox & Mason Reference Kao, Cox and Mason1977).
In this paper, we analyse the fluid pathlines around a torque-free neutrally buoyant spheroid of an arbitrary aspect ratio (regarded as a canonical anisotropic particle) in simple shear flow. Keeping in mind the scenario for suspensions of spherical particles detailed above, we expect our findings to, on one hand, shed some light on the transport of heat or mass from a single spheroid immersed in a shearing flow. On the other hand, our findings will also be relevant to the nature of pair-spheroid interactions, and thence, to the pair-level microstructure of Stokesian suspensions of interacting anisotropic particles. In any event, analysing the fluid motion around a single spheroid is a natural first step towards an understanding of the more complicated two-body problem. The results reported here suggest that the diffusion-limited scalar transport at 
$O(\phi )$, and the indeterminate rheology at 
$O(\phi ^{2})$, encountered for the case of spherical particles above, are both likely to be resolved, to a significant extent, by the effects of shape anisotropy.
We characterize the topology of the fluid pathlines around a single prolate or oblate spheroid in an ambient simple shear flow. The exact disturbance velocity field required for this purpose is available in closed form from earlier efforts (Dabade, Marath & Subramanian Reference Dabade, Marath and Subramanian2015, Reference Dabade, Marath and Subramanian2016; Marath & Subramanian Reference Marath and Subramanian2018) that have used a spheroidal harmonics formalism for this purpose. As is well known, a neutrally buoyant torque-free spheroid in simple shear rotates along Jeffery orbits, and the disturbance velocity field is therefore time dependent in all cases except when the spheroid axis is aligned with the ambient vorticity (the log-rolling or spinning mode). Our results show that the departure from sphericity has a profound effect on the nature of the fluid pathlines around a spheroid in a precessional (non-spinning) orbit. The obvious changes happen for the open pathlines which are still open, but unlike the case of a sphere, are no longer fore–aft symmetric. Such pathlines undergo a lateral displacement in both the gradient and vorticity directions as they proceed from upstream to downstream infinity. The unexpected aspect concerns the originally closed pathlines. A fraction of these pathlines open up, but unlike the pathlines above, the resulting open pathlines do not head directly downstream. Instead, these pathlines loop around the rotating spheroid with the number of loops being extremely (indeed, infinitely) sensitive to the upstream coordinates of the particular pathline. The sensitivity with regard to the upstream coordinates has the signatures of chaotic scattering (Aref & Pomphrey Reference Aref and Pomphrey1980, Reference Aref and Pomphrey1982; Aref Reference Aref1983; Bleher, Grebogi & Ott Reference Bleher, Grebogi and Ott1990; Aref & Stremler Reference Aref and Stremler1999). The presence of a chaotic saddle suggested by the aforementioned scattering pattern implies that both closed (periodic) and bounded pathlines continue to exist for a spheroid in a generic precessional orbit, but only constitute a set of vanishingly small measure, and one that is therefore numerically inaccessible.
The paper is organized as follows. In § 2 we discuss the problem formulation. Calculation of the fluid pathlines requires an expression for the disturbance velocity field due to a torque-free spheroid in simple shear flow, and as mentioned above, this is obtained in terms of the appropriate vector spheroidal harmonics (Kushch Reference Kushch1997, Reference Kushch1998, Reference Kushch2013), after resolving the ambient simple shear into simpler canonical components in body-fixed coordinates (Dabade et al. Reference Dabade, Marath and Subramanian2015, Reference Dabade, Marath and Subramanian2016; Marath & Subramanian Reference Marath and Subramanian2018). In § 3 we first analyse the fluid pathlines (which are the same as streamlines) for the simpler steady case of a spinning spheroid. In this case, the topology of the streamline pattern is shown to be identical to that around a torque-free sphere and a circular cylinder (Torza et al. Reference Torza, Henry, Cox and Mason1971; Kao et al. Reference Kao, Cox and Mason1977; Powell Reference Powell1983), with fore–aft symmetric open and closed streamlines being demarcated by a (non-axisymmetric) separatrix surface. For a prolate spheroid, the streamline pattern (including the separatrix in particular) transitions from that for a sphere to a cylinder, with increasing spheroid aspect ratio. In §§ 4.1 and 4.2, we investigate the topology of the fluid pathlines for the general time dependent case when the spheroid is not in its spinning orbit. Section 4.1 concerns regular open pathlines that do not loop around the spheroid, while heading from upstream to downstream infinity. Unlike both a sphere and a spinning spheroid, these pathlines are fore–aft asymmetric. This asymmetry is best characterized in terms of the net gradient and vorticity displacements suffered by the pathline (Da Cunha & Hinch Reference Da Cunha and Hinch1996; Subramanian & Brady Reference Subramanian and Brady2006). The numerically determined displacements are then compared to the far-field analytical predictions, as a function of the upstream pathline coordinates, for the spheroid rotating over a range of Jeffery orbits. In § 4.2 we analyse the second category of pathlines which includes the singular open pathlines that, although open, loop around the spheroid before heading off to downstream infinity. For this case, we first show that the number of loops around the spheroid, and therefore, the residence time of the pathline in a sufficiently large neighbourhood of the spheroid, is sensitively dependent on the initial pathline coordinates (chosen to lie on the negative flow axis). Plots of the residence time reveal an extremely fine-scaled fractal-like structure, suggesting the existence of a chaotic saddle in the vicinity of the spheroid. We examine the nature of this sensitive dependence as a function of the spheroid aspect ratio, and more importantly, as a function of the orbit constant (C) for a given aspect ratio; the latter reveals the singular nature of the approach to the integrable limit (
$C = 0$) of a log-rolling or a spinning spheroid. In § 5, we analyse the fluid pathlines from the chaotic scattering perspective that has been gainfully used in other scenarios (Aref Reference Aref1983; Bleher et al. Reference Bleher, Grebogi and Ott1990; Jánosi et al. Reference Jánosi, Tél, Wolf and Gallas1997; Faisst & Eckhardt Reference Faisst and Eckhardt2004), by focussing on the nature of the pathlines as a function of their upstream coordinates and the Jeffery phase angle. Herein, we examine the transition from regular to chaotic scattering with varying spheroid aspect ratio, and again, as a function of the Jeffery orbit constant. There exists a critical upstream offset corresponding to the onset of chaotic scattering, and for offsets less than this value, intervals of regular and chaotic scattering appear interlaced down to the smallest scales; the interlaced pattern is strongly dependent on aspect ratio. This critical offset also varies with the Jeffery phase angle, and one may therefore characterize the scattering in terms of a residence-time surface plotted as a function of the upstream offset and the Jeffery phase angle, with a critical curve that separates regions of regular and irregular dependence. The extent of chaotic scattering is nevertheless shown to decrease for the largest aspect ratios, with the relevant signatures being eventually undetectable (numerically). We present the fractal nature of the stable and the unstable manifold of the chaotic saddle along with the uncertainty dimension plots for spheroids of varying aspect ratios, which shows the necessary power law scaling for the underlying saddle. We also examine the probability density of residence times, the emphasis being on the functional form of the tail, which are determined by the large residence times associated with the singular pathlines; rather unexpectedly, the tails appear to be algebraic, suggesting that the underlying (chaotic) invariant set has a non-hyperbolic character. Section 6 presents a summary of our findings, a discussion of the implications of these findings in the context of transport in multiphase systems, and the expected nature of the fluid pathlines in a more general setting – the one-parameter family of hyperbolic linear flows.
2. Single spheroid in a simple shear flow
2.1. Jeffery orbits
The motion of spheroidal particles in linear flows in the absence of inertia is known since the work of Jeffery (Reference Jeffery1922). In simple shear flow, a single torque-free spheroid rotates in any one of a single parameter family of closed trajectories now known as Jeffery orbits. The motion along a Jeffery orbit may be characterized in terms of the polar 
$(\theta _j)$ and azimuthal 
$(\phi _j)$ angles of the spheroid axis (see figure 1), and these are given by the following functions of time:
\begin{equation} \phi_j=\cot^{{-}1}\left[\kappa\tan\left(\frac{\kappa\dot{\gamma}t}{1+\kappa^{2}}\right)\right],\quad \theta_j=\tan^{{-}1}\left[\frac{C\kappa}{\sqrt{\kappa^{2} \sin^{2}\phi_j+\cos^{2}\phi_j}}\right]. \end{equation}
Here, 
$\dot {\gamma }$ is the shear rate, 
$\kappa$ is the spheroid aspect ratio and 
$C$ is the orbit constant, ranging from 
$0$ to 
$\infty$, that parameterizes the Jeffery orbits; 
$\kappa >1({<}1)$ for prolate (oblate) spheroids. The generic orbit is a spherical ellipse, with its major axis along the flow (gradient) direction for a prolate (an oblate) spheroid. The limiting cases 
$C = 0$ and 
$C = \infty$ correspond to a log-rolling motion of the prolate spheroid (spinning for the oblate one) about the vorticity axis, and a tumbling motion in the flow–gradient plane, respectively. The period of rotation is independent of 
$C$, being equal to 
$T_j = 2{\rm \pi} /\dot {\gamma }(\kappa + \kappa ^{-1})$. Figure 2 shows the Jeffery orbits for prolate and oblate spheroids of different aspect ratios. The eccentricity of the orbits can be seen to increase with increasing (decreasing) aspect ratio of the prolate (oblate) spheroid. For slender fibres (
$\kappa = \infty$) and flat disks (
$\kappa = 0$), the Jeffery orbits have a meridional character with the end points of the meridians, the poles, corresponding to the intersections of the flow and gradient axes with the unit sphere, respectively.
A prolate and an oblate spheroid in simple shear flow. The primed (
$X',Y',Z'$) and unprimed (
$X,Y,Z$) coordinates denote the body and the space-fixed axes, respectively. Here, 
$\theta _j$ and 
$\phi _j$ are polar and azimuthal Jeffery angles, respectively.
Jeffery orbits for prolate (top) (
$\kappa =1.25, 2, 4$) and oblate spheroids (
$\kappa =0.8, 0.5, 0.25$).
2.2. The spheroid velocity field in simple shear flow
The determination of the velocity field due to an inertialess torque-free spheroid in an ambient simple shear flow involves solving the (dimensionless) Stokes equations in spheroidal coordinates (Kushch Reference Kushch1997, Reference Kushch1998, Reference Kushch2013)
\begin{equation} \nabla^{2}\boldsymbol{u} = \boldsymbol{\nabla} p, \end{equation}with the following boundary conditions in place:
\begin{equation} \boldsymbol{u}(\xi=\xi_0,\eta,\phi)=\boldsymbol{\omega}\times\boldsymbol{x}\quad\textrm{and}\quad \boldsymbol{u}(\xi\to\infty,\eta,\phi)=\boldsymbol{\varGamma}\boldsymbol{\cdot}\boldsymbol{x}. \end{equation}
Here, u and 
$p$ are the velocity and pressure fields, 
$\boldsymbol {x} \equiv (x,y,z)$ is the position vector and 
$\boldsymbol {\varGamma }$ is the (transpose of the) velocity gradient tensor; 
$\boldsymbol {\varGamma } = \boldsymbol {1}_{x}\boldsymbol {1}_{y}$ in the space-fixed coordinate system (
$X,Y,Z$) shown in figure 1, with 
$\boldsymbol {1}_{x}$, 
$\boldsymbol {1}_{y}$ and 
$\boldsymbol {1}_{z}$ denoting the unit vectors along the flow, gradient and vorticity directions, respectively. The spheroid angular velocity, 
$\boldsymbol {\omega }$, in (2.3a,b) is determined from a torque-free constraint. The component of 
$\boldsymbol {\omega }$ normal to the spheroid axis, is obtained from the Jeffery relations (2.1a,b), while that along the axis is just half the projected ambient vorticity. In (2.2) and (2.3a,b), (
$\xi ,\eta ,\phi$) are prolate spheroidal coordinates which, for a prolate spheroid, are related to the Cartesian ones as: 
$x'=\textrm {d}\bar {\xi }\bar {\eta }\cos \phi$, 
$y'=\textrm {d}\bar {\xi }\bar {\eta }\sin \phi$, 
$z'=\textrm {d}\xi \eta$, where 
$\bar {\xi }=\sqrt {\xi ^{2}-1}$, and 
$\bar {\eta }=\sqrt {1-\eta ^{2}}$ (for an oblate spheroid, the corresponding relations are obtained by setting 
$d\rightarrow -\mathrm {i}d$ and 
$\xi \rightarrow \mathrm {i}\bar {\xi }$). Here, the constant-
$\xi$ surfaces represent a family of confocal prolate (oblate) spheroids with the interfocus distance (diameter of the focal circle) being 
$2d$; 
$\xi =\xi _0$ denotes the surface of the spheroidal particle. The constant-
$\eta$ surfaces denote a family of confocal two-sheeted (single-sheeted) hyperboloids, and the constant-
$\phi$ surfaces are planes passing through the spheroidal axis of symmetry.
The total velocity field 
$\boldsymbol {u}$ may be written as 
$\boldsymbol { u} = \boldsymbol {\varGamma }\boldsymbol {\cdot }\boldsymbol {x} + \boldsymbol {u'}$ where 
$\boldsymbol {u'}$ is the disturbance velocity field that is a function of the instantaneous spheroid orientation and vanishes in the far field. The solution for 
$\boldsymbol {u'}$ is best accomplished in a body-fixed coordinate system. Defining (
$X',Y',Z'$) as the body-fixed coordinate system (see figure 1), the polar angle 
$\theta _j$ denoting the angle between the 
$Z$ and 
$Z'$ axes, and the azimuthal angle 
$\phi _j$ denoting the dihedral angle between the 
$X$–
$Z$ and 
$X'$–
$Z'$ planes, with the 
$Y'$ axis constrained to lie in the flow gradient (
$X$–
$Y$) plane, the rate of strain tensor, 
$\boldsymbol{\mathsf{E}} = (\boldsymbol {\varGamma } + \boldsymbol {\varGamma }^{T})/2$, is given by
\begin{equation} {\boldsymbol{\mathsf{E}}}' = \left[ \begin{array}{@{}ccc@{}} \dfrac{\cos^{2}\theta_j\sin 2\phi_j}{2} & \dfrac{\cos\theta_j\cos 2\phi_j}{2} & \dfrac{\sin 2\theta_j\sin 2\phi_j}{4} \\ \dfrac{\cos\theta_j\cos 2\phi_j}{2} & -\dfrac{\sin 2\phi_j}{2} & \dfrac{\sin\theta_j\cos 2\phi_j}{2} \\ \dfrac{\sin 2\theta_j\sin 2\phi_j}{4} & \dfrac{\sin\theta_j\cos 2\phi_j}{2} & \dfrac{\sin^{2}\theta_j\sin 2\phi_j}{2} \end{array} \right] \end{equation}in body-fixed coordinates. As shown by Subramanian & Koch (Reference Subramanian and Koch2006b), in the context of a nearly spherical particle, the rate of strain tensor above can be resolved into five elementary components, each corresponding to a canonical linear flow. The corresponding disturbance velocity fields may be expressed in terms of the appropriate decaying vector spheroidal biharmonics. Following Dabade et al. (Reference Dabade, Marath and Subramanian2015, Reference Dabade, Marath and Subramanian2016), the disturbance velocity field for an arbitrary aspect ratio spheroid may similarly be written as
\begin{equation} \boldsymbol{u'}(\boldsymbol{x'};\theta_j,\phi_j) = \sum_{i=1}^{5}\boldsymbol{u}_{is}(\boldsymbol{x'};\theta_j,\phi_j), \end{equation}
where the 
$\boldsymbol {u}_{is}(\boldsymbol {x'};\theta _j,\phi _j)$ values refer to the component disturbance fields, and are defined in table 1.
Canonical Stokesian velocity fields, associated with an arbitrarily oriented spheroid, in a body-fixed coordinate system.
Here, 
$\boldsymbol {u}_{1s}$ corresponds to an axisymmetric extension with the axis of symmetry coincident with the spheroid axis; 
$\boldsymbol {u}_{2s}$ and 
$\boldsymbol {u}_{3s}$ correspond to a pair of planar extensions in the plane orthogonal to the spheroid axis, with their principal axes rotated by an angle of 
${\rm \pi} /4$ relative to each other, and 
$\boldsymbol {u}_{4s}$ and 
$\boldsymbol {u}_{5s}$ correspond to planar extensions in a pair of orthogonal planes containing the spheroid axis. The decaying spheroidal biharmonics in the expressions for the disturbance velocity fields in table 1 are of the general form 
$\boldsymbol {S}^{(3)}_{ts}$, and are defined as
\begin{align} \boldsymbol{S}^{(3)}_{ts}&=\boldsymbol{e}_1\left\{-(x'-\mathrm{i}y')D_2F_{t-1}^{s-1}-\bar{\xi}_0^{2}\, \textrm{d}D_1F_{t}^{s}+(t+s-1)(t+s)\beta_{-(t+1)}F_{t-1}^{s-1}\right\} \nonumber\\ &\quad + \boldsymbol{e}_2\left\{(x'+\mathrm{i}y')D_1F_{t-1}^{s+1}- \bar{\xi}_0^{2}\,\textrm{d}D_2F_{t}^{s}-(t-s-1)(t-s)\beta_{-(t+1)}F_{t-1}^{s+1}\right\} \nonumber\\ &\quad + \boldsymbol{1}_{z'}\left\{z'D_3F_{t-1}^{s}-\xi_0^{2}\,\textrm{d}D_3F_t^{s}-C_{-(t+1),s}F_{t-1}^{s}\right\}, \end{align}where,
\begin{gather} F_t^{s}(\boldsymbol{r},d)=Q_t^{s}(\xi)Y_t^{s}(\eta,\phi), \quad Y_t^{s}(\eta,\phi)=\frac{(t-s)!}{(t+s)!}P_t^{s}(\eta)\exp(\mathrm{i}\,s\phi) \end{gather}
\begin{gather}\beta_t=\frac{t+3}{(t+1)(2t+3)}, \quad C_{t,s}=(t+s+1)(t-s+1)\beta_t, \quad |s| \leqslant t \end{gather}
\begin{gather}\boldsymbol{e}_1=\frac{\boldsymbol{1}_{x'}+\mathrm{i}\,\boldsymbol{1}_{y'}}{2}, \quad \boldsymbol{e}_2=\frac{\boldsymbol{1}_{x'}-\mathrm{i}\,\boldsymbol{1}_{y'}}{2} \end{gather}
\begin{gather}D_1=(\partial_{x'}-\textrm{i}\partial_{y'}), \quad D_2=(\partial_{x'}+\textrm{i}\partial_{y'}), \quad D_3=\partial_{z'}. \end{gather}
Here, 
$Y_t^{s}$ are the scalar surface harmonics, with 
$P_t^{s}$ and 
$Q_t^{s}$ being the associated Legendre functions of the first and second kind, respectively; the index 
$t$ in 
$\boldsymbol {S}^{(3)}_{ts}$ denotes the rapidity of the decay for large 
$\xi$ (or 
$r$), with 
$\lim _{r\to \infty }\boldsymbol {S}^{(3)}_{ts}\propto r^{-t}$. Thus, all of the 
$\boldsymbol {S}^{(3)}_{2s}$ values appearing in the disturbance velocity fields (table 1) decay as 
$1/r^{2}$, as is appropriate since the force-free freely rotating spheroid appears as a fluctuating force dipole in the far field.
The complete dynamical system governing the evolution of the fluid pathlines in terms of the position vector (
$\boldsymbol {x}$) and Jeffery angles (
$\theta _j,\phi _j$) may now be written as
\begin{gather} \frac{\textrm{d}\boldsymbol{x}}{\textrm{d}t}=\boldsymbol{\varGamma}\boldsymbol{\cdot}\boldsymbol{x} + \boldsymbol{\mathsf{R}}^{T}(\theta_j,\phi_j)\cdot\sum_{i=1}^{5}\boldsymbol{u}_{is}(\boldsymbol{x'};\theta_j,\phi_j), \end{gather}
\begin{gather}\frac{\textrm{d}\phi_j}{\textrm{d}t}={-}\frac{1}{\kappa^{2}+1}(\kappa^{2}\sin^{2}{\phi_j}+\cos^{2}{\phi_j}), \end{gather}
\begin{gather}\frac{\textrm{d}\theta_j}{\textrm{d}t}=\frac{1}{4}\left(\frac{\kappa^{2}-1}{\kappa^{2}+1}\right)\sin{2\theta_j}\sin{2\phi_j}. \end{gather}
Here, 
$\boldsymbol{\mathsf{R}}(\theta _j,\phi _j)$ is the rotation tensor relating the space and body-fixed coordinates.
3. A spheroid in the spinning or log-rolling orbit (
$C=0$)
To begin with, we analyse the simpler case of a spinning (log-rolling) spheroid. This also serves as a partial validation for our numerical integration, based on prior knowledge of the streamline patterns for the limiting cases of a cylinder and a sphere. The spinning (log-rolling) orbit corresponds to 
$\theta _j = 0$, and is obtained from (2.1a,b) by setting the orbit constant 
$C = 0$. For this configuration, the only non-zero component in the general disturbance field is 
$\boldsymbol {u}_{3s}$, corresponding to a transverse planar extension with the extensional axis oriented at 45 degrees in the flow–gradient plane.
Numerically integrating (2.11) with only the disturbance velocity field 
$\boldsymbol {u}_{3s}$ one may obtain the streamlines around a spinning spheroid. Figures 3 and 4 show the topology of the streamlines both within and outside of the flow–gradient plane. Figure 3 depicts the in-plane streamline patterns for a sphere, a cylinder and a prolate spheroid with 
$\xi _0=1.05 (\kappa = 3.28)$, the qualitative resemblance between the patterns being readily evident. The only perceptible difference is an increase in the size of the closed-streamline region as one moves from the sphere to the cylinder, this being consistent with the increasing strength of the velocity disturbance, as is also evident from the far-field behaviour of the disturbance field in the two cases – the 
$O(1/r^{2})$ decay for the sphere as opposed to the 
$O(1/r)$ decay for the circular cylinder. Figure 4 shows the off-plane streamlines, which intersect the gradient–vorticity plane at 
$z = 0.5$. The three-dimensional view (figure 4c), and the projections of these streamlines on the flow–gradient (figure 4a) and gradient–vorticity (figure 4b) planes are included. The separatrix surface is non-axisymmetric for all cases except that of a sphere, and this is evident from the differing projections of the separatrix envelope in the aforementioned planes; note that the envelope in the flow–gradient plane is a single streamline, but that in the gradient–vorticity plane is constructed from multiple separatrix streamlines (also shown in red in figure 4b). We have verified that, although non-axisymmetric, the separatrix surface extends to infinity in the vorticity direction for all 
$\kappa$.
The in-plane fore–aft symmetric streamline pattern for (a) a sphere, (b) a cylinder and (c) a prolate spheroid with 
$\kappa =3.28 (\xi _0=1.05)$. In each case, fore–aft open streamlines are separated from the closed ones by separatrices (shown in red).
The three-dimensional streamline pattern around a log-rolling prolate spheroid with 
$\kappa = 3.28 (\xi _0 = 1.05)$ corresponding to 
$z = 0.5$ (on the gradient–vorticity plane). The panels include projections on the flow–gradient plane (a), the gradient–vorticity plane (b) and a three-dimensional view (c). The projection of the separatrix surface, and the constituent separatrix streamlines, in each of these figures are shown in red. In figure 4(c), black dots denote the intersection of the separatrix envelope with the gradient–vorticity plane.
Figures 5(a) and 5(b) show the shapes of both an open and a closed streamline around log-rolling prolate spheroids with different aspect ratios (
$\kappa =10.04,3.25,1.67,1.15$); also included are limiting cases of a sphere (
$\kappa = 1$) and a cylinder (
$\kappa = \infty$). The open streamline deviates the most (least) from the corresponding ambient streamline, and the closed streamline is largest (smallest) in size, for 
$\kappa = \infty (1)$. Figure 5(c) shows the time periods for closed streamlines, starting from a given location on the flow axis, increasing monotonically as one goes from a cylinder to a sphere. The limiting values for a sphere and cylinder are available for earlier efforts (Torza et al. Reference Torza, Henry, Cox and Mason1971; Powell Reference Powell1983). All of these features are again consistent with the larger disturbance field associated with longer prolate spheroids (the cross-sectional diameter is fixed as 
$\kappa$ varies), with the largest disturbance corresponding to a cylinder. In summary, the streamline topology associated with a spinning spheroid is identical to that around a sphere and a circular cylinder, with the various measures associated with the streamline pattern around a log-rolling prolate spheroid, including the time periods of the closed streamlines, the width of the closed streamline region etc., bounded between those for a sphere and a cylinder.
The shapes of the open (a) and closed (b) streamlines for prolate spheroids with 
$\kappa (\xi _0)=10.04(1.005),3.25(1.05),1.67(1.25),1.15(2)$; the limiting cases of a sphere (
$\kappa =1$) and a cylinder (
$\kappa =\infty$) are shown. The open streamlines coincide at upstream and downstream infinity, while the closed streamlines start from the same location on the flow axis (
$-3,0.0$); (c) A comparison of closed-streamline time periods for the aforementioned aspect ratios, including the limiting sphere and cylinder cases.
Finally, figure 6 shows the separatrices for log-rolling prolate spheroids of different aspect ratios, and compares their far-field forms with analytical predictions (shown in dotted lines) derived in appendix A. For sufficiently large aspect ratios (
$\kappa \rightarrow \infty$ or 
${\xi _0 \rightarrow 1}$), the separatrix exhibits an intermediate asymptotic form, 
$(1/\sqrt {2}) x^{-1}$, in the interval 
$1 \ll x \ll \kappa$, corresponding to that of a circular cylinder, before eventually transitioning to the more rapid 
$O(x^{-3/2})$ decay for 
$x \gg \kappa$.
Separatrices in the flow–gradient plane, for spinning prolate spheroids (
$\kappa =100.04,31.63,10.04,3.28$) in a simple shear flow. The dashed lines denote the analytical far-field asymptotes derived in appendix A.
4. A spheroid in non-spinning orbits ($C>0$)
4.1. Regular open pathlines
The results presented thus far show that the streamline pattern around spinning spheroids is steady, and analogous to that around a sphere. In this section, we consider spheroids in orbits other than the spinning one (
$C > 0$). For all such cases, the velocity field is time dependent. The focus from hereon will therefore be on the fluid pathlines rather than the (time periodic) sequence of instantaneous streamline patterns associated with the rotating spheroid. In our examination of fluid pathlines, it is natural to start from a spheroid tumbling in the flow–gradient plane. The polar Jeffery angle is now a constant (
$\theta _j = {\rm \pi}/2$), with the time dependence of the velocity field arising only due to 
$\phi _j$ varying with time. Thus, one only needs to integrate (2.11) and (2.12) in order to characterize the pathlines. Further, the disturbance velocity field only involves three of the five canonical components in the disturbance velocity fields in table 1 viz. the axisymmetric extension (
$\boldsymbol {u}_{1s}$), the first transverse extension (
$\boldsymbol {u}_{2s}$) and the second longitudinal extension (
$\boldsymbol {u}_{5s}$). An additional symmetry-induced simplification is that fluid elements in the flow–gradient plane are confined to this plane for all times, despite the time dependence of the flow field. This makes it easier to visualize the changes in the pathline topology relative to the case of a spheroid rotating in a generic precessional orbit.
The changing orientation of the tumbling spheroid implies that, unlike the spinning case, one needs to specify an ‘initial’ orientation in order to fix the fluid pathlines. A convenient choice is the spheroid orientation at the instant that the fluid element crosses the gradient–vorticity plane (
$x =0$); the coordinates of the element at this instant may be taken as (
$y_0,z_0$). This orientation is fixed by the azimuthal Jeffery angle 
$\phi _{j0}$. The fore–aft symmetry of the spheroid implies that 
$\phi _{j0}$ may be restricted to the interval (
$0,{\rm \pi}$). Further, from reversibility considerations, one only need consider 
$\phi _{j0}$ values in the interval (
$0,{\rm \pi} /2$), since the pathline configuration for 
${\rm \pi} -\phi _{j0}$ may be obtained from that for 
$\phi _{j0}$ via a reflection transformation with respect to the gradient–vorticity plane. Each 
$\phi _{j0}$ in (
$0,{\rm \pi} /2$) leads to a distinct configuration of pathlines, an individual pathline being obtained by integrating forward and backward in time starting from the initial position of the fluid element in the gradient–vorticity plane. The pathline configurations for fluid elements in the flow–gradient plane, and for five different values of 
$\phi _{j0}$ in (
$0,{\rm \pi} /2$), are shown in figure 7. While the pathlines for the cases 
$\phi _{j0} = 0$ and 
${\rm \pi} /2$ are evidently fore–aft symmetric, all other values of 
$\phi _{j0}$ lead to fore–aft asymmetric open pathlines that suffer a net displacement (
$\Delta y$) in the gradient direction as they head from upstream to downstream infinity. By analogy with the spinning spheroid, the ‘limiting’ open pathlines are again shown in red for all configurations in figure 7. The expectation is that these red curves now serve as fore–aft asymmetric separatrices (for 
$\phi _{j0} \neq 0,{\rm \pi} /2$), asymptoting to the flow axis far upstream or downstream. These separatrices would therefore seem to demarcate pathlines that extend to infinity in the upstream and downstream directions from those that do not, being forced to cross the flow axis at a finite 
$x$, either positive or negative (in a manner similar to the spinning spheroid or sphere above). It will be seen in § 5, where we consider pathlines that loop around the spheroid as a function of their upstream coordinates, that a subset of these separatrices mark the onset of chaotic scattering. This subset must correspond to the interval 
$\phi _{j0} \in ({\rm \pi} /2,{\rm \pi} )$ since, as implied by figure 7, it is separatrices with 
$\phi _{j0}$ values in the above interval that begin from finite offsets at negative 
$x$ values, and are therefore accessible from upstream infinity.
The configuration of the in-plane regular pathlines for a tumbling prolate spheroid with 
$\kappa =3.28 (\xi _0=1.05)$. The different configurations correspond to 
$\phi _{j0}$ values of (a) 
$0$, (b) 
${\rm \pi} /6$, (c) 
${\rm \pi} /4$, (d) 
${\rm \pi} /3$ and (e) 
${\rm \pi} /2$. All configurations are bounded below by separatrices shown in red. The 
$\Delta y$ values marked in (b) and (d) denote the lateral displacement in the gradient direction.
Figure 8 shows the off-plane pathline configurations, again for a tumbling prolate spheroid. As for the in-plane case, the pathlines correspond to fluid elements that cross the gradient–vorticity plane (
$x= 0$) at a fixed vorticity offset of 
$z_0 = 0.25$; the 
$\phi _{j0}$ values are the same as in figure 7. Expectedly, the off-plane pathlines are three-dimensional trajectories, and their fore–aft asymmetry manifests as net displacements in both the gradient and vorticity directions; these displacements (
$\Delta y$ or 
$\Delta z$) are indicated in the relevant pathline projections. Note that the ‘separatrices’ in particular have small-scale wiggles superposed on a slower large-scale variation, on account of the short-time tumbling dynamics of the spheroid.
The configurations of the off-plane (vorticity offset 
$z = 0.25$) regular open pathlines projected onto the flow–gradient (a–e) and flow–vorticity (f–j) planes for a tumbling prolate spheroid with 
$\kappa =3.28(\xi _0=1.05)$; off-plane separatrices are shown in red. The different configurations correspond to 
$\phi _{j0}$ values of (a,f) 
$0$, (b,g) 
${\rm \pi} /6$, (c,h) 
${\rm \pi} /4$, (d,i) 
${\rm \pi} /3$ and (e,j) 
${\rm \pi} /2$. Here, 
$\Delta y$ and 
$\Delta z$ denote the lateral displacements in the gradient and vorticity directions, respectively.
4.1.1. Lateral displacements of the regular open pathlines: $0 < C \leqslant \infty$
The upstream and downstream coordinates, along the gradient and vorticity directions, are identical for all open streamlines associated with the spinning spheroid or a sphere. An obvious measure of fore–aft asymmetry, as already implied in figure 7 and 8, is therefore the difference between the upstream and downstream gradient and vorticity coordinates of a given open pathline. These lateral displacements in the gradient (
$\Delta y$) and vorticity (
$\Delta z$) directions have been used earlier to characterize the asymmetry of finite-
$St$ pair-sphere trajectories (
$St$ here being the Stokes number, and a measure of particle inertia; see Subramanian & Brady (Reference Subramanian and Brady2006) and to calculate the shear-induced diffusivity arising from irreversible pair interactions in spherical particle suspensions (Da Cunha & Hinch Reference Da Cunha and Hinch1996). The definition of the lateral displacement in the gradient direction (
$\Delta y$) is indicated in figure 9(a), along with a schematic in figure 9(b) highlighting the coordinates used for a general off-plane pathline. Both 
$\Delta y$ and 
$\Delta z$ may be calculated numerically, in a straightforward manner, for all pathlines down to the ‘separatrices’ shown in figure 7 and 8. For pathlines that approach the spheroid to within a distance of order its major axis, the lateral displacements will be comparable to the respective upstream coordinates: that is, 
$\Delta y \sim O(y_{-\infty })$ and 
$\Delta z \sim O(z_{-\infty })$. For far-field pathlines, however, the displacements are much smaller, allowing for an analytical approach involving a perturbation about the upstream coordinates. The analysis, detailed in appendix B, leads to the following expressions for the lateral displacements associated with the far-field pathlines:
\begin{align} {\rm \Delta} y &= \frac{K_{ap/o}}{8y_{-\infty}^{2}}\int_{-\infty}^{\infty}\frac{\textrm{d}t}{(t^{2}+\sec^{2}\varTheta)^{5/2}}[\sin ^{2}\theta_j (\tan ^{2}\varTheta (3 \cos 2\theta_j+1) \sin 2\phi_j \nonumber\\ &\quad +\,6 \tan \varTheta \sin 2\theta_j \sin 2\phi_j (\sin \phi_j+t \cos \phi_j)\nonumber\\ &\quad +\,12 \sin ^{2}\theta_j \sin \phi_j \cos \phi_j (t ^{2} \cos ^{2}\phi_j+\sin ^{2}\phi_j)\nonumber\\ &\quad +6 t \sin ^{2}\theta_j \sin ^{2}2\phi_j-2 (t ^{2}+1) \sin 2\phi_j)] \nonumber\\ &\quad +\frac{K_{tp/o}}{8y_{-\infty}^{2}}\int_{-\infty}^{\infty}\frac{\textrm{d}t}{(t^{2}+\sec^{2}\varTheta)^{5/2}}[\sin 2\theta_j (({-}2 \tan ^{2}\varTheta+t ^{2}+1) \sin 2\theta_j \sin 2\phi_j\nonumber\\ &\quad +4 \tan \varTheta (t \sin \phi_j+\cos \phi_j))-8 \tan \varTheta \sin ^{3}\theta_j \cos \theta_j \sin 2\phi_j (\sin \phi_j+t \cos \phi_j)\nonumber\\ &\quad +\sin ^{4}\theta_j (2 t \cos 4\phi_j-(t ^{2}-1) \sin 4\phi_j)+t \sin ^{2}\theta_j (\cos 2\theta_j+3)] \nonumber\\ &\quad + \frac{K_{lp/o}}{8y_{-\infty}^{2}}\int_{-\infty}^{\infty}\frac{\textrm{d}t}{(t^{2}+\sec^{2}\varTheta)^{5/2}}[\sin ^{2}\theta_j \sin 2\phi_j (\tan ^{2}\varTheta (\cos 2\theta_j+3)-2 (t ^{2}+1))\nonumber\\ &\quad +4 \tan \varTheta \sin ^{3}\theta_j \cos \theta_j \sin 2\phi_j (\sin \phi_j+t \cos \phi_j)-4 \tan \varTheta \sin 2\theta_j (t \sin \phi_j+\cos \phi_j)\nonumber\\ & \quad +2 \sin ^{4}\theta_j \sin 2\phi_j (\sin \phi_j+t \cos \phi_j)^{2}+8 t \cos ^{2}\theta_j], \end{align}
\begin{gather} \Delta z=\Delta y\tan\varTheta, \end{gather}
where 
$K_{ap/o}$, 
$K_{tp/o}$ and 
$K_{lp/o}$ are aspect-ratio-dependent stresslet coefficients defined in table 2. The integral in (4.1) is evaluated numerically, and is a function of the spheroid aspect ratio (
$\kappa$), the orbit constant (
$C$) and 
$\phi _{j0}$; the case of a tumbling spheroid may be obtained by setting 
$\theta _j = {\rm \pi}/2(C = \infty )$. On account of the radial nature of the dipolar field, the dominant displacement of the pathline occurs when the fluid element is at a distance of O(
$y_{-\infty }$) or smaller from the gradient–vorticity plane. Since this corresponds to a time interval of O(
$\dot {\gamma }^{-1}$), independent of 
$y_{-\infty }$, the far-field lateral displacements exhibit a decay of (
$y_{-\infty }^{-2}$) which is the same as that of the dipolar field. In light of this decay, in figure 10, we plot the rescaled lateral displacements, 
$y_{-\infty }^{2}\Delta y$, for pathlines in the flow–gradient plane, for tumbling prolate and oblate spheroids. Owing to the simple relationship between the two lateral displacements in (4.1), in figure 11 we plot the re-scaled total lateral displacement, 
$y_{-\infty }^{2} \Delta r = y_{-\infty }^{2}\sqrt {\Delta y^{2} + \Delta z^{2}} = y_{-\infty }^{2}\Delta y\sec \varTheta$, for pathlines outside the flow–gradient plane (corresponding to 
$\varTheta = {\rm \pi}/9$), again for tumbling prolate and oblate spheroids. Figure 12 shows the re-scaled total lateral displacements for off-plane fluid elements (
$\varTheta ={\rm \pi} /9$), and for prolate and oblate spheroids in a non-tumbling orbit (
$C = 10$). Note that for this generic case, it is necessary to consider 
$\phi _{j0}$ in the larger interval (
$0,{\rm \pi}$), with the displacements in the interval (
${\rm \pi} ,2{\rm \pi}$) again obtained via a reflection transformation. In all cases (figure 10, 11 and 12), the numerical displacement curves approach the analytical far-field ones in the limit 
$(y_{-\infty }^{2} + z_{-\infty }^{2})^{1/2} \rightarrow \infty$.
Panel (a) shows the lateral displacement of a regular open fluid pathline in the flow–gradient plane. In (b) 
$\varTheta$ characterizes the initial upstream location of a fluid element outside the flow–gradient plane.
Lateral displacements in the gradient direction for fluid elements in the flow–gradient plane, as a function of the initial spheroid orientation (
$\phi _{j0}$). (a) A prolate spheroid (
$\xi _0=1.01(\kappa =7.12$)), and (b) an oblate spheroid (
$\xi _0=1.01(\kappa =0.14$)); both spheroids are in the tumbling orbit. The lateral displacement for 
$\phi _{j0} \in ({\rm \pi} /2,{\rm \pi} )$ may be obtained as 
$\Delta y(\phi _{j0}) = -\Delta y({\rm \pi} -\phi _{j0})$.
The total lateral displacements for off-plane fluid elements (
$\varTheta ={\rm \pi} /9$) as a function of the initial spheroid orientation (
$\phi _{j0}$). (a) A prolate spheroid (
$\xi _0=1.01(\kappa =7.12$)), and (b) an oblate spheroid (
$\xi _0=1.01(\kappa =0.14$)); both spheroids are in the tumbling orbit. The lateral displacement for 
$\phi _{j0} \in ({\rm \pi} /2,{\rm \pi} )$ may be obtained as 
$\Delta y(\phi _{j0}) = -\Delta y({\rm \pi} -\phi _{j0})$.
The total lateral displacements for off-plane fluid elements (
$\varTheta ={\rm \pi} /9$) as a function of the initial spheroid orientation (
$\phi _{j0}$). (a) A prolate spheroid (
$\xi _0=1.01(\kappa =7.12$)), and (b) an oblate spheroid (
$\xi _0=1.01(\kappa =0.14$)). Both spheroids are in a Jeffery orbit with 
$C=10$. The lateral displacement for 
$\phi _{j0} \in ({\rm \pi} ,2{\rm \pi} )$ may be obtained as 
$\Delta y(\phi _{j0}) = -\Delta y(2{\rm \pi} -\phi _{j0})$.
The aspect-ratio-dependent stresslet coefficients that appear in the far-field stresslet disturbance (see (4.1) in the text).
4.1.2. The separatrix of a tumbling spheroid
We end this subsection on the regular pathlines by analysing the approach of the separatrices, marked in red in figure 7, towards the flow axis – the approach occurs in both the upstream and downstream directions for 
$\phi _{j0} = 0,{\rm \pi} /2$, and in either of these two directions for other 
$\phi _{j0}$ values. For the spinning spheroid, as seen from the analytical asymptote in figure 6, the rate of approach of the far-field separatrix (
$y \propto x^{-3/2}$) remains identical to that of a sphere, with only a pre-factor that is a function of aspect ratio. The analysis for the analogously defined separatrices in figure 7 is more involved owing to the spheroid rotation. In the far field (
$r \gg 1$), one may, however, exploit a separation of time scales, with the fast scale corresponding to the Jeffery period (of 
$O(\dot {\gamma }^{-1})$ for 
$\kappa \sim O(1)$) and the slow time scale characterizing the rate of change of position of the fluid element. Note that, unlike the regular pathlines for large 
$y_{-\infty }$ analysed above, where this latter time scale remains 
$O(\dot {\gamma }^{-1})$, for the separatrices, the approach towards the flow axis implies that 
$y$ is asymptotically small and the time scale of convection by the ambient shear is therefore correspondingly large. Note also that the magnitude of this slow time scale, being a function of 
$y$, is not known a priori. Keeping this in mind, we write the governing equations as
\begin{gather} \frac{\textrm{d}\boldsymbol{x}}{\textrm{d}T}=\boldsymbol{u}(\boldsymbol{x},\boldsymbol{p}(\tau)), \end{gather}
\begin{gather}\frac{\textrm{d}\boldsymbol{p}}{\textrm{d}\tau}=\boldsymbol{\omega}\boldsymbol{\cdot} \boldsymbol{p}+\frac{\kappa^{2}-1}{\kappa^{2}+1}\left({\boldsymbol{\mathsf{E}}}\boldsymbol{\cdot} \boldsymbol{p}-\boldsymbol{p}\left({\boldsymbol{\mathsf{E}}}\boldsymbol{\colon}\boldsymbol{pp}\right)\right), \end{gather}
where we have (formally) identified 
$\tau$ and 
$T$ as the fast and slow time variables, respectively. The identification implies that 
$\boldsymbol {p}$ only evolves on the fast time scale, while 
$\boldsymbol {x}$ is a function of both 
$\tau$ and 
$T$. One may expand 
$\boldsymbol {x}$ in the usual manner 
$\boldsymbol {x}(T,\tau )=\boldsymbol {X}(T,\tau )+\tilde {\boldsymbol {x}}(\boldsymbol {X},T,\tau )$ with the fast variable satisfying 
$1/T_j \int _0^{T_j}\tilde {\boldsymbol {x}}\,\textrm {d}\tau =0$, with 
$T$ fixed, and 
$T_j$ being the Jeffery period. The slow and fast contributions of 
$\boldsymbol {x}$ may further be expanded in the following manner: 
$\boldsymbol {X}(T,\tau )=\boldsymbol {X}_0(T) + \boldsymbol {X}_1(T,\tau )+\cdots$ and 
$\tilde {\boldsymbol {x}}= \boldsymbol {\tilde {x}}_0+\boldsymbol {\tilde {x}}_1+\cdots$, respectively, where we have anticipated the leading-order position to only evolve on the slow time scale. Since there is not a readily available estimate of the slow time scale a priori, we do not introduce 
$\varepsilon$ above, but nevertheless assume the higher-order terms in the above expansions to be asymptotically smaller; thus, 
$\tilde {\boldsymbol {x}} \ll \boldsymbol {X}, \boldsymbol {X}_1 \ll \boldsymbol {X}_0$, and so on. We restrict ourselves to calculating the leading-order slow contribution 
$\boldsymbol {X}_0(T)$ below (although, note that 
$\tilde {x}_0$ is responsible for the small-scale wiggles on the separatrices seen in figure 8 and 13.)
The far-field form may now be obtained by recalling the approximations used in the context of the spinning spheroid in appendix A. Thus, in the far field at leading order, we only retain the ambient flow in the 
$x$-component, and the octupolar contribution in the 
$y$-component. One may now write the scalar components of the two-time-scale equations as
\begin{gather} \frac{\textrm{d}X_0}{\textrm{d}T}+\frac{\partial\tilde{x}_0}{\partial \tau}=Y_0+\tilde{y}_0\end{gather}
\begin{gather} \frac{\textrm{d}Y_0}{\textrm{d}T}+\frac{\partial\tilde{y}_0}{\partial \tau}={-}\frac{15}{4}d^{5}\xi_0^{3}L^{2}_c\kappa^{2}(\kappa^{2}-1)\left[(K_{ap/o}-2K_{tp/o}+K_{lp/o})\kappa^{2}\vphantom{\left(\frac{\kappa\tau}{\kappa^{2}+1}\right)}\right. \nonumber\\ \quad \left. +\, 2(K_{tp/o}-K_{ap/o})\cot^{2}\left(\frac{\kappa\tau}{\kappa^{2}+1}\right)\right] \left[\kappa^{2}+\cot^{2}\left(\frac{\kappa\tau}{\kappa^{2}+1}\right)\right]^{{-}3} X_0^{{-}4}, \end{gather}
where 
$L_c=\overline {\xi _0}(\xi _0)$ for a prolate (an oblate) spheroid. Averaging (4.6) and (4.7) over 
$\tau$, with 
$T$ fixed, and using the periodicity of the fast contributions, one obtains the equations governing 
$X_0$ and 
$Y_0$ as
\begin{gather} \frac{\textrm{d}X_0}{\textrm{d}T} = Y_0 \end{gather}
\begin{gather} \frac{\textrm{d}Y_0}{\textrm{d}T} ={-}\frac{d^{5}\xi_0^{3}L_c^{2}}{32X_0^{4}(1+\kappa)^{2}} \kappa[K_{ap/o} ({-}3 - 30\kappa + 9 \kappa^{2}) \nonumber\\ \quad +\, K_{lp/o} (7 - 5 \kappa(2 + \kappa)) - 4 K_{tp/o} (5 + \kappa ({-}2 + 5 \kappa))]. \end{gather}The analytical far-field form may now be obtained by taking the ratio of (4.7) and (4.8). Reverting to the original notation, one obtains the far-field separatrices (without the superimposed small-scale wiggles) in the following form:
\begin{align} y_{p/o}=4\kappa^{f_{p/o}}\sqrt{\frac{(K_{ap/o} ({-}9 \kappa^{2}\!+\!30 \kappa \!+\!3)+ 4 K_{tp/o} (\kappa (5 \kappa\!-\!2)\!+\!5)+K_{lp/o} (5 \kappa (\kappa +2)-7))}{3(\kappa +1)^{2}x^{3}}} \end{align}
again exhibiting a far-field decay of 
$x^{-3/2}$, with 
$f_{p/o}=-1/2(1/2)$. Figure 13 shows the comparison between both fore–aft symmetric (for 
$\phi _{j0} = {\rm \pi}/2$) and asymmetric (for 
$\phi _{j0} = {\rm \pi}/3({\rm \pi} /6$ for oblate)) separatrices obtained numerically, and the analytical predictions above. The numerical separatrices are obtained by choosing an initial location (
$y_0$) on the gradient axis, and progressively refining this based on the final outcome of the integrated trajectory (escape to 
$x \rightarrow \pm \infty$, or crossing of the flow axis at a finite 
$x$). Notwithstanding the small-scale wiggles that decrease in amplitude with increasing 
$x$, a good comparison is obtained with (4.9).
Comparison between the numerical and analytical separatrix branches for a tumbling prolate spheroid with 
$\xi _0=1.15 (\kappa =2.025)$: (a) 
$\phi _{j0} = {\rm \pi}/2$ (fore–aft symmetric); (b) 
$\phi _{j0} ={\rm \pi} /3$ (fore–aft asymmetric) and for an oblate spheroid with 
$\xi _0=1.15(\kappa =0.49)$; (c) 
$\phi _{j0}={\rm \pi} /2$ (fore–aft symmetric), (d) 
$\phi _{j0}={\rm \pi} /6$ (fore–aft asymmetric). For the asymmetric cases, one of the separatrix branches approaches a finite gradient offset at infinity.
4.2. The singular pathlines
Thus far, we have dealt with the regular open pathlines, which are (except for 
$\phi _{j0} = 0,{\rm \pi} /2$) fore–aft asymmetric generalizations of the open streamlines around a sphere. That these do not constitute the entire set of pathlines may be seen from the fact that all of the in-plane separatrices in figure 7 intersect the gradient axis at an ordinate value (say, 
$y^{sep}_0$) greater than 
$y_{min}$, the value corresponding to the intersection of the spheroid surface with 
$x = 0$; the 
$y_{min}$ value depends on 
$\phi _{j0}$, and ranges between 
$\bar {\xi }_0/\xi _0$ and unity, these limiting values corresponding to flow and gradient-aligned prolate spheroids, respectively. In-plane pathlines spanning the interval (
$y_{min},y^{sep}_0$) have therefore not been accounted for; a similar argument may be made for the off-plane pathlines. The pathlines in these intervals are the analogues of closed streamlines for the sphere case, and the naive expectation is perhaps that the strictly periodic closed streamlines for a sphere would transform to aperiodic but bounded pathlines for a non-spinning spheroid. Interestingly, as will be seen below, this is not the case.
To examine pathlines in the interval (
$y_{min}, y^{sep}_0$), we consider the case 
$\phi _{j0} = {\rm \pi}/2 (y_{min} = 1)$. The regular (in-plane) pathline configuration analogous to that shown in figure 7, but for a prolate spheroid with 
$\xi _0=2(\kappa = 1.15)$, is shown in figure 14(a). The regular pathlines are bounded below by a fore–aft symmetric separatrix that asymptotes to the flow axis at upstream and downstream infinity, and intersects the gradient axis at 
$y_0^{sep} \approx 1.11621$. As seen in figure 14(b–e), the pathlines ‘below’ the separatrix that intersect the gradient axis in the interval (
$1,y_0^{sep}$) are again open, in that they eventually asymptote to finite 
$y$ values in the upstream and downstream directions, but loop around the spheroid a certain number of times before doing so. Importantly, there appears to be no pattern to the number of loops; this number varies in a seemingly random (but sensitive) fashion as 
$y_0$ is decreased below 
$y^{sep}_0$ even by a very small amount. For instance, the number of loops for the pathlines in figure 14(b) and 14(c) is 2, but that in figure 14(d) is 23, while that in figure 14(e) is again 2; the dramatic increase in the number of loops in figure 14(d) occurs despite a change in 
$y_0$ that is 
$O(10^{-5})$! The fact that the pathlines crossing the gradient axis below the separatrix nevertheless asymptote to finite gradient offsets at upstream and downstream infinity obviously implies a crossing of the separatrix at some point. This crossing does not violate uniqueness since the spheroid orientation, at the time of crossing, is different for the separatrix and any of these looped trajectories (which we term the singular open pathlines); an argument along the same lines also shows that such looped trajectories, asymptoting to different upstream and downstream gradient offsets, is not inconsistent with the reversibility of the Stokes equations either.
The open pathline configuration for a tumbling prolate spheroid (
$\xi _0 =2.0(\kappa =1.15$)) with 
$\phi _{j0} = {\rm \pi}/2$. (b–e) Correspond to singular open pathlines that cross the gradient axis just below the separatrix (shown in red in graph (a), with 
$y_0=y_0^{sep}$) at four different 
$y_0$ values: (b) 
$y_0 = y_0^{sep}-10^{-5}$; (c) 
$y_0 = y_0^{sep}-4\times 10^{-5}$; (d) 
$y_0 = y_0^{sep}-10\times 10^{-5}$; (e) 
$y_0 = y_0^{sep}-20\times 10^{-5}$.
4.2.1. The residence-time distributions
To characterize this seemingly random dependence of the number of loops of the pathlines above, we define an appropriate residence time (
$\mathcal {D}t$), taken to be the total time that a fluid element spends in a certain neighbourhood of the spheroid. Unless otherwise specified, this neighbourhood is taken to be the 
$x$-interval (
$-20,20$) in the figures to follow; the qualitative nature of the findings reported below is independent of the particular choice of neighbourhood, provided the 
$x$-value chosen is large enough compared to the spheroid dimensions. Figures 15(a) and 15(b) depict the variation in the residence time defined above. The residence times are plotted as a function of 
$y_0$, over a range that brackets 
$y_0^{sep}$, for both a tumbling prolate spheroid (with 
$\phi _{j0} = {\rm \pi}/2$) and a spinning one. In the former case, we see a discontinuous change in the nature of the residence-time curve as one moves across the separatrix. The value of 
$\mathcal {D}t$ exhibits an initial smooth increase for the regular pathlines corresponding to 
$y_0 > y_0^{sep}$, but there is an abrupt transition to a seemingly random dependence for the singular pathlines with 
$y_0 < y_0^{sep}$. Figure 15(b) shows the analogous residence-time plot for the spinning spheroid, in which case the residence time varies in a smooth manner with 
$y_0$, attaining a maximum at the closed streamline (corresponding to 
$y'_0$ in the inset) that exactly spans the 
$x$-interval under consideration. The subsequent decrease for 
$y_0 < y_0'$ corresponds to choosing 
$\mathcal {D}t$ as half the period of the closed streamlines that no longer span the chosen 
$x$-interval, but instead intersect the flow axis at a pair of points, within the interval 
$(-20,20)$, and symmetrically placed about the origin.
Variation of the residence time of a fluid element for a prolate spheroid (
$\xi _0=2.0(\kappa =1.15$)): (a) in the tumbling orbit with 
$\phi _{j0} = {\rm \pi}/2$, (b) in the spinning mode; the ordinate of the separatrix 
$y_0^{sep}$ is shown by the dashed line; 
$y_0'$ corresponds to the closed orbit spanning the 
$x$-interval (
$-20,20$).
We now examine the irregular variation of the residence time encountered in figure 15(a) in greater detail. For this purpose, we consider the same spheroid as in figure 14, again in its tumbling orbit (
$\theta _j = {\rm \pi}/2$). The near-unity aspect ratio (
$\xi _0=2; \kappa = 1.15$) chosen serves the additional purpose of illustrating the singular nature of the spherical particle limit from the perspective of the streamline/pathline topology. Figure 16(a) shows the variation of the residence time of the singular open pathlines which pass through the negative 
$x$-axis within the interval 
$x_0 \in (-9.5,-3)$. The residence time is obtained in the manner defined above, via forward and backward time integrations starting from a specific initial point in the aforementioned interval, and until the pathline leaves the interval 
$(-20,20)$. Similar to figure 15(a), the extremely irregular dependence of the residence time on 
$x_0$ is readily apparent (note that each of the pathlines intersects the 
$x$-axis multiple times, and in principle, it is possible that a pair of chosen initial conditions end up being part of the same singular pathline, although this should, in principle, occur with an infinitesimal probability; we have verified that this is not the case by monitoring the upstream and downstream coordinates, and ensuring that these are different for each of the initial points chosen). The plots that follow in figure 16(b–e) show the variation in the residence time at progressively finer resolutions by focussing on increasingly small subsets of the 
$x_0$-interval considered above. These figures correspond to 1000 initial conditions in the 
$x_0$-intervals 
$(-4,-3)$, 
$(-3.1,-3)$, 
$(-3.01,3)$ and 
$(-3.001,3)$, respectively. Thus, in going from figure 16(a to 16e), the length scale of the interval under consideration has decreased by four orders of magnitude, but the irregular dependence of the residence time clearly persists down to the smallest scales. This fractal dependence of the residence-time measure is a signature of chaotic scattering, and is suggestive of the existence of a chaotic saddle in the region around the spheroid (Aref & Pomphrey Reference Aref and Pomphrey1980, Reference Aref and Pomphrey1982; Bleher, Ott & Grebogi Reference Bleher, Ott and Grebogi1989; Bleher et al. Reference Bleher, Grebogi and Ott1990; Aref & Stremler Reference Aref and Stremler1999; Faisst & Eckhardt Reference Faisst and Eckhardt2004).
The sensitive dependence of the residence time of a fluid element, for a tumbling prolate spheroid (
$\xi _0=2.0(\kappa =1.13$)) with 
$\phi _{j0}=0$, plotted as a function of the ‘initial condition’ (defined as a point on the negative 
$x$-axis). Panel (a) depicts the dependence in the entire interval 
$(-9.5,-3)$, while panels (b–e) show the dependence of the residence time in the intervals 
$(-4,-3)$, 
$(-3.1,3)$, 
$(-3.01,3)$ and 
$(-3.001,3)$, respectively.
For the hyperbolic case, the chaotic saddle above would include an infinite number of unstable periodic orbits (and bounded aperiodic ones), whose stable and unstable manifolds intersect to form a Cantor-like set. The existence of single-looped periodic orbits may, for instance, be inferred by considering two neighbouring pathlines, integrated over the approximate duration of a single loop, as shown in figure 17(a) and 17(b). The pair of initial points in the figures correspond to 
$x_0=-3.2217088$ and 
$x_0=-3.2214089$ on the negative 
$x$-axis, and one observes the subsequent intersections with the (negative) 
$x$-axis to occur opposite sides of the initial point in the two cases, implying that the pair of pathlines shown bracket a single-loop closed orbit between them. This argument may be extended in an obvious manner to show the existence of closed pathlines with multiple loops. Note that, from the perspective of a Poincaré map (which gives the dynamics strobed at intervals of 
$T_j$), these closed pathlines are examples of bounded aperiodic orbits since their periods are not integer multiples of 
$T_j$.
Existence of periodic pathlines for the tumbling prolate spheroid (
$\xi _0=4.0(\kappa =1.03$) for a single-loop periodic orbit with an approximate period of 
$40 T_j$: (a) final point (
$x_f$) displaced to the left of the initial point 
$(-3.2214088,0)$ and (b) final point displaced to the right of the initial point 
$(-3.2214089,0)$.
Having established the signatures of chaotic scattering in the residence time plot for a tumbling spheroid, in figure 18(a–d), we consider analogous residence-time plots for the spheroid in other finite-
$C$ precessional orbits. The spheroid is the same as that in figure 16, and the orbits chosen correspond to 
$C = \infty$ (already shown), 5, 0.5, 0.2 and 0.05, respectively, with the initial point again restricted to the interval 
$(-9.5,-3)$ on the negative 
$x$-axis. The plots appear to show a progressive decrease in the extent of irregularity with decreasing 
$C$. One reason for this decrease is that, for lower 
$C$ values, there appear regular intervals between those that correspond to chaotic scattering, and this interlacing behaviour will be seen in more detail in the next section; note that this interlacing tendency is, in fact, already present for 
$C = \infty$, although not evident owing to the much smaller scale (see the inset of figure 18a) which presents a magnified view). The decrease in irregularity is also because, for the smaller 
$C$ values, an increasing fraction of pathlines do not open for the duration of the numerical integration (
$t_{max} = 8000 T_j$); examples of such bounded pathlines are shown in figure 18(d,e). Such bounded pathlines, arising as an artifact of the finite integration duration, begin to populate an upper horizontal plateau with residence times equal to the maximum integration time; in fact, figure 18(d) for 
$C$ = 0.05 consists entirely of pathlines that do not open over the duration of the integration. However, as shown in figure 19, the fraction of initial conditions leading to bounded pathlines for the case 
$C = 0.2$ decreases perceptibly for an integration duration (
$t_{max} = 32\,000 T_j$) four times the one chosen for figure 18. Thus, the fraction of bounded pathlines appears to depend on the duration of the numerical integration. Although impossible to ascertain by numerical means, the limiting integrable case of a spinning spheroid (
$C$ = 0) appears to be a singular one, with a spheroid rotating in an orbit with 
$C$ small but finite nevertheless leading to chaotic scattering; the associated large residence times are, however, inaccessible to numerics.
Sensitive dependence of the residence time of fluid pathlines, on initial conditions, for a prolate spheroid (
$\xi _0=2.0(\kappa =1.13$)) in finite-
$C$ precessional orbits; (a) 
$C=\infty$, (b) 
$C=5.0$, (c) 
$C=0.5$, (d) 
$C=0.05$ and (e) 
$C=0.2$. The insets in (e) shows that the highest residence times (plateau points) correspond to bounded pathlines, while those off the plateau lead to pathlines that open up over times shorter than the integration time (
$8000 T_j$); the insets in (d) confirm that all of the pathlines for this case are bounded.
The residence time of fluid pathlines for a prolate spheroid (
$\xi _0=2(\kappa =1.13$)), in a Jeffery orbit corresponding to 
$C=0.2$, for two different integration times: (a) 
$t_{max}=8000T_j$, (b) 
$t_{max}=32\,000T_j$. The fraction of singular open pathlines increases significantly for 
$t_{max} = 32\,000 T_j$; sample pathlines are shown for points on the plateau as well as for ones off it.
4.2.2. Local graphs of the invariant manifolds and the chaotic saddle
Graphical representations of the unstable and stable manifolds of the chaotic saddle are another qualitative signature of chaos. Here, we numerically evolve an initially circular blob of ‘dye’ located along the negative 
$x$-axis on account of being advected by the fluid motion induced by the rotating spheroid. The numerical protocol is motivated by an experiment demonstrating chaotic scattering, albeit in a different geometric setting (Gouillart et al. Reference Gouillart, Dauchot, Thiffeault and Roux2009). A forward integration in time must trace out the unstable manifold of the chaotic saddle for sufficiently long times, while a backward-in-time integration yields the stable manifold. Figure 20 compares the evolution of the aforementioned circular blob (made up of 
$10^{6}$ randomly distributed initial points) for a sphere and for a (nearly spherical) prolate spheroid with 
$\xi _0 = 4 (\kappa = 1.03)$; the spheroid is chosen to rotate in a Jeffery orbit with 
$C = 20(\phi _{j0} = 0)$. Figures 20(a) and 20(b) show that the blob for the case of a sphere, although significantly distorted for large times owing to differential convection, has, nevertheless a finite spatial extent – the distorted blob lies within the interval 
$(-6,6)$ for all time. This is because the initial blob has been chosen to lie entirely within the region of closed streamlines, and is therefore always bounded by the largest closed streamline that passes through one of the initial points on its periphery. In contrast, figure 20(c) and 20(d) show that the blob, for a spheroid, is sheared out to arbitrarily large distances in the flow direction as time goes to plus and minus infinity on account of the singular open pathlines discussed above (note the differing horizontal scales in figure 20(a,b), and in figure 20(c,d)). While figure 20 highlights the singular role of non-sphericity with regard to the long-time temporal evolution of a blob, which in turn has implications for the rate of mass/heat transport from the spheroid (see conclusions section), it is figure 21 that highlights the striking differences between the sphere and spheroid cases even for finite times. The magnified views in figure 21 clearly show the fine-scaled features (wiggles) that emerge in the evolving blob, only for the spheroid, on account of the convoluted nature of the underlying unstable manifold of the chaotic saddle. In fact, in figure 21(f), the finite initial distribution of points has been sheared out to an extent that the small-scale wiggles associated with the underlying chaotic saddle are no longer well resolved.
The evolution of an initially circular blob of 
$\textrm {radius} = 0.03$, located at 
$(-5,0,0)$ (shown in the inset) for a sphere (a,b) and a prolate spheroid (c,d) with 
$\kappa = 1.03(\xi _0=4.0)$ rotating in a Jeffery orbit with 
$C = 20$: (a,c) forward-in-time integration; (b,d) backward-in-time integration. For (a) 
$t = 2000T_j$, (b) 
$t = -2000T_j$, (c) 
$t = 2000T_j$, (d) 
$t = -2000T_j$.
A comparison of the evolution of an initially circular blob for a sphere and a prolate spheroid (
$\xi _0=4(\kappa =1.03$)) in a precessional orbit with 
$C$ = 20. The evolving blob, for the spheroid, traces out the unstable manifold, of the underlying chaotic saddle, for large times; (a) 
$t=50 T_j$, (b) 
$t=68 T_j$, (c) 
$t=200 T_j$, (d) 
$t=50 T_j$, (e) 
$t=68 T_j$, (f) 
$t=200 T_j$.
Figure 22 shows evolving blobs for spheroids of different aspect ratios, and here, the scale of the aforementioned wiggles is seen to become larger with increasing aspect ratio. Further, the horizontal extent of the sheared blob, for large times, also increases substantially for high aspect ratio spheroids, implying a decrease in the ‘loopiness’ of the singular open pathlines. In the next section we will see a similar behaviour emerge from a different perspective, wherein the dependence of the residence time (of fluid pathlines) is analysed from the point of view of an ‘impact parameter’, defined as the gradient and/or vorticity coordinates of a fluid element far upstream. Finally, in figure 23, we give a representation of the invariant chaotic saddle, its projection onto the flow–gradient plane and in the neighbourhood of the rotating spheroid, obtained in the following manner. The representation corresponds to a near sphere, with 
$\xi _0 = 4 (\kappa = 1.03)$, and it therefore makes sense to relate this representation to an invariant set for a sphere; an example of the latter is a region contained within any of the closed streamlines. With this in mind, we initialize 34 million fluid elements uniformly distributed in a rectangular domain spanning the intervals 
$(-5,5)$ and 
$(-1,1)$ in the 
$x$ and 
$y$-directions, and integrate these over a duration of 
$1000T_j$ both forward and backward in time. Those initial points that remain within a closed streamline that spans the interval 
$(-5,5)$ are plotted in the said figure. The fine-scaled structure evident in figure 23 is a representation of the invariant set responsible for the chaotic scattering signatures detailed above. Further, the emergent non-trivial foliated structure is also indicative of the singular nature of the spherical limit. For a sphere, the invariant set would be the closed streamline above with the region between the streamline and sphere uniformly filled in, since fluid elements within this streamline remain so regardless of the duration of integration. Finally, note that figure 23 denotes (approximately) the intersection of the saddle with the region inside the aforementioned closed streamline. As evident from figure 21 and 22, the saddle will also have a fine-scaled structure ‘radiating’ out of the closed streamline, but resolving this appears to require excessively long integration times.
A comparison of evolving blobs for prolate spheroids in a precessional orbit with 
$C = 20$: (a–c) 
$\xi _0=8.0(\kappa =1.008)$, (d–f) 
$\xi _0=6.0(\kappa =1.01)$ and (g–i) 
$\xi _0=4.0(\kappa =1.03)$. The scale of the wiggles, and the horizontal spread of the distorted blob, increase with increasing aspect ratio. For (a) 
$t=50 T_j$, (b) 
$t=68 T_j$, (c) 
$t=200 T_j$, (d) 
$t=50 T_j$, (e) 
$t=68 T_j$, (f) 
$t=200 T_j$, (g) 
$t=50 T_j$, (h) 
$t=68 T_j$, (i) 
$t=200 T_j$.
A representation of the intersection of the invariant chaotic saddle with the region contained within a closed streamline (for a sphere); the graphs shown are for a prolate spheroid with 
$\xi _0=4(\kappa =1.03)$ in a Jeffery orbit with 
$C=$ 20. (a) Full view; (b) and (c) correspond to magnified views of the left- and right-hand portions, respectively. The dashed curves denote the closed streamline of choice (note that the central part of the saddle is truncated at the top and bottom, being limited by the width of the initial rectangular region).
4.2.3. The uncertainty dimension of the chaotic saddle
The chaotic saddle may be characterized quantitatively in terms of its fractal dimension. While there are quite a few ways by which one may obtain such a dimension, herein we adopt the methodology of Bleher et al. (Reference Bleher, Grebogi and Ott1990) by determining the uncertainty dimension. In this method a large number (
$N_0$) of initial conditions are chosen randomly, and each of these conditions is perturbed by a small amount (
$\epsilon$), and is marked uncertain if the corresponding output states are markedly different. In the present case, we choose 5000 pairs of original and 
$\epsilon$-perturbed upstream gradient offsets in the interval 
$0 \leqslant y_{-\infty } < y^{sep}_{-\infty }$, with 
$x = -20$, 
$z_{-\infty }=0$ and 
$\phi _{j,-\infty } = 0$, the latter corresponding to the spheroid phase, in a particular Jeffery orbit, at 
$x = - 20$. We then check if the integrated pair of final states have similarly signed 
$y$-components, implying that both initial members of the chosen pair asymptote to either upstream (negative 
$y$) or downstream (positive 
$y$) infinity. If not, the pair is deemed uncertain, and one then calculates the uncertainty ratio as 
$f(\epsilon ) = N/N_0$, 
$N$ being the number of uncertain pairs. For the number of initial pairs chosen above, 
$N$ is 
$O(10^{3})$ even for the smallest 
$\epsilon$ values examined, ensuring statistical accuracy. Figures 24(a) and 24(b) show the variation of the uncertainty ratio with aspect ratio for tumbling prolate spheroids, and with orbit constant (
$C$) for a given prolate spheroid with 
$\xi _0 = 1.05 (\kappa = 3.28)$, respectively. An algebraic scaling, 
$f(\epsilon ) \sim \epsilon ^{\alpha }$, for sufficiently small 
$\epsilon$, is observed in all cases over the range of 
$\epsilon$ values examined (down to 
$O(10^{-8})$). The uncertainty dimension is then defined as 
$D_c = D - \alpha$, where we take 
$D = 3$ since fractality is only observed is along the gradient direction; in § 6, it is shown that the non-trivial interlacing of the regular and chaotic scattering intervals persists for all values of the vorticity coordinate. For a spinning spheroid, as already shown (see figure 15b) there is a unique value of the separatrix offset. For this case alone, one expects 
$f(\epsilon ) \propto \epsilon$, so that 
$\alpha = 1$ and 
$D_c = 2$, consistent with the fact that the separatrix is a regular surface embedded in three dimensions. However, as already implied by the results in figure 18, the spinning spheroid appears to be a singular limit, and this is reinforced by the results of figure 24(b), where 
$D_c$ does not approach 2 even as 
$C$ decreases to 0.01. Finally, it needs mention that despite the algebraic scaling obtained in figure 24, as will be seen later, the probabilities densities of residence times show algebraic rather than the exponential tails familiar from much of the literature (Bleher et al. Reference Bleher, Ott and Grebogi1989, Reference Bleher, Grebogi and Ott1990; Jánosi et al. Reference Jánosi, Tél, Wolf and Gallas1997; Faisst & Eckhardt Reference Faisst and Eckhardt2004; Eckhardt et al. Reference Eckhardt, Schneider, Hof and Westerweel2007). The implications of this finding are discussed briefly in the next section.
Plot of the uncertainty ratio 
$f(\epsilon )$ of the chaotic saddle for (a) prolate spheroids of different aspect ratios in the precessional orbit of 
$C$=20 and (b) a prolate spheroid of 
$\xi _0=1.05$ in different Jeffery orbits (
$C=0.01$, 0.07, 0.5 and 1).
5. Transition from regular to chaotic pathlines: boundary of the chaotic saddle
Having established the signatures of chaotic scattering for the fluid pathlines around any non-spinning spheroid, we now proceed to analyse the pathlines from the perspective used in earlier efforts (for instance, see Bleher et al. Reference Bleher, Ott and Grebogi1989, Reference Bleher, Grebogi and Ott1990). This is done by plotting the output variable, which is still the residence time, as a function of an appropriate input one, the focus being on the singular regions corresponding to apparent divergences of the residence time. Such regions would appear scale dependent in a numerical investigation, on account of the true divergences being localized on a Cantor set of vanishing measure (Jung, Tél & Ziemniak Reference Jung, Tél and Ziemniak1993; Ziemniak, Jung & Tél Reference Ziemniak, Jung and Tél1994; Skufca, Yorke & Eckhardt Reference Skufca, Yorke and Eckhardt2006; Eckhardt et al. Reference Eckhardt, Schneider, Hof and Westerweel2007). The input variable must correspond to the analogue of the impact parameter used in the chaotic scattering literature (Bleher et al. Reference Bleher, Grebogi and Ott1990). In our case, this is the upstream gradient offset (
$y_{-\infty }$) of a fluid element in the plane of shear; or, both the gradient and vorticity offsets (
$y_{-\infty }$ and 
$z_{-\infty }$) of an element off the flow–gradient plane. In addition, one needs to specify the particular Jeffery orbit, and the phase of rotation (within a given orbit) at the initial instant; the latter is again done via the azimuthal angle which is now denoted as 
$\phi _{j-\infty }$ (already used in § 4.2.3 above). Note that these upstream offsets have already been used to analyse the regular open pathlines above. The difference here is that the offsets are used as a common basis to analyse both the regular and singular pathlines; the latter group of pathlines was examined earlier, only as a function of their points of intersection with the flow axis, in order to establish the signatures associated with chaotic scattering. One may now analyse the transition from the regular to the singular open pathlines, as a function of the aforementioned upstream offsets, and thereby, analyse the boundary of the chaotic saddle (what has been termed ‘the edge of chaos’ in the literature; see Skufca et al. Reference Skufca, Yorke and Eckhardt2006). The latter has been the subject of earlier investigations in the context of two-dimensional Hamiltonian problems (Jung & Ziemniak Reference Jung and Ziemniak1992).
In figure 25 and 26 we plot the residence time as a function of 
$y_{-\infty } (z_{-\infty } = 0)$ for tumbling prolate spheroids with 
$\xi _0=2.0(\kappa =1.13)$ and 
$\xi _0=1.05(\kappa =3.28$), respectively; 
$\phi _{j,-\infty } = 0$ for both cases. The residence-time curves exhibit a small-scale modulation superposed on an underlying monotonic increase with decreasing 
$y_{-\infty }$. The reason for the latter increase is obvious. The modulation leading to a non-monotonic dependence on 
$y_{-\infty }$ arises owing to the spheroid orientation being a function of time. For 
$y_{-\infty }$ values not too different, a fluid element crossing the gradient–vorticity plane, at the instant that the spheroid is close to a vertical orientation, will do so at a higher 
$y_0$, leading to a reduced residence time. With decreasing 
$y_{-\infty }$, the residence time, for a given change in 
$y_{-\infty }$, changes by a greater amount relative to the fixed Jeffery period. This leads to the wavelength of the modulation in figure 25 and 26 (and the ones thereafter) decreasing with 
$y_{-\infty }$. Note that the absolute wavelength scales in proportion to the Jeffery period, and is therefore greater for the larger aspect ratio spheroid (
$\xi _0$ = 1.05). The red vertical line in each plot denotes the ordinate of the separatrix (
$y^{sep}_{-\infty }$) that demarcates the regular regions from those that exhibit intervals of chaotic scattering; this demarcation is evident from the magnified views in figure 25(b) and 26(b), which clearly contrast the smooth peak(s) for 
$y_{-\infty } > y^{sep}_{-\infty }$ with the irregular dependence that ensues for 
$y < y^{sep}_{-\infty }$. Note that it is 
$\phi _{j,-\infty }$ that is now fixed at 0, and the value of 
$\phi _{j0}$ emerges during the course of the integration; one finds 
$\phi _{j0} = 111.8^{\circ }$ and 
$168.6^{\circ }$ for the separatrices corresponding to 
$\xi _0=2.0$ and 
$\xi _0=1.05$, respectively. Thus, the separatrices marked by the red vertical lines in figure 25 and 26 correspond to fore–aft asymmetric separatrices with finite upstream offsets – keeping in mind the pathline configurations in figure 7, these separatrices correspond to 
$\phi _{j0} \in ({\rm \pi} /2,{\rm \pi} )$, and are obtained from those in figure 7 by reflecting about the gradient axis. Similar to figure 15(a), the dependence for 
$y_{-\infty }$ values less than 
$y^{sep}_{-\infty }$ is not entirely irregular. Instead, there are intervals of regular dependence separated by chaotic ‘bursts’, with the relative sizes of the regular and chaotic intervals being clearly sensitive to the aspect ratio. Figures 25(b,c) and 26(b,c) correspond to magnified views, with each view corresponding to 
$10^{4}$ initial conditions (
$y_{-\infty }$ values), and highlight the interlacing of regular and chaotic intervals down to the smallest (numerically) resolved scales; in other words, what appears as an interval of chaotic dependence at a given coarse resolution contains smaller intervals of regular intervals within, and this appears to continue ad infinitum. Figures 25(c) and 26(c) suggest qualitatively different transitions from regular to chaotic scattering for the two aspect ratios. The transition for the spheroid with a near-unity aspect ratio (
$\xi _0=2$) in figure 25(c) appears to be discontinuous, with 
$\mathcal {D}t$ approaching a finite limit for 
$y_{-\infty } \rightarrow y^{sep+}_{-\infty }$, but appearing to diverge for 
$y_{-\infty } \rightarrow y^{sep-}_{-\infty }$; in contrast, the transition for the larger aspect ratio spheroid appears to involve a divergence of the residence time regardless of the direction of approach towards the separatrix offset. Despite this difference in appearance, one expects 
$\lim _{y_{-\infty } \rightarrow y^{sep+}_{-\infty }}\mathcal {D}t$ to be finite, since this corresponds to the finite residence times of the fore–aft asymmetric separatrices mentioned above in the interval 
$(-20,20)$ – the inset in figure 26(c) shows this to be the case.
The residence time for a tumbling prolate spheroid (
$\xi _0=2.0(\kappa =1.13$)), with 
$\phi _{j,-\infty }=0$, as a function of the upstream gradient offset of the fluid pathline (‘upstream’ here corresponds to 
$x = -20$); panels (a–c) correspond to increasing levels of magnification. The red vertical line in all figures denotes the separatrix (
$y^{sep}_{-\infty }$) that separates the regular (
$y_{-\infty } > y^{sep}_{-\infty }$) region from the one that includes intervals of chaotic scattering (
$y_{-\infty } < y^{sep}_{-\infty }$).
The residence time for a tumbling prolate spheroid (
$\xi _0=1.05(\kappa =3.28$)), with 
$\phi _{j,-\infty }=0$, as a function of the upstream gradient offset of the fluid pathline (‘upstream’ here corresponds to 
$x = -20$); panels (a–c) correspond to increasing levels of magnification. The red vertical line in all figures denotes separatrix (
$y^{sep}_{-\infty }$) that separates the regular (
$y_{-\infty } > y^{sep}_{-\infty }$) region from the one that includes intervals of chaotic scattering (
$y_{-\infty } < y^{sep}_{-\infty }$).
The separatrix offsets in figure 25 and 26 arise for 
$\phi _{j,-\infty } = 0$. A different offset will result for another choice of 
$\phi _{j,-\infty }$. To access all possible separatrices, one needs to choose 
$\phi _{j,-\infty }$ values spanning the interval (
$0,{\rm \pi}$) for a given large negative 
$x$. The maximum among all such separatrix offsets would then determine the onset of chaotic scattering in the flow–gradient plane. The restriction of a large negative 
$x$ is necessary so that the aforementioned maximum separatrix offset is independent of 
$x$; although, this maximum will occur at an 
$x$-dependent 
$\phi _{j,-\infty }$. Thus, in order to characterize the onset of chaotic scattering in the flow–gradient plane, one ought to have a plot, as depicted in figure 27, of the residence-time surface as a function of both 
$y_{-\infty }$ and 
$\phi _{j,-\infty }$. The onset of chaotic scattering is demarcated by a critical curve 
$y^{sep}_{-\infty }(\phi _{j,-\infty })$ in this plane. The residence-time surface is smooth on one side of this separatrix curve, but one expects a jagged irregular distribution of peaks on the other side; the latter irregular dependence has indeed been observed in the context of the laminar–turbulent transition (Eckhardt et al. Reference Eckhardt, Schneider, Hof and Westerweel2007). Note that, provided 
$x$ is sufficiently large, the set of 
$y^{sep}_{-\infty }$ values is invariant to a change in 
$x$. Instead, and as mentioned above, each 
$y^{sep}_{-\infty }$ now corresponds to a different 
$\phi _{j,-\infty }$, the relation between the two 
$\phi _{j,-\infty }$ values being determined by the angle through which the spheroid rotates as the fluid element translates from the first to the second value of 
$x$ (this translation is almost entirely due to the ambient simple shear). This ‘juggling’ of 
$y^{sep}_{-\infty }$ values will, of course, change the separatrix curve in the 
$y_{-\infty }-\phi _{j,-\infty }$ plane. In both figure 27(a) and 27(b), the separatrix curve has a step discontinuity; as will be seen below, this jump arises due to the abrupt appearance of a chaotic burst on a previously regular peak.
Residence-time surfaces on the 
$y_{-\infty }-\phi _{j,-\infty }$ plane, for tumbling spheroids with (a) 
$\xi _0=1.05(\kappa =3.28)$) and (b) 
$\xi _0=2(\kappa =1.15)$); ‘
$-\infty$’ here corresponds to 
$x =$ - 20. The red separatrix curves mark the transition from regular to chaotic scattering in both cases.
Having characterized the residence-time surface, we return to the specific case of 
$\phi _{j,-\infty } = 0$, and now compare the manner in which the signatures of chaotic scattering disappear as one approaches the limit of an integrable system (note that the term ‘integrable’ is used here in a loose sense, to denote the regular nature of the trajectories, and the resulting smooth dependence on initial conditions; the problem of a spheroid in a simple shear flow is evidently non-Hamiltonian). The two limiting integrable cases are that of a spinning spheroid of any aspect ratio, the streamline topology for which has already been discussed in detail in § 3 (and appendix A), and an infinitely slender prolate spheroid (
$\kappa \rightarrow \infty$ or 
$\xi _0 \rightarrow 1$) in any of the finite-
$C$ Jeffery orbits. We examine the latter limit first. The integrability in this case arises because viscous slender body theory (Batchelor Reference Batchelor1970) shows that the time dependent disturbance velocity field associated with a slender fibre (responsible for the non-integrability of the governing system (2.11)–(2.13)), in a region around it of order its length, scales as the inverse of the logarithm of its (large) aspect ratio. Thus, at leading logarithmic order, a fluid element is merely convected by the ambient simple shear flow. Such an approximation has been earlier used to analyse the effect of polymeric stresses on fibre motion in simple shear flow; in particular, to analyse the viscoelasticity-induced drift across (meridional) Jeffery orbits for large Deborah numbers (Harlen & Koch Reference Harlen and Koch1993). In this limit, the polymeric stress at a point is approximated based on the integrated effect of velocity gradients experienced by a fluid element (containing polymer molecules) along an ambient flow streamline that convects the element up to the point of interest. Figure 28(a–e) show a comparison of the residence-time distributions for tumbling prolate spheroids of increasing aspect ratio, approaching the slender fibre limit; 
$\kappa$ ranges from 2.4 (
$\xi _0$ = 1.1) in figure 28(a) to 70.7 (
$\xi _0$ = 1.0001) in figure 28(e). Note that the length scale characterizing the undulations in the residence-time curves, induced by the spheroid orientation dynamics, increases with increasing 
$\kappa$. As mentioned earlier, this is due to an increase in the Jeffery period. Each of the residence-time curves is overlaid on a smooth curve that plots the residence time estimated based on convection by the ambient shear alone and that is therefore proportional to 
$1/y^{-\infty }$. For the smaller aspect ratios, notwithstanding the chaotic bursts, the ambient-flow-based estimate, ends up overestimating the residence time. This is because accounting for the disturbance velocity field leads to the fluid element being advected to larger 
$y$ in the vicinity of the rotating spheroid, in turn leading to a faster (local) convection by the ambient flow. For the largest aspect ratios (figure 28d,e), the ambient-flow-based estimate is quite accurate, validating to some extent a slender-body-theory-based analysis that neglects of the disturbance velocity field as far as the fluid pathlines is concerned. From figure 28(d) and 28(e), it is also evident that the separatrix offset, 
$y^{sep}_{-\infty }$, that marks the onset of the chaotic bursts, eventually decreases with increasing aspect ratio, again emphasizing the dominance of the ambient flow, and thence, the approach towards an integrable limit 
$(\kappa = \infty )$. This approach is, however, a singular one, with the intermittent chaotic bursts corresponding to large departures of the residence time from the 
$O(1/y^{-\infty })$ ambient-flow-based estimate even for the largest 
$\kappa$ shown. As shown in figure 28(e), these departures correspond to a ‘trapping’ of the fluid pathline for a long time in the vicinity of the rotating spheroid. Such trapped pathlines execute a large number of loops around the spheroid; in contrast, the pathlines above the separatrix, and those in the regular intervals for 
$y_{-\infty } < y^{sep}_{-\infty }$, are open. Interestingly, for the largest 
$\kappa$, the dips of the residence-time curve below the ambient-flow-based estimate correspond to regular reversing pathlines that do not loop around the spheroid (see inset in figure 28e). Finally, note that the signatures of chaotic scattering are also not readily apparent for the largest 
$\kappa$. Thus, although we have marked a separatrix offset in this case, it is possible that the irregular dependence of the residence time, in a manner resembling the smaller 
$\kappa$ values, now occurs on scales smaller than those resolved.
Variation of the residence time for tumbling prolate spheroids with 
$\xi _0=1.1,1.05,1.01,1.001,1.0001(\kappa =2.4,3.28,7.12,22.68,70.7)$) with 
$\phi _{j,-\infty }=0$. The red vertical line separates the regular and chaotic regions, the magenta curve denotes the 
$1/y_{-\infty }$ estimate based on the ambient simple shear flow. The insets in (e) show sample pathlines for upstream offsets less and greater than that of the separatrix (
$y^{sep}_{-\infty }$); those below the separatrix, and within the chaotic burst intervals, loop around the spheroid a large number of times (the two insets on the left), while the other pathlines are open (the insets on the right).
Based on the residence-time distributions shown in figure 28, and those for other aspect ratios, figure 29(a) plots the upstream (gradient) offset of the separatrix (
$y^{sep}_{-\infty }$), marking the onset of chaotic scattering, as a function of the spheroid aspect ratio. Although the separatrix offset appears to eventually decrease monotonically for the largest aspect ratios with 
$y^{sep}_{-\infty } \sim \kappa ^{-4.5}$, for 
$\kappa \rightarrow \infty$, the variation for moderately large aspect ratios is non-monotonic on account of a series of seemingly discontinuous jumps in 
$y^{sep}_{-\infty }$ at certain aspect ratios (similar to the jumps in the separatrix curves seen in figure 27a,b). To examine this further, we plot in figure 29(b) the residence-time distributions for a pair of aspect ratios on either side of a particular 
$y^{sep}_{-\infty }$ jump. It is seen from the distributions that the jump arises from the sudden appearance of a chaotic burst in a previously smooth peak (located just above the separatrix) with increasing aspect ratio. Note that, even for a slender fibre, there would still be fluid elements that approach the fibre to within a distance of order its diameter, and which would then suffer a strong interaction (since the disturbance field is not small in this region). But, the fraction of fluid elements that undergo such a strong interaction must approach zero as 
$\kappa ^{-1}$, owing to the smallness of the ‘collisional’ cross-section which scales as the product of the fibre length and diameter. Thus, the upstream gradient offset of fluid elements that undergo a strong interaction is expected to scale as 
$\kappa ^{-1}$. While 
$y^{sep}_{-\infty }$ in figure 29(a) asymptotes to zero more rapidly, as 
$\kappa ^{-4.5}$, the numerical values of the offsets remain larger than the spheroid diameter for the largest aspect ratios examined. We have verified that the separatrix offsets for oblate spheroids exhibit a similar series of jumps for 
$\kappa$ decreasing to zero (not shown). Unlike the prolate case, the oblate spheroid disturbance field remains finite in the limit 
$\kappa \rightarrow 0$, and one might therefore expect the separatrix offset for a flat disk to also remain finite. Over the range of aspect ratios examined, however, we find separatrix offsets for oblate and prolate spheroids to remain comparable in magnitude.
(a) The upstream gradient offset of the separatrix, 
$y^{sep}_{-\infty }$, for tumbling spheroids, as a function of the spheroid aspect ratio, (b) residence-time plots for the (circled) aspect ratios 
$\kappa =20.18$ and 20.099, in the close vicinity and on either side of a 
$y^{sep}_{-\infty }$ jump.
Next, in figure 30, we present a comparison of the residence-time distributions as a function of the orbit constant 
$C$, for a prolate spheroid of 
$\xi =1.05(\kappa =3.28)$. The sequence of orbit constants chosen is 
$C=1$, 0.5, 0.1, 0.04, 0.01, and thus approaches the other integrable limit (a spinning spheroid: 
$C = 0$) mentioned above. Similar to figure 28, and for purposes of comparison, we plot the residence-time distribution corresponding to a spinning spheroid, of the same aspect ratio, in each of the panels. For a spinning spheroid, and for the finite 
$x$-interval 
$(-20,20)$ under consideration, the residence time remains finite regardless of 
$y_{-\infty }$. Since closed streamlines correspond to an infinite residence time, the spinning-spheroid curve shown terminates at the separatrix (
$y^{sep(spin)}_{-\infty }$); this is unlike figure 15(b) where it was continued to the closed orbits within, by taking the residence time to be equal to the orbital half-periods. Although for the aspect ratio examined in figure 30, there remains a large difference between the offsets of the actual separatrix (
$y^{sep}_{-\infty }$) and the spinning-spheroid separatrix (
$y^{sep(spin)}_{-\infty }$) down to the smallest 
$C$ values, the location of the separatrix for the smallest 
$C$ values nevertheless begins to correlate with a rather abrupt increase in the residence times, in a manner resembling the spinning-spheroid distribution. Unlike the slender fibre limit examined above, figure 31 shows that the separatrix offset starts from an order-unity value (specific to 
$\kappa = 1.05$) for large 
$C$, and decreases monotonically for 
$C$ going to zero, implying a similar decrease of the size of the chaotic region as one approaches the spinning-spheroid limit.
Variation of the residence time, for prolate spheroid with 
$\xi _0=1.05(\kappa =3.28)$ in Jeffrey orbits 
$C=1.0,0.5,0.07,0.04,0.01$ with 
$\phi _{j,-\infty }=0$. The red vertical line in each figure separates the regular and chaotic regions, the magenta curve denotes the residence time for the spinning (
$C=0$) spheroid, with the magenta vertical line denoting the separatrix for the spinning orbit: (a) 
$C=1, t=4000T_j$; (b) 
$C=0.5, t=4000T_j$; (c) 
$C=0.07, t=10\,000T_j$; (d) 
$C=0.04, t=10\,000T_j$; (e) 
$C=0.01, t=20\,000T_j$.
The gradient offset of the separatrix, 
$y^{sep}_{-\infty }$, plotted as a function of the Jeffery orbit constant 
$C$, for a prolate spheroid with 
$\xi _0 = 1.05 (\kappa =3.28)$.
Figure 32 shows the residence-time probability densities for tumbling prolate spheroids of decreasing aspect ratios: 
$\kappa (\xi _0)=22.38(1.001)$, 7.12(1.01), 3.28(1.05), 2.4(1.1), 1.15(2) and 1.03(4). Rather surprisingly, and in contrast to much of the literature on chaotic scattering, the probability densities for all cases examined appear to have algebraic rather than exponential tails (Bleher et al. Reference Bleher, Ott and Grebogi1989, Reference Bleher, Grebogi and Ott1990; Jánosi et al. Reference Jánosi, Tél, Wolf and Gallas1997; Faisst & Eckhardt Reference Faisst and Eckhardt2004; Eckhardt et al. Reference Eckhardt, Schneider, Hof and Westerweel2007). The lower figure on the right shows the invariance of the algebraic tail to the interval of upstream offsets considered. The top figure on the right shows three probability densities corresponding to regular pathlines in (
$y^{sep}_{-\infty },y_{-\infty }$), singular pathlines in (
$0,y^{sep}_{-\infty }$) and all of the pathlines in (
$0,y_{-\infty }$). The tail, corresponding to asymptotically long residence times, is seen to be entirely controlled by the singular pathlines; an increase in 
$y_{-\infty }$ only amounts to considering a greater number of regular pathlines, with shorter residence times, leaving the tails unaffected. An important implication of the algebraic tails in figure 32 is the absence of a characteristic escape time, and the related possibility of a scale dependent structure of the underlying chaotic saddle, in turn implying that the spheroid may be a non-hyperbolic scatterer (Lau, Finn & Ott Reference Lau, Finn and Ott1991). Despite the aforementioned evidence of non-hyperbolicity, however, we find algebraic scaling behaviour for the uncertainty ratio, as seen in figure 24, and as is expected for a hyperbolic scatterer. Further, successively magnified views of the residence-time plots do not reveal an obviously denser distribution of chaotic scattering intervals on the smaller scale. Thus, ascertaining the nature of the scatterer (non-hyperbolic or otherwise) requires a more detailed numerical investigation with an emphasis on scales (the 
$\epsilon$ values in figure 24) much smaller than those examined here.
The main figure shows the residence-time probability densities for tumbling prolate spheroids of different aspect ratios; the slopes of the algebraic tails appear within brackets. The top figure on the right shows the total probability density for 
$\xi _0 = 2$, and the component densities corresponding to regular (
$y^{sep}_{-\infty },y_{-\infty }$) and singular (
$0,y^{sep}_{-\infty }$) pathlines. The bottom figure on the right demonstrates the insensitivity of the tail to the interval of upstream offsets considered.
6. Discussion and conclusions
6.1. Summary of main results
In this study, we have examined the topology of the fluid pathlines induced by freely rotating neutrally buoyant spheroids, both prolate and oblate, and of an arbitrary aspect ratio, in simple shear flow. The fluid pathline configuration is analysed as a function of the particular spheroid orbit, as defined by the orbit constant 
$C$. In the spinning or log-rolling (
$C = 0$) orientations, one has a steady scenario, and the fluid pathlines are the same as streamlines. Further, the topology of the streamline configuration is identical to that already known for the limiting cases of a sphere or a cylinder. There are two distinct groups of streamlines, open and closed, and these are separated by a surface of limiting streamlines, termed separatrices, that ‘closes at infinity’. For spinning prolate spheroids, various measures relating to the streamline configuration, including the size of the closed-streamline region, the orbital periods of the individual closed streamlines, etc, are bounded between those for a sphere and a cylinder (see figure 5).
On the other hand, when the spheroid is in any of the precessional orbits, including the tumbling one (
$0 < C \leqslant \infty$), almost all pathlines are open, in that they come from and eventually go to infinity. Nevertheless, the open pathlines may be divided into two groups. The first group are the regular open pathlines which come from and go to infinity without looping around the precessing spheroid. The residence time for these pathlines, in a certain neighbourhood of the spheroid, is a smooth function of their initial coordinates defined as their offsets, far upstream, in the gradient and vorticity directions. One may therefore regard these regular pathlines as trivial generalizations of the open pathlines around a sphere or a cylinder. In general, these pathlines are fore–aft asymmetric (with respect to the gradient–vorticity plane), and may therefore be characterized via the net displacements in the gradient or vorticity directions (see figure 7). For the most slender spheroids, there is the appearance of regular reversing pathlines that are absent for a torque-free sphere or a cylinder.
The most important results of this manuscript concern the second group of open pathlines, the singular ones, that loop around the precessing spheroid on their way from upstream to downstream infinity. These differ in a profound manner from the closed pathlines around a freely rotating cylinder or a sphere. It is observed that the number of loops, or alternatively, the residence time of these pathlines in a certain neighbourhood of the spheroid, has an extremely sensitive, seemingly random, dependence on the initial offset, with this sensitivity persisting until the smallest (numerically) resolved scales. Such a fractal dependence on initial coordinates is a defining characteristic of chaotic scattering as investigated earlier in other scenarios that include the laminar-to-turbulent transition in plane Couette flow (Skufca et al. Reference Skufca, Yorke and Eckhardt2006), pipe flow (Faisst & Eckhardt Reference Faisst and Eckhardt2004; Eckhardt et al. Reference Eckhardt, Schneider, Hof and Westerweel2007; Schneider, Eckhardt & Yorke Reference Schneider, Eckhardt and Yorke2007), dynamics of point vortices (Aref & Pomphrey Reference Aref and Pomphrey1980; Aref Reference Aref1983; Eckhardt Reference Eckhardt1988), to name a few. The infinitely sensitive variation of the residence time found here implies the existence of a chaotic saddle, and thence, of a set of singularities on a Cantor-like set in the vicinity of the rotating spheroid; for the particular case of simple shear flow examined here (see final paragraph below), this singular set might extend to infinity along the flow axis. We have investigated the nature of the unstable manifold by numerically mimicking the experiment of Gouillart et al. (Reference Gouillart, Dauchot, Thiffeault and Roux2009), whereby the underlying fractal structure is readily evident (figure 22); a partial representation of the invariant chaotic saddle, for a near sphere, is shown in figure 23. We have also investigated in detail the transition from regular to chaotic scattering as a function of the underlying physical parameters that include the upstream gradient offset, the orbit constant and the spheroid aspect ratio. This included looking at the nature of approach to the known integrable limits – that of an arbitrary aspect ratio spheroid in the spinning orbit (figure 31), and that of an infinitely slender prolate spheroid in any orbit (figure 29). An aspect not investigated in detail is the variation in the spatial extent of the chaotic saddle as one moves away from the flow–gradient plane, that is, as a function of the upstream vorticity offset of the pathlines. Figure 33 presents the results of a limited investigation along these lines, and shows that (i) the scenario, for a fixed non-zero 
$z_{-\infty }$, remains similar to the in-plane case, and (ii) that the separatrix offset 
$y^{sep}_{-\infty }$ decreases with increasing 
$z_{-\infty }$.
The gradient offset of the separatrix, 
$y^{sep}_{-\infty }$, plotted as a function of the vorticity offset 
$z_{-\infty }$, for a prolate spheroid with 
$\xi _0 = 1.05 (\kappa =3.28)$, with the sample residence time distributions for 
$z_{-\infty }=$ 0.1 and 1.0
One needs to emphasize the distinction between Hamiltonian and non-Hamiltonian settings here. In the former case, one can relate the onset of chaos to the onset of bound orbit sequences that correspond to infinite yet non-recurring sequences of reflection from the Hamiltonian peaks (as indicated by Eckhardt Reference Eckhardt1988; Eckhardt & Aref Reference Eckhardt and Aref1988). It is also not difficult to identify a critical energy level, based on the governing Hamiltonian, that signals the onset of the aforementioned bound sequences. There does not seem to exist any such analogy in our case, and thus, no obvious way to predict the offset of the separatrices that mark the onset of chaotic scattering. The latter is true for the other non-Hamiltonian problems too (the laminar–turbulent transition, for example). Finally, it is important to mention that we have found the probability densities of residence times, for all cases examined, to have algebraic rather than exponential details, suggesting that the underlying scatterer has a non-hyperbolic character. The non-hyperbolic character has previously been attributed to a no-slip boundary in the vicinity of the chaotic saddle (Gouillart et al. Reference Gouillart, Dauchot, Thiffeault and Roux2009), but is more difficult to rationalize in the present problem, since the spheroid rotates at a finite rate at all times.
6.2. Implications for transport in multiphase systems
As indicated in the Introduction, the findings of the present investigation might have profound implications for transport problems in disperse multiphase systems. In what follows, we briefly discuss three instances, in roughly increasing order of complexity. The first is the transport of heat or mass from particles in shearing flows. As originally shown by Acrivos (Reference Acrivos1971) and Poe & Acrivos (Reference Poe and Acrivos1976), the rate of scalar transport from a freely rotating sphere in a planar linear flow has an exceptional character; planar linear flows form a one-parameter subset of the general family of linear flows (see discussion in the final paragraph below, for the values that this parameter takes), with the limiting values of the parameter corresponding to planar extension and solid-body rotation. The dimensionless rate of transport, as characterized by the Nusselt number (
$Nu$), does not increase indefinitely with the Péclet number (Pe). Such an increase would occur if, at large Pe corresponding to the convectively dominant limit, the transport were to occur via a thin boundary layer, as is the case for a translating sphere; see Acrivos & Goddard (Reference Acrivos and Goddard1965). Instead, Nu saturates at an order-unity value dependent on flow type, the saturation arising due to the transport eventually being limited by diffusion across closed streamlines that surround the sphere in all of the planar linear flows (except planar extension). The singular alteration of the trajectory topology owing to a deviation of the particle from sphericity, found here, must manifest as qualitatively different Nu v/s Pe curves for a sphere and spheroid in a given planar linear flow. In contrast to a sphere, for large enough Pe, one expects the transport from a freely rotating spheroid to be dominated by the set of singular open pathlines which convect heat or mass away in an efficient manner. This should avoid the diffusion limitation for a sphere, leading to a continued growth of Nu with Pe, with the nature of this growth dependent on the residence-time distributions of fluid elements obtained here (for instance, those shown in figure 25); any scaling exponent characterizing the large-Pe asymptotics of Nu, would be a sensitive function of aspect ratio via the residence time distributions. For spheroids with near-unity aspect ratios (the weakly non-integrable cases), in particular, one expects Nu to exhibit an intermediate plateau on account of diffusion limitation, similar to a sphere. For sufficiently large Pe, however, Nu will increase in a manner determined by the fraction of fluid pathlines with residence times shorter than a characteristic diffusion time through the (hypothetical) envelope of closed streamlines. Recent efforts (Subramanian & Koch Reference Subramanian and Koch2005, Reference Subramanian and Koch2006b; Krishnamurthy & Subramanian Reference Krishnamurthy and Subramanian2018) have examined the singular effect of weak inertia-induced convection on the originally diffusion-limited behaviour of Nu in the Stokesian limit. The above discussion suggests that departure from sphericity again have a singular effect, manifesting as a convective enhancement for large Pe. The convective flow in the inertial case may be derived using a straightforward perturbation expansion in the particle Reynolds number, but that for the non-spherical case is non-trivial, and would require the notion of lobe dynamics. The latter has been used earlier to describe the efficient (non-diffusive) exchange of material between the irrotational volume entrained by a vortex ring, or a pair of counter-rotating point vortices, and the ambient (Rom-Kedar, Leonard & Wiggins Reference Rom-Kedar, Leonard and Wiggins1990).
A second problem concerns the rheology of dilute viscoelastic suspensions of anisotropic (axisymmetric) particles which sensitively depends on the orientation dynamics of these particles in an imposed shear. The degenerate nature of the orientation distribution, arising from the existence of closed (Jeffery) orbits in the Stokesian limit (Leal & Hinch Reference Leal and Hinch1971; Dabade et al. Reference Dabade, Marath and Subramanian2015; Marath & Subramanian Reference Marath and Subramanian2017) is well known, and implies that the orientation distribution, in a viscoelastic suspending medium, crucially depends on a viscoelasticity-induced drift across orbits. This drift is governed by the polymeric stress field. For small Deborah numbers (De), the polymeric stress may be determined using an ordered fluid expansion, and is therefore only a function of the local velocity gradient associated with the Stokesian velocity field. But, for large De, the relaxation times are long compared to the characteristic residence time of a fluid element in the vicinity of the rotating particle, and one expects the polymeric stresses to therefore arise from the integrated effects of the velocity gradients seen by a polymer molecule as it is convected along a fluid pathline. For small polymer concentrations, these pathlines will have a near-Newtonian (Stokesian) character, and the non-trivial pathline topology identified here implies that the polymeric stress field, and thence the particle drift, might be crucially influenced by the stresses that develop along the singular open pathlines. In fact, on account of the non-integrable nature of the pathline topology, one expects the spatial variation of polymeric stresses to have a singular character; that is to say, the polymeric stresses at points located arbitrarily close to each other might differ by a finite amount, on account of an analogous difference in the residence times of fluid pathlines arriving at these points. This in turn points to a probabilistic (rather than deterministic) formulation of the particle–polymer–molecule interaction where, rather than determining the polymer stress field at each point in the domain, one attempts to instead determine the probability distribution of polymeric stresses that develops in the neighbourhood of a single point (within the domain influenced by the chaotic saddle). For large aspect ratios, the above problem admits a well known simplification, on account of the disturbance velocity field being logarithmically small (as discussed in § 5). For finite De, this simplification allows one to approximate the polymeric stresses by integrating a functional of the Stokesian velocity gradients along an ambient streamline instead of the actual fluid pathline. This simplification has been exploited in earlier calculations (for instance, see Harlen & Koch Reference Harlen and Koch1993). While such a simplification is certainly valid for steady velocity fields, as for instance, the flow around a sedimenting spheroid with an unchanging orientation in the Stokesian limit (see § 5 in Dabade et al. (Reference Dabade, Marath and Subramanian2015) and Li, Thomases & Guy (Reference Li, Thomases and Guy2019)), the scenario in an ambient shear flow might be more complicated even in the slender fibre limit. As shown in figure 28, the approach to this limit has a singular character – while the residence time, based on the convection by the ambient shear alone, does give an accurate estimate for the regular intervals, for sufficiently large aspect ratios, the occasional chaotic bursts nevertheless lead to very pronounced departures of the residence time from the ambient-shear-based envelope. It would be of interest to examine quantitatively the consequences of such departures for the polymeric stresses that develop along the heavily looped pathlines.
The pathline topology found here is also relevant to inertio-elastic focussing of anisotropic particles in microfluidic channels with applications towards shape sorting (Yang et al. Reference Yang, Kim, Lee, Lee and Kim2011; Amini, Lee & Di Carlo Reference Amini, Lee and Di Carlo2014). The focusing arises due to transverse migration driven by elasto-inertial lift forces. For small Reynolds number (
$Re$) and weakly elastic fluids, the contributions from the elastic and inertial nonlinearities will be additive. The inertial lift force arises from the flow on length scales of order the channel cross-sectional dimension, while the dominant contributions to the elastic lift force arise from smaller length scales of order the particle size. Scaling arguments based on this may be used to show that the inertial lift force is 
$O(Re\,\mu U a^{3}/H^{2})$, while the elastic lift force must be 
$O(De \,\mu U a^{4}/H^{3})$ (see Ho & Leal Reference Ho and Leal1974, Reference Ho and Leal1976). As a result, a non-trivial inertio-elastic equilibrium arising from a balance of the two lift force contributions requires 
$De$ to be 
$O(Re\,H/a)$ where 
$H/a \gg 1$ is the ratio of channel to particle size. Thus, for 
$Re \sim O(1)$ or greater, calculating the elastic lift force will involve determining polymeric stresses, for large 
$De$, along fluid pathlines along the lines mentioned above, and this calculation will beset by the same issues.
The third problem which, in fact, motivated the examination of the fluid pathline topology in the first place concerns pair-spheroid interactions. As mentioned in the introduction, one expects the fluid pathline topology around a single particle to have some bearing on the nature of pair-particle interactions. This is certainly true for a spherical particle, where the streamline and pair-particle-pathline topologies are identical. Thus, the findings here imply that pair-spheroid trajectories are likely to have a similar singular character. As indicated in Marath, Dwivedi & Subramanian (Reference Marath, Dwivedi and Subramanian2017), one would ideally want to characterize pair-spheroid interactions via a phase-averaged scattering kernel that would relate the pre- and post-interaction orbit constants. The findings here would, however, imply that even an infinitesimal change in the input orbit constant would likely lead to an order-unity change in the post-interaction orbit constant. In other words, the scattering kernel would be singular (in fact, such singular kernels in the context of reacting molecules constituted the very first pieces of evidence of chaotic scattering, see Gottdiener Reference Gottdiener1975; Noid, Gray & Rice Reference Noid, Gray and Rice1986). Therefore, similar to the polymer-spheroid problem above, instead of the exact pair-interaction kernel, one would again want to describe even individual pair interactions in a probabilistic framework.
Although a canonical flow-type from the fluid mechanical viewpoint, owing to its significance to rheology, for instance, the choice of simple shear flow does appear inconvenient from the dynamical systems perspective. As is well known (for instance, see Leal & Hinch Reference Leal and Hinch1971; Poe & Acrivos Reference Poe and Acrivos1976), simple shear flow may be regarded as a member of the aforementioned one-parameter family of planar linear flows. The parameter here is 
$\alpha$ (say), with 
$(1+\alpha )/(1-\alpha )$ being proportional to the ratio of extension to vorticity. Simple shear, corresponding to 
$\alpha = 0$, is then the threshold flow, with equal magnitudes of extension and vorticity (leading to straight streamlines), separating the hyperbolic planar linear flows with open streamlines (
$\alpha > 0$), from the elliptic ones with closed streamlines (
$\alpha < 0$). If one now considers a spinning spheroid in a hyperbolic linear flow, from the results of § 3, one expects, similar to a sphere, a trajectory topology where separatrices demarcate the closed from the open streamlines; for instance, see figure 34(a) above (figure 6 in Subramanian & Koch (Reference Subramanian and Koch2006c) which sketches the trajectory topology for a sphere). Importantly, for a hyperbolic linear flow, the closed streamline region, projected onto the flow–gradient plane, has a finite spatial extent for any spheroid aspect ratio, and the aforementioned separatrices are now heteroclinic trajectories connecting a pair of saddle points (see figure 34a). If the analogy with a sphere is exact, then the closed-streamline envelope has a finite volume, and the saddle points will correspond to the points of intersection, of a fixed-point curve, with the flow–gradient plane. All of these fixed points have a saddle-like character, with the curve reducing to a fixed-point circle for the case of a sphere. The existence of such heteroclinic connections highlights the structurally unstable trajectory topology for a spinning spheroid. One may now interpret the deviation of the spheroid, from the spinning orbit to a precessional one corresponding to a small but finite 
$C$, to constitute a perturbation that renders the system non-integrable. Importantly, one expects the non-integrability of the perturbed system to manifest as transverse intersections of the stable and unstable manifolds associated with the two saddle points (see figure 34b), and the existence of such intersections can be rigorously demonstrated via the construction of a Melnikov function (Wiggins Reference Wiggins2003), as has been done earlier in the fluid mechanical context (Kawakami & Funakoshi Reference Kawakami and Funakoshi1999; Angilella, Vilela & Motter Reference Angilella, Vilela and Motter2014). As sketched in figure 34, the implication of these intersections is that intervals of chaotic scattering must replace the original separatrix for the spheroid rotating in any precessional orbit. The trajectory topology in figure 34, although less symmetric, nevertheless resembles the pathline topology analysed by Rom-Kedar et al. (Reference Rom-Kedar, Leonard and Wiggins1990), wherein the non-integrability arose from the time dependence imposed by an oscillatory planar extension acting on a counter-rotating vortex pair. Although related, the emphasis in this work was not on the chaotic scattering aspects, but instead on the exchange of fluid across the (infinitely) convoluted intersecting manifolds. Figure 34 suggests that, from a dynamical systems viewpoint, the non-integrability for the simple shear flow examined here is perhaps best interpreted as the singular limiting case, for 
$\alpha \rightarrow 0^{+}$, of a hyperbolic linear flow, when the saddle points recede to infinity. The nature of chaotic scattering in hyperbolic linear flows, and its relation to the chaotic scattering in the simple shear flow configuration examined here, will be taken up in a future communication.
Smooth and chaotic separatrix in the hyperbolic flow over a sphere (a) and a spheroid (b).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Separatrix of a spinning spheroid
The analytical form of the separatrix for a sphere is known and is given by 
$y=\sqrt {1/3}x^{-3/2}$ (Batchelor & Green Reference Batchelor and Green1972a); an implication of the rather slow approach of the separatrix towards the flow axis is that an infinite volume of fluid is contained within the axisymmetric separatrix envelope. To derive the far-field separatrix for an arbitrary aspect ratio spheroid, we start from the multipole expansion for the velocity field induced by a spheroid in an ambient linear flow, given by (Brenner Reference Brenner1964a,Reference Brennerb, Reference Brenner1966):
\begin{align} \boldsymbol{u}(\boldsymbol{x})&=\boldsymbol{\varGamma}\boldsymbol{\cdot}\boldsymbol{x}+ {\boldsymbol{\mathsf{S}}}\boldsymbol{\cdot}\boldsymbol{\nabla} \frac{3}{D}\frac{\partial}{\partial D}\left[\frac{\sinh D}{D}\right] \frac{\boldsymbol{\mathcal{G}}(\boldsymbol{x})}{8{\rm \pi}\mu} \nonumber\\ &=\boldsymbol{\varGamma}\boldsymbol{\cdot}\boldsymbol{x}+ {\boldsymbol{\mathsf{S}}}\boldsymbol{\cdot}\boldsymbol{\nabla} \left[1+\frac{3.4}{5!}D^{2}+\frac{3.6}{7!}D^{4}+\cdots\right]\frac{\boldsymbol{\mathcal{G}}(\boldsymbol{x})}{8{\rm \pi}\mu}, \end{align}
where 
$D^{2}=d^{2}\overline {\xi _0}^{2}[\nabla ^{2}+(\kappa ^{2}-1)\boldsymbol {pp}\boldsymbol {\colon }\boldsymbol {\nabla }\boldsymbol {\nabla } ](D^{2}=d^{2}\xi _0^{2}[\nabla ^{2}+(\kappa ^{2}-1)\boldsymbol {pp}\boldsymbol {\colon } \boldsymbol {\nabla }\boldsymbol {\nabla }])$ for a prolate (an oblate) spheroid. In (A1), 
$\boldsymbol {p}=\boldsymbol {1}_z$, 
$\boldsymbol {\mathcal {G}}$ is the Oseen–Burger's tensor given by 
$\mathcal {G}_{ij}=\delta _{ij}/r+x_ix_j/r^{3}$, and 
${\boldsymbol{\mathsf{S}}}$ is the stresslet induced by a force and torque-free spinning spheroid in an ambient shear flow, which may be written as:
\begin{equation} {\boldsymbol{\mathsf{S}}}=\frac{20{\rm \pi}\mu d^{3}\xi_0^{3}}{3}K_{tp/o} \left[(\delta_{il}-p_ip_l)E_{lk}(\delta_{kj}-p_kp_j)+\frac{1}{2}(E_{kl}p_kp_l)(\delta_{ij}-p_ip_j)\right]. \end{equation}
In (A2), the stresslet coefficient 
$K_{tp/o}$ is a different function of aspect ratio for prolate and oblate spheroids (see table 2 in the main text), and characterizes the magnitude of the stresslet response for a transverse planar extensional flow.
The dipole (
$\boldsymbol {u}^{d}$) and octupole (
$\boldsymbol {u}^{o}$) contributions correspond to the first and second terms in (A1) and, for a prolate spheroid, are given by
\begin{align} u^{d}_i&=\frac{1}{8{\rm \pi}\mu}S_{jk}\partial_k\mathcal{G}_{ij}={-}\frac{3}{8{\rm \pi}\mu}\frac{S_{jk}n_jn_k}{r^{2}}n_i, \end{align}
\begin{align} u^{o}_i&= \frac{d^{2}\bar{\xi}_0^{2}}{80{\rm \pi}\mu r^{4}}[{-}12S_{ij}n_j+30n_i(S_{jk}n_jn_k)+(\kappa^{2}-1) \{{-}3(4p_i(S_{jk}n_jp_k)+2n_i(S_{jk}p_kp_k))\nonumber\\ &\quad +15(4n_i(p_ln_lS_{jk}n_jp_k)+2p_i(p_ln_lS_{jk}n_jn_k)+n_i(S_{jk}n_jn_k))\nonumber\\ &\quad -105n_i(p_ln_l)^{2}(S_{jk}n_jn_k)\}], \end{align}
where the unit radial vector 
$\boldsymbol {n}=\boldsymbol {x}/r$.
To derive the asymptotic form of the separatrix, one needs to account for the higher-order octupole contribution, owing to the dipole field being purely radial, and thereby, lacking a transverse (
$y$) component along the flow axis. Using the above velocity field in 
$\textrm {d}\boldsymbol {x}/\textrm {d}t = \boldsymbol {\varGamma }\boldsymbol {\cdot }\boldsymbol {x} + \boldsymbol {u}^{d} + \boldsymbol {u}^{o}$, and considering the limit 
$y \rightarrow 0$, 
$x \rightarrow \infty$, one obtains the following equations at leading order:
\begin{gather} \frac{\textrm{d} x}{\textrm{d}t} = \lim_{y\to 0}u_x=y, \end{gather}
\begin{gather}\frac{\textrm{d} y}{\textrm{d}t} =\lim_{y\to 0}u_y={-}\frac{5y^{2}d^{3}\xi_0^{3}\,K_p}{2x^{4}}-\frac{d^{5}\xi_0^{3}\overline{\xi_0}^{2}}{2x^{4}}K_{p}. \end{gather}
It is evident that the 
$x$-component of the velocity in (A5) is just the ambient flow (
$y$) at leading order, while the dominant contribution in (A6) is the octupolar contribution; that this term of 
$O(1/x^{4})$ is dominant over the dipole contribution of 
$O(y^{2}/x^{4})$ is easily seen since 
$y$ for the separatrix is asymptotically small in the far field. With the far-field forms in place, the form of the separatrix may be obtained by taking the ratio of (A5) and (A6), and integrating the resulting differential equation. One obtains 
$y = 1/\sqrt {3}\,d^{5/2}\xi _0^{3/2}\overline {\xi _0}K_p^{1/2}x^{-3/2}$ which exhibits the same far-field decay as the sphere separatrix. For sufficiently large aspect ratios (
$\kappa \rightarrow \infty$ or 
$\xi _0 \rightarrow 1$), the asymptotic form is readily obtained using scaling arguments. 
$u_y \sim O(1/x^{3})$, the octupole contribution in two dimensions, with 
$u_x$ again of 
$O(y)$, so that 
${\textrm {d} y}/{\textrm {d} x} \sim 1/(yx^{3})$, which gives 
$y \sim x^{-1}$.
Appendix B. Separatrix of a tumbling spheroid
The analysis of the separatrix of a tumbling spheroid again requires the expression for the far-field disturbance velocity field due to a spheroid rotating in a precessional orbit. In contrast to a spinning spheroid discussed in appendix A, a torque-free spheroid rotating in a generic Jeffery orbit appears, at leading order, as a time-dependent stresslet. The stresslet is now the sum of three distinct contributions, and may be written in the form
\begin{equation} {\boldsymbol{\mathsf{S}}}=\frac{20{\rm \pi}\mu d^{3}\xi_0^{3}}{3}(K_{ap/o}\,{\boldsymbol{\mathsf{E}}}_{{\boldsymbol{\mathsf{a}}}} +K_{tp/o}\,{\boldsymbol{\mathsf{E}}}_{{\boldsymbol{\mathsf{t}}}}+K_{lp/o}\, {\boldsymbol{\mathsf{E}}}_{{\boldsymbol{\mathsf{l}}}}), \end{equation}
where the subscripts 
$p$ and 
$o$ denote prolate and oblate spheroids, respectively. Here, 
${\boldsymbol{\mathsf{E}}}_{{\boldsymbol{\mathsf{a}}}}$, 
${\boldsymbol{\mathsf{E}}}_{{\boldsymbol{\mathsf{t}}}}$ and 
${\boldsymbol{\mathsf{E}}}_ {\boldsymbol{\mathsf{l}}}$ are the rate of strain tensors associated with axisymmetric extension, transverse and longitudinal planar extensions, respectively, that, in invariant form, are given by
\begin{gather} E_{a\,ij}=(E_{kl}p_kp_l)\left[p_ip_j-\frac{\delta_{ij}}{3}\right], \end{gather}
\begin{gather}E_{t\,ij}=(\delta_{il}-p_ip_l)E_{lk}(\delta_{kj}-p_kp_j)+\tfrac{1}{2}(E_{kl}p_kp_l)(\delta_{ij}-p_ip_j), \end{gather}
\begin{gather}E_{l\,ij}=p_ip_lE_{lk}(\delta_{kj}-p_kp_j)+(\delta_{il}-p_ip_l)E_{kl}p_kp_j. \end{gather}
There are only three distinct stresslet coefficients, despite the five independent canonical velocity fields in table 1, since the responses of the spheroid are identical for the pair of planar and transverse longitudinal extensions (owing to its circular cross-section; see Marath & Subramanian Reference Marath and Subramanian2017). The stresslet coefficients 
$K_{ap/o}$, 
$K_{tp/o}$ and 
$K_{lp/o}$ are functions of aspect ratio, and the relevant expressions for both prolate and oblate spheroids appear in table 2. Note that the time dependence of the stresslet in (B1) arises from the spheroid orientation vector 
$(\boldsymbol {p})$ being time dependent; 
$\boldsymbol {p} = \sin {\theta _j}\cos {\phi _j}\,\boldsymbol {1}_x+\sin {\theta _j}\sin {\phi _j}\, \boldsymbol {1}_y+\cos {\theta _j}\,\boldsymbol {1}_z$, where the time evolution of the Jeffery angles is given by (2.12) and (2.13). The lateral displacements may now be determined, using (B1), as
\begin{equation} \Delta y=\int_{-\infty}^{\infty}\boldsymbol{1}_y\boldsymbol{\cdot}\boldsymbol{u}^{d}\,\textrm{d}t, \quad \Delta z = \int_{-\infty}^{\infty}\boldsymbol{1}_z\boldsymbol{\cdot}\boldsymbol{u}^{d}\,\textrm{d}t, \end{equation}
where 
$\boldsymbol {u}^{d} = -(3/(8{\rm \pi} \mu )){\boldsymbol{\mathsf{S}}}\boldsymbol {:}\boldsymbol {nnn}/r^{2}$ is the leading-order dipolar term in the multipole expansion (see (A3)), with 
${\boldsymbol{\mathsf{S}}}$ being defined by (B1). The distance 
$r$, between the fluid element and the spheroid centre, is a function of time, and may be approximated as 
$r = [(y_{-\infty } t)^{2} + y_{-\infty }^{2} (1 + \tan ^{2}\varTheta )]^{1/2}$, where the gradient and vorticity coordinates have been replaced by their upstream values, and we have used that 
$z_{-\infty } = y_{-\infty } \tan {\varTheta }$(see figure 9b). Thus, 
$r$ changes at leading order only due to convection by the ambient shear, with the instant 
$t = 0$ corresponding to the fluid element being in the gradient–vorticity plane, with 
$\phi _j = \phi _{j0}$ (see (2.12) and (2.13)). The unit radial vector is given by 
$\boldsymbol {n}=(t\,\boldsymbol {1}_x+\boldsymbol {1}_y+\tan {\varTheta }\,\boldsymbol {1}_z)/( 1+t^{2}+\tan ^{2}{\varTheta })^{1/2}$ and its projections along the gradient and vorticity directions, appearing in (B5a,b), may be written as 
$n_y=\boldsymbol {n}\boldsymbol {\cdot }\boldsymbol {1}_y=1/(1+t^{2}+\tan ^{2}{\varTheta })^{1/2}$ and 
$n_z = n_y\tan \varTheta$, respectively.
