1. Introduction
Active particles, which include both biological organisms and synthetic particles, have the capability to convert stored energy to directed movement (Schweitzer Reference Schweitzer2007). A large number of active particles can form a dynamic system commonly referred to as active matter. The active constituents in active matter can span a wide range of scales, from nanorobots and microswimmers to larger organisms like birds, fish and even humans (Vicsek & Zafeiris Reference Vicsek and Zafeiris2012; Marchetti et al. Reference Marchetti, Joanny, Ramaswamy, Liverpool, Prost, Rao and Simha2013). In this study, we focus on micron-sized active particles. The widespread existence of microorganisms in natural settings, combined with substantial advancements in microfluidic experimental techniques, has led to an explosion of research focusing on the motion of small active particles, both biological and synthetic, in viscous fluids (Lauga & Powers Reference Lauga and Powers2009; Elgeti, Winkler & Gompper Reference Elgeti, Winkler and Gompper2015; Bechinger et al. Reference Bechinger, Di Leonardo, Löwen, Reichhardt, Volpe and Volpe2016).
Active particles often exist within gradients of a variety of physical quantities, such as heat, light (Jékely Reference Jékely2009) or chemicals (Moran & Posner Reference Moran and Posner2017), and often respond by reorienting themselves to swim up or down these gradients, a behaviour known as taxis. For instance, E. coli is found to display chemotaxis in gradients of oxygen, galactose, glucose, aspartic acid, threonine or serine (Adler Reference Adler1966). Meanwhile, the photophobic behaviour of E. coli can be used to ‘paint’ with bacteria by selective exposure to light (Arlt et al. Reference Arlt, Martinez, Dawson, Pilizota and Poon2018). Here, we focus on taxis due to environments that are mechanically inhomogeneous, specifically where the viscosity is spatially varying. Viscosity gradients can be found in nature when properties of the fluid such as temperature, salinity or even suspended substances are spatially varying. As an example, numerous coral species secrete mucus that builds up on the sea's surface, leading to areas with differing viscosities where marine microorganisms navigate (Wild et al. Reference Wild, Huettel, Klueter, Kremb, Rasheed and Jørgensen2004; Guadayol et al. Reference Guadayol, Mendonca, Segura-Noguera, Wright, Tassieri and Humphries2021). It has also been shown that the movement and distribution of intestinal bacteria is influenced by viscosity variations in the mucus layer (Swidsinski et al. Reference Swidsinski, Sydora, Doerffel, Loening-Baucke, Vaneechoutte, Lupicki, Scholze, Lochs and Dieleman2007).
Previous experimental studies have observed that several microorganisms demonstrate apparent viscotaxis. For example, Leptospira and Spiroplasma are observed to propel up viscosity gradients (Kaiser & Doetsch Reference Kaiser and Doetsch1975; Petrino & Doetsch Reference Petrino and Doetsch1978; Daniels, Longland & Gilbart Reference Daniels, Longland and Gilbart1980; Takabe et al. Reference Takabe, Tahara, Islam, Affroze, Kudo and Nakamura2017). In contrast, E. coli have been observed to swim down the viscosity gradients (Sherman, Timkina & Glagolev Reference Sherman, Timkina and Glagolev1982). Chlamydomonas reinhardtii, a type of green microalgae, demonstrates complex behaviour in viscosity gradients: it accumulates in high-viscosity regions when gradients are weak, but reorients towards low-viscosity regions in strong gradients (Stehnach et al. Reference Stehnach, Waisbord, Walkama and Guasto2021). When interacting with sharp viscosity gradients, this same alga displays dynamics analogous to the refraction of light, as observed experimentally (Coppola & Kantsler Reference Coppola and Kantsler2021) and modelled theoretically (Gong, Shaik & Elfring Reference Gong, Shaik and Elfring2023). Experiments on the effects of viscosity differences on synthetic swimmers are as of yet limited to helical swimmers crossing perpendicular to a viscosity interface and thus not displaying reorientation (Esparza López et al. Reference Esparza López, Gonzalez-Gutierrez, Solorio-Ordaz, Lauga and Zenit2021).
Recently, it was demonstrated by Liebchen et al. (Reference Liebchen, Monderkamp, ten Hagen and Löwen2018) that a purely hydrodynamic mechanism can lead to viscotaxis. In that work, active particles were modelled as interconnected spheres propelled by a fixed thrust in weak viscosity gradients. These particles were shown to display positive viscotaxis due to an imbalance in viscous drag acting on different spheres. Later work included the effect of viscosity variations on thrust using the spherical squirmer model where the particle activity responsible for generating thrust is represented as a surface slip velocity (Lighthill Reference Lighthill1952; Blake Reference Blake1971). It was shown that hydrodynamic interactions between the active slip conditions on the squirmer's surface and the fluid with spatially varying viscosity generally lead to negative viscotaxis (Datt & Elfring Reference Datt and Elfring2019; Shaik & Elfring Reference Shaik and Elfring2021; Gong et al. Reference Gong, Shaik and Elfring2023). The dynamics of a spherical squirmer in spatially varying viscosity that results from non-uniform distribution of nutrients has also been explored (Shoele & Eastham Reference Shoele and Eastham2018). And recently, the scallop theorem (Purcell Reference Purcell1977) was shown to hold in viscosity gradients (Esparza López & Lauga Reference Esparza López and Lauga2023).
While previous work has focused on spherical squirmers, the influence of particle shape on viscotaxis has yet to be investigated. Previous studies using a two-dimensional swimming sheet have shown that speed increases when it moves either along or against gradients (Dandekar & Ardekani Reference Dandekar and Ardekani2020). More recently, it was demonstrated that viscosity gradients can introduce new forces on slender bodies, offering potential ways to control their orientation and drift (Kamal & Lauga Reference Kamal and Lauga2023). Sedimenting spheroids were also shown to reorient in viscosity gradients, unlike in homogeneous fluids (Anand & Narsimhan Reference Anand and Narsimhan2024).
In order to understand the impact of particle shape on swimming in viscosity gradients, in this paper we use a prolate spheroid squirmer as a model microswimmer. Spheroidal squirmers can be used to represent ciliates with non-spherical bodies (such as Tetrahymena thermophila and Paramecium). The model was first proposed by Keller & Wu (Reference Keller and Wu1977), who showed that the streamlines predicted by their model aligned closely with experimental streak photographs of freely swimming and inertly sedimenting Paramecium caudatum. Later, other researchers modified the model by adding a force-dipole mode to represent various types of swimmers, such as pushers or pullers, to examine the behaviour of a single or pair of spheroidal squirmers moving in a narrow slit (Theers et al. Reference Theers, Westphal, Gompper and Winkler2016). More recent work explored the dynamics, power dissipation and swimming efficiency of a spheroidal squirmer in shear-thinning fluids (van Gogh et al. Reference van Gogh, Demir, Palaniappan and Pak2022) using the reciprocal theorem, an approach similar to that which we employ in this work.
We organize this paper as follows. In § 2, we provide the essential mathematical details of an active prolate spheroid swimming in constant viscosity gradients. We then use the reciprocal theorem and asymptotic analysis to derive expressions for the translational and rotational velocity of the particles in § 3. In § 4, we give an analytical expression for the mobility tensor of passive particles subject to an external force and/or torque. In § 5, we calculate the swimming dynamics of active prolate spheroids, and compare our results with those of a spherical squirmer. In § 6, we discuss the effect of disturbance viscosity, and § 7 concludes the paper.
2. Prolate spheroids in viscosity gradients
We consider a prolate spheroid particle in an otherwise quiescent Newtonian fluid. A prolate spheroid has two equatorial semi-axes of equal length, and one polar longer semi-axis (see figure 1 for a schematic). We label the semi-major axis length 
$a$, and the semi-minor axis length 
$b$ (
$b \leq a$). The eccentricity 
$e = \sqrt {1 - (b/a)^2}$ is a measure of the slenderness of the particle, 
$e=0$ being spherical, while 
$e=1$ is infinitely slender. The orientation of the prolate spheroid is defined as the direction 
$\boldsymbol {p}$ along its major axis.
Figure 1. Sketch of a prolate spheroidal active particle swimming in a constant viscosity gradient. Here, 
$a$ and 
$b$ are the lengths of the semi-major and semi-minor axes. The background colour variations depict the viscosity variations.
The viscosity of the fluid 
$\eta (\boldsymbol {x})$ is taken to be non-uniform due to spatial differences in some physical property of the fluid, such as temperature or salinity. Here, we assume a constant viscosity gradient
\begin{equation} \boldsymbol{\nabla} \eta = \frac{\eta_{\infty}}{L}\,\boldsymbol{d}, \end{equation}
where 
$\eta _\infty /L$ is the magnitude and 
$\boldsymbol {d}$ the direction of the viscosity gradient. The size of the particle is assumed to be small compared with the macroscopic length scale of the variation of viscosity, 
$L$, so we introduce a small parameter 
$\varepsilon = a/L \ll 1$. The viscosity gradient can then be written as 
$\boldsymbol{\nabla} \eta = \varepsilon ({\eta _{\infty }}/{a}) \boldsymbol {d}$.
The fluid surrounding the particle is assumed to be incompressible and Newtonian. In the limit of zero Reynolds number, the governing equations for the flow induced by the particle are
$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\sigma} = \boldsymbol{0}, \end{gather}$$
where 
$\boldsymbol {u}$ is the velocity field, and 
$\boldsymbol {\sigma }$ is the stress tensor. The stress tensor can be written in the form
\begin{align} \boldsymbol{\sigma} ={-} p \boldsymbol{I} + \eta_{\infty} \dot{\boldsymbol{\gamma}} + \boldsymbol{\tau}_{NN}, \end{align}
\begin{align} \boldsymbol{\tau}_{NN} = (\eta (\boldsymbol{x}) - \eta_{\infty}) \dot{\boldsymbol{\gamma}}, \end{align}
where 
$p$ is the pressure, 
$\dot {\boldsymbol {\gamma }} = \boldsymbol{\nabla} \boldsymbol {u} + (\boldsymbol{\nabla} \boldsymbol {u})^{\rm T}$, and 
$\boldsymbol {\tau }_{NN}$ is the extra deviatoric stress due to viscosity differences (from an arbitrary constant viscosity 
$\eta _\infty$).
The boundary conditions on the velocity field 
$\boldsymbol {u}$ are as follows. The disturbance flow caused by the particle should diminish in the far field,
\begin{equation} \boldsymbol{u} \rightarrow \boldsymbol{0} \quad \text{as} \ |\boldsymbol{r}| \rightarrow \infty, \end{equation}
where 
$\boldsymbol {r} = \boldsymbol {x} - \boldsymbol {x}_c$, with 
$\boldsymbol {x}_c$ the centre of the spheroid. And the fluid velocity should satisfy no-slip conditions on the surface of the particle 
$S_p$:
\begin{equation} \boldsymbol{u}(\boldsymbol{x}\in S_p) = \boldsymbol{U} + \boldsymbol{\varOmega} \times \boldsymbol{r} + \boldsymbol{u}^s. \end{equation}
The surface velocity 
$\boldsymbol {u}^s$ arises from activity such as deformation or slip, while the unknown translational and rotational velocities 
$\boldsymbol {U}$ and 
$\boldsymbol {\varOmega }$ are found by enforcing the dynamic conditions on the particle.
We use the prolate spheroidal squirmer model to represent non-spherical active swimmers in this paper. This model is a reasonable representation of ciliates such as Paramecium caudatum that utilize synchronized beating cilia to facilitate movement. The original spheroidal squirmer model developed by Keller & Wu (Reference Keller and Wu1977) includes only one swimming mode, 
$\boldsymbol {u}^s = -B_1 (\boldsymbol {s} \boldsymbol {\cdot } \boldsymbol {p}) \boldsymbol {s}$, where 
$\boldsymbol {s}$ is the unit tangent vector to the surface of the spheroidal microswimmer. Subsequent studies have incorporated the contribution of a force-dipole into this model as a second mode. Following Theers et al. (Reference Theers, Westphal, Gompper and Winkler2016) and van Gogh et al. (Reference van Gogh, Demir, Palaniappan and Pak2022), the slip velocity in our model is expressed as
\begin{equation} \boldsymbol{u}^s ={-}B_1 (\boldsymbol{s} \boldsymbol{\cdot} \boldsymbol{p}) \boldsymbol{s} - B_2 \left(\frac{\boldsymbol{r}}{a} \boldsymbol{\cdot} \boldsymbol{p} \right) (\boldsymbol{s} \boldsymbol{\cdot} \boldsymbol{p}) \boldsymbol{s}. \end{equation}
The sign of squirming ratio 
$\beta = B_2 / B_1$ can be used to divide the swimmers into three types: pushers (
$\beta < 0$), pullers (
$\beta > 0$) and neutral swimmers (
$\beta = 0$). Pushers, such as E. coli, generate propulsion from the back. Chlamydomonas reinhardtii, on the other hand, is categorized as a puller because it uses its flagella to pull fluid from the front. Finally, neutral squirmers produce a flow corresponding to a source dipole. The two-mode spheroidal squirmer model simplifies to the spherical squirmer model in the case of zero eccentricity.
Recent research offers a more general representation of the flow field around a spheroidal squirmer, accounting for an infinite number of squirming modes (Pöhnl, Popescu & Uspal Reference Pöhnl, Popescu and Uspal2020). The swimming speed and stresslet of such a squirmer are influenced by more than just the 
$B_1$ and 
$B_2$ modes. However, these additional modes significantly affect the outcome only when the particle is notably slender (Pöhnl et al. Reference Pöhnl, Popescu and Uspal2020), so the two-mode prolate squirmer model is generally sufficient to depict swimming behaviour (Theers et al. Reference Theers, Westphal, Gompper and Winkler2016; Qi et al. Reference Qi, Annepu, Gompper and Winkler2020; Chi et al. Reference Chi, Gavrikov, Berlyand and Aranson2022; van Gogh et al. Reference van Gogh, Demir, Palaniappan and Pak2022). For simplicity, we use only two modes in our calculations.
Finally, in the absence of inertia, the net force and torque on the particle must be zero, i.e.
\begin{equation} {{\boldsymbol{\mathsf{F}}}}_{ext}+{{\boldsymbol{\mathsf{F}}}} = {{\boldsymbol{\mathsf{0}}}}. \end{equation}
Here, 
${{\boldsymbol{\mathsf{F}}}} = [\boldsymbol {F} \ \boldsymbol {L} ]^{\rm T}$ is a six-dimensional vector including both hydrodynamic force and torque, respectively, i.e.
$$\begin{gather} \boldsymbol{F} = \int_{S_p} \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\sigma} \, {\rm d} S, \end{gather}$$
$$\begin{gather}\boldsymbol{L} = \int_{S_p} \boldsymbol{r} \times (\boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\sigma}) \, {\rm d} S, \end{gather}$$
and 
$\boldsymbol {n}$ is the unit normal vector to the surface of the spheroidal particle, whereas 
${{\boldsymbol{\mathsf{F}}}}_{ext} = [\boldsymbol {F}_{ext}\ \boldsymbol {L}_{ext} ]^{\rm T}$ represents any external forces and torques acting on the particle. Enforcing this dynamic condition sets the particle's translational and rotational velocities.
3. Reciprocal theorem
Rather than solving the velocity field due to the spheroid directly, we instead use the reciprocal theorem to project onto operators from a known auxiliary flow in order to obtain the hydrodynamic force and torque. Following the approach outlined by Elfring (Reference Elfring2017), active particle dynamics in a fluid of arbitrary rheology can be written as
\begin{equation} {{\boldsymbol{\mathsf{U}}}} = \hat{{{\boldsymbol{\mathsf{R}}}}}^{{-}1}_{{{\boldsymbol{\mathsf{FU}}}}} \boldsymbol{\cdot} ({{\boldsymbol{\mathsf{F}}}}_{ext} + {{\boldsymbol{\mathsf{F}}}}_{s} + {{\boldsymbol{\mathsf{F}}}}_{NN}), \end{equation}
where 
${{\boldsymbol{\mathsf{U}}}} = [\boldsymbol {U} \ \boldsymbol {\varOmega } ]^{\rm T}$ is a six-dimensional vector including translational and rotational velocities.
The term
\begin{equation} {{\boldsymbol{\mathsf{F}}}}_{s} = \int_{S_p} \boldsymbol{u}^s \boldsymbol{\cdot} (\boldsymbol{n} \boldsymbol{\cdot} \hat{{{\boldsymbol{\mathsf{T}}}}}_{{{\boldsymbol{\mathsf{U}}}}}) \, {\rm d} S \end{equation}
represents the propulsive force and torque exerted by the particle due to the slip velocity 
$\boldsymbol {u}^s$, in a homogeneous Newtonian fluid, while the term
\begin{equation} {{\boldsymbol{\mathsf{F}}}}_{NN} ={-} \int_{\mathcal{V}} \boldsymbol{\tau}_{NN} : \hat{{{\boldsymbol{\mathsf{E}}}}}_{{{\boldsymbol{\mathsf{U}}}}} \, {\rm d} V \end{equation}
accounts for the additional force and torque stemming from the extra deviatoric stress 
$\boldsymbol {\tau }_{NN}$ in the fluid volume 
$\mathcal {V}$ where the squirmer is immersed.
The terms denoted with a hat are linear operators associated with the auxiliary flow solution of rigid-body motion of a body of the same shape in a homogeneous Newtonian fluid of viscosity 
$\eta _{\infty }$. The tensors 
$\hat {{{\boldsymbol{\mathsf{T}}}}}_{ {{\boldsymbol{\mathsf{U}}}}}$ and 
$\hat {{{\boldsymbol{\mathsf{E}}}}}_{{{\boldsymbol{\mathsf{U}}}}}$ are spatially dependent functions that map velocities of the particle 
$\hat {{{\boldsymbol{\mathsf{U}}}}}$ to the stress 
$\hat {\boldsymbol {\sigma }} = \hat {{{\boldsymbol{\mathsf{T}}}}}_{{{\boldsymbol{\mathsf{U}}}}} \boldsymbol {\cdot } \hat {{{\boldsymbol{\mathsf{U}}}}}$ and rate of strain 
$\hat {\dot {\boldsymbol {\gamma }}} = 2 \hat {{{\boldsymbol{\mathsf{E}}}}}_{{{\boldsymbol{\mathsf{U}}}}} \boldsymbol {\cdot } \hat {{{\boldsymbol{\mathsf{U}}}}}$, respectively, while 
$\hat {{{\boldsymbol{\mathsf{R}}}}}_{{{\boldsymbol{\mathsf{FU}}}}}$ is the (
$6\times 6$) resistance tensor. These operators are well known for prolate spheroids (see Appendix A for further details).
The extra stress 
$\boldsymbol {\tau }_{NN}$ due to small viscosity variations is parametrized by 
$\varepsilon$, so we expand all flow quantities in regular perturbation series in 
$\varepsilon$:
\begin{align} \{\boldsymbol{u},
\boldsymbol{\sigma}, \boldsymbol{\tau}_{NN},
\dot{\boldsymbol{\gamma}}, \boldsymbol{U},
\boldsymbol{\varOmega}\} &= \{ \boldsymbol{u}_0,
\boldsymbol{\sigma}_0, \boldsymbol{0},
\dot{\boldsymbol{\gamma}}_0, \boldsymbol{U}_0,
\boldsymbol{\varOmega}_0 \}\nonumber\\
&\quad + \varepsilon \{
\boldsymbol{u}_1, \boldsymbol{\sigma}_1,
\boldsymbol{\tau}_{NN,1}, \dot{\boldsymbol{\gamma}}_1,
\boldsymbol{U}_1, \boldsymbol{\varOmega}_1 \} + O
(\varepsilon^2).
\end{align} At leading order, we have a homogeneous Newtonian fluid of viscosity 
$\eta _{\infty }$. Viscosity variations are captured at the next order, 
$O(\varepsilon )$, where the extra stress is
\begin{equation} \boldsymbol{\tau}_{NN} = (\eta (\boldsymbol{x}) - \eta_{\infty}) \dot{\boldsymbol{\gamma}}_0 + O(\varepsilon^2), \end{equation}
and 
$\dot {\boldsymbol {\gamma }}_0$ is the strain rate of the flow of an active particle in the leading-order homogeneous fluid. Upon substitution of (3.5) in (3.3), we see that calculation of the extra force and torque
\begin{equation} {{\boldsymbol{\mathsf{F}}}}_{NN} ={-} \int_{\mathcal{V}} (\eta (\boldsymbol{x}) - \eta_{\infty}) \dot{\boldsymbol{\gamma}}_0 : \hat{{{\boldsymbol{\mathsf{E}}}}}_{{{\boldsymbol{\mathsf{U}}}}} \, {\rm d} V+ O (\varepsilon^2), \end{equation}
due to spatial variations of viscosity, up to 
$O(\varepsilon )$, requires only the integration of known Stokes flow solutions 
$\dot {\boldsymbol {\gamma }}_0$ and 
$\hat {{{\boldsymbol{\mathsf{E}}}}}_{{{\boldsymbol{\mathsf{U}}}}}$ from the auxiliary resistance problem. Analytical evaluation of the integral is performed most easily in a particle-aligned spheroidal coordinate system, with details given in Appendix B.
We note that the functional form of the viscosity 
$\eta (\boldsymbol {x})$ is restricted insofar as we have assumed a regular perturbation expansion of all flow quantities and 
$\eta =\eta _\infty$ in the limit 
$\varepsilon \rightarrow 0$. When dealing with a linearly varying viscosity field such that 
$(\eta (\boldsymbol {x}) - \eta _{\infty }) \sim \varepsilon x$, the expansion maintains regularity only for 
$x \sim o(1/\varepsilon )$; however, the far-field contribution of a squirmer at distances 
$r \sim O(1/\varepsilon )$ is 
$O(\varepsilon ^2)$ with respect to the non-Newtonian force 
$\boldsymbol {F}_{NN}$, and 
$O(\varepsilon ^3)$ with respect to the non-Newtonian torque 
$\boldsymbol {L}_{NN}$. The velocity field of a passive spheroid decays slower than that of a squirmer, but in constant viscosity gradients, the far-field contribution to the integrals at 
$O(\varepsilon )$ is exactly zero (due to symmetry), making these systems suitable for analysis using a regular perturbation scheme.
4. Passive spheroids
Before examining the dynamics of an active particle, we first derive the mobility of a passive prolate spheroid subject to an external force and torque 
${{\boldsymbol{\mathsf{F}}}}_{ext}$ in a viscosity gradient. For a passive spheroid, there is no active slip, 
$\boldsymbol {u}^s = \boldsymbol {0}$, thus 
${{\boldsymbol{\mathsf{F}}}}_{s}={{\boldsymbol{\mathsf{0}}}}$.
At leading order, 
${{\boldsymbol{\mathsf{F}}}}_{NN}={{\boldsymbol{\mathsf{0}}}}$, and from (3.1), we simply obtain the dynamics of a passive spheroid in a homogeneous Newtonian fluid of viscosity 
$\eta _{\infty }$, under an external force and torque 
${{\boldsymbol{\mathsf{F}}}}_{ext}$, which satisfies the usual mobility relationship (Kim & Karilla Reference Kim and Karilla1991)
\begin{equation} {{\boldsymbol{\mathsf{U}}}}_0 = \hat{{{\boldsymbol{\mathsf{R}}}}}^{{-}1}_{{{\boldsymbol{\mathsf{F}}}} {{\boldsymbol{\mathsf{U}}}}} \boldsymbol{\cdot} {{\boldsymbol{\mathsf{F}}}}_{ext}. \end{equation}
The flow field at this order is identical to the auxiliary flow field in the previous section (see Appendix A), thus the strain rate 
$\dot {\boldsymbol {\gamma }}_0 = 2 \hat {{{\boldsymbol{\mathsf{E}}}}}_{{{\boldsymbol{\mathsf{U}}}}} \boldsymbol {\cdot } {{\boldsymbol{\mathsf{U}}}}_0$ can be written as
\begin{equation} \dot{\boldsymbol{\gamma}}_0 = 2 \hat{{{\boldsymbol{\mathsf{E}}}}}_{{{\boldsymbol{\mathsf{U}}}}} \boldsymbol{\cdot} \hat{{{\boldsymbol{\mathsf{R}}}}}^{{-}1}_{{{\boldsymbol{\mathsf{F}}}} {{\boldsymbol{\mathsf{U}}}}} \boldsymbol{\cdot} {{\boldsymbol{\mathsf{F}}}}_{ext}. \end{equation}At first order, substitution of (4.2) into (3.6) yields
\begin{equation} {{\boldsymbol{\mathsf{F}}}}_{NN} ={-}\boldsymbol{{{\boldsymbol{\mathsf{R}}}}}_{NN} \boldsymbol{\cdot} \hat{{{\boldsymbol{\mathsf{R}}}}}^{{-}1}_{{{\boldsymbol{\mathsf{F}}}} {{\boldsymbol{\mathsf{U}}}}} \boldsymbol{\cdot} {{\boldsymbol{\mathsf{F}}}}_{ext}, \end{equation}where for convenience we have defined the tensor
\begin{equation} \boldsymbol{{{\boldsymbol{\mathsf{R}}}}}_{NN} =\int_{\mathcal{V}}2(\eta (\boldsymbol{x}) -\eta_{\infty})\hat{{{\boldsymbol{\mathsf{E}}}}}_{{{\boldsymbol{\mathsf{U}}}}}:\hat{{{\boldsymbol{\mathsf{E}}}}}_{{{\boldsymbol{\mathsf{U}}}}}\, {\rm d}V. \end{equation}Using (3.1), we obtain the translational and rotational velocity of a passive prolate spheroid at first order:
\begin{equation} \varepsilon{{\boldsymbol{\mathsf{U}}}}_1 ={-}\hat{{{\boldsymbol{\mathsf{R}}}}}^{{-}1}_{{{\boldsymbol{\mathsf{F}}}} {{\boldsymbol{\mathsf{U}}}}} \boldsymbol{\cdot} \boldsymbol{{{\boldsymbol{\mathsf{R}}}}}_{NN} \boldsymbol{\cdot} \hat{{{\boldsymbol{\mathsf{R}}}}}^{{-}1}_{{{\boldsymbol{\mathsf{F}}}} {{\boldsymbol{\mathsf{U}}}}} \boldsymbol{\cdot} {{\boldsymbol{\mathsf{F}}}}_{ext}. \end{equation}
\begin{equation} {{\boldsymbol{\mathsf{U}}}} = {{\boldsymbol{\mathsf{U}}}}_0 + \varepsilon {{\boldsymbol{\mathsf{U}}}}_1 = ( \hat{{{\boldsymbol{\mathsf{R}}}}}^{{-}1}_{{{\boldsymbol{\mathsf{F}}}} {{\boldsymbol{\mathsf{U}}}}} -\hat{{{\boldsymbol{\mathsf{R}}}}}^{{-}1}_{{{\boldsymbol{\mathsf{F}}}} {{\boldsymbol{\mathsf{U}}}}} \boldsymbol{\cdot} \boldsymbol{{{\boldsymbol{\mathsf{R}}}}}_{NN} \boldsymbol{\cdot} \hat{{{\boldsymbol{\mathsf{R}}}}}^{{-}1}_{{{\boldsymbol{\mathsf{F}}}} {{\boldsymbol{\mathsf{U}}}}}) \boldsymbol{\cdot} {{\boldsymbol{\mathsf{F}}}}_{ext}, \end{equation}
we obtain the mobility 
${{\boldsymbol{\mathsf{M}}}}_{{{\boldsymbol{\mathsf{U}}}} {{\boldsymbol{\mathsf{F}}}}} = \hat {{{\boldsymbol{\mathsf{R}}}}}^{-1}_{{{\boldsymbol{\mathsf{F}}}}{{\boldsymbol{\mathsf{U}}}}}-\hat {{{\boldsymbol{\mathsf{R}}}}}^{-1}_{{{\boldsymbol{\mathsf{F}}}} {{\boldsymbol{\mathsf{U}}}}} \boldsymbol {\cdot } \boldsymbol {{{\boldsymbol{\mathsf{R}}}}}_{NN} \boldsymbol {\cdot } \hat {{{\boldsymbol{\mathsf{R}}}}}^{-1}_{{{\boldsymbol{\mathsf{F}}}} {{\boldsymbol{\mathsf{U}}}}}$, connecting the particle velocities 
${{\boldsymbol{\mathsf{U}}}}$ to the external force and torque 
${{\boldsymbol{\mathsf{F}}}}_{ext}$, valid to first order in 
$\varepsilon$, where
\begin{equation} {{\boldsymbol{\mathsf{M}}}}_{{{\boldsymbol{\mathsf{U}}}} {{\boldsymbol{\mathsf{F}}}}} =\begin{pmatrix} {{\boldsymbol{\mathsf{M}}}}_{\boldsymbol{U} \boldsymbol{F}} & {{\boldsymbol{\mathsf{M}}}}_{\boldsymbol{U} \boldsymbol{L}} \\ {{\boldsymbol{\mathsf{M}}}}_{\boldsymbol{\varOmega} \boldsymbol{F}} & {{\boldsymbol{\mathsf{M}}}}_{\boldsymbol{\varOmega} \boldsymbol{L}} \end{pmatrix}, \end{equation}
and 
${{\boldsymbol{\mathsf{M}}}}_{\boldsymbol {U} \boldsymbol {L}} = {{\boldsymbol{\mathsf{M}}}}_{\boldsymbol {\varOmega } \boldsymbol {F}}^{\rm T}$. In homogeneous fluids, the mobility is determined solely by the shape and orientation of the particle, specified by the eccentricity 
$e$ and the orientation vector 
$\boldsymbol {p}$. In viscosity gradients, the mobility also depends on 
$\boldsymbol{\nabla} \eta$. The expressions for the force-translational velocity coupling, 
${{\boldsymbol{\mathsf{M}}}}_{\boldsymbol {U} \boldsymbol {F}}$, and the torque-angular velocity coupling, 
${{\boldsymbol{\mathsf{M}}}}_{\boldsymbol {\varOmega } \boldsymbol {L}}$, are essentially identical to when the viscosity is constant,
$$\begin{gather} {{\boldsymbol{\mathsf{M}}}}_{\boldsymbol{UF}} = \frac{1}{6 {\rm \pi}\eta |_{\boldsymbol{x}_c} a} \left[ \frac{1}{\mathcal{X}^{A}}\,\boldsymbol{p} \boldsymbol{p} + \frac{1}{\mathcal{Y}^{A}}\,({{\boldsymbol{\mathsf{I}}}} - \boldsymbol{p} \boldsymbol{p})\right], \end{gather}$$
$$\begin{gather}{{\boldsymbol{\mathsf{M}}}}_{\boldsymbol{\varOmega} \boldsymbol{L}} = \frac{1}{8 {\rm \pi}\eta |_{\boldsymbol{x}_c} a^3} \left [ \frac{1}{\mathcal{X}^{C}}\,\boldsymbol{p} \boldsymbol{p} + \frac{1}{\mathcal{Y}^{C}}\,({{\boldsymbol{\mathsf{I}}}} - \boldsymbol{p} \boldsymbol{p}) \right ], \end{gather}$$
except that the viscosity is now evaluated at the instantaneous particle centre 
$\boldsymbol {x}_c$. Here, 
$\eta |_{{\boldsymbol {x}}_c} = \eta _{\infty } + {\boldsymbol {x}}_c\boldsymbol {\cdot }\boldsymbol {\nabla } \eta$, hence 
$\eta |_{{\boldsymbol {x}}_c} = \eta _{\infty } + O( \varepsilon )$ as long as the reference viscosity 
$\eta _{\infty }$ is defined close to the particle, 
$| {\boldsymbol {d}} \boldsymbol {\cdot } {\boldsymbol {x}}_c | /a =O(1)$. The coefficients 
$\mathcal {X}^{A}$, 
$\mathcal {Y}^{A}$, 
$\mathcal {X}^{C}$, 
$\mathcal {Y}^{C}$ are functions of eccentricity 
$e$, and their expressions are given in Appendix A.
Unlike in homogeneous Newtonian fluids, in viscosity gradients there arises a torque-translational velocity (and force-angular velocity) coupling
\begin{equation} {{\boldsymbol{\mathsf{M}}}}_{\boldsymbol{U} \boldsymbol{L}} = {{\boldsymbol{\mathsf{M}}}}_{\boldsymbol{\varOmega} \boldsymbol{F}}^{\rm T} = \frac{\varepsilon}{6 {\rm \pi}\eta_{\infty} a^2}\,[ \varLambda_1 ( \boldsymbol{d} \times {{\boldsymbol{\mathsf{I}}}}) + \varLambda_2 (\kern 0.06em \boldsymbol{p} \boldsymbol{\cdot} \boldsymbol{d}) (\kern 0.06em \boldsymbol{p} \times {{\boldsymbol{\mathsf{I}}}}) + \varLambda_3 \boldsymbol{p}(\boldsymbol{d} \times \boldsymbol{p}) ],\end{equation}where
$$\begin{gather} \varLambda_1 = \frac{3[ 2e - (1 - e^2) \mathcal{L}_e]}{16 e^3}, \end{gather}$$
$$\begin{gather}\varLambda_2 = \frac{3[ 2e({-}3 + e^2) + (3 - 2 e^2 + 3e^4) \mathcal{L}_e]}{32 e^3 (2 - e^2)}, \end{gather}$$
$$\begin{gather}\varLambda_3 = \frac{3[ 2e(9 - 5 e^2) + ({-}9 + 8 e^2 + e^4) \mathcal{L}_e]}{32 e^3 (2 - e^2)}, \end{gather}$$
and 
$\mathcal {L}_e = \ln ({(1+e)}/{(1-e)})$.
In the spherical limit (
$e \rightarrow 0$), (4.10) simplifies to
\begin{equation} {{\boldsymbol{\mathsf{M}}}}_{\boldsymbol{U} \boldsymbol{L}} = {{\boldsymbol{\mathsf{M}}}}_{\boldsymbol{\varOmega} \boldsymbol{F}}^{\rm T} = \frac{\varepsilon}{24 {\rm \pi}\eta_{\infty} a^2}\,\boldsymbol{d} \times {{\boldsymbol{\mathsf{I}}}},\end{equation}
and the inverse of the mobility agrees with the resistance tensor for a sphere reported previously by Datt & Elfring (Reference Datt and Elfring2019). Applying a constant external force 
$\boldsymbol {F}_{ext}$ (say due to gravity), the angular velocity due to a viscosity gradient would be 
$\boldsymbol {\varOmega } = (\varepsilon /24{\rm \pi} \eta _\infty a^2) \boldsymbol {F}_{ext}\times \boldsymbol {d}$ or 
$\boldsymbol {\varOmega } = (1/4) \boldsymbol {U}_0\times \boldsymbol {\nabla }(\eta /\eta _\infty )$, where 
$\boldsymbol {U}_0 = \boldsymbol {F}_{ext}/(6{\rm \pi} \eta _\infty a)$.
We have also made a comparison between our calculations and the results for an elongated prolate spheroid sedimenting in viscosity gradients by Kamal & Lauga (Reference Kamal and Lauga2023). At large aspect ratios 
$\lambda = a/b \rightarrow \infty$, the mobilities can be written as
$$\begin{gather} {{\boldsymbol{\mathsf{M}}}}_{\boldsymbol{UF}} \sim \frac{1}{6 {\rm \pi}\eta |_{\boldsymbol{x}_c} a} \left[\frac{3 \ln \lambda}{2}\,\boldsymbol{p} \boldsymbol{p} + \frac{3 \ln \lambda}{4}\,({{\boldsymbol{\mathsf{I}}}} - \boldsymbol{p} \boldsymbol{p}) \right], \end{gather}$$
$$\begin{gather}{{\boldsymbol{\mathsf{M}}}}_{\boldsymbol{UL}} \sim \frac{\varepsilon}{6 {\rm \pi}\eta_{\infty} a^2} \left[ \frac{3}{8}\,( \boldsymbol{d} \times {{\boldsymbol{\mathsf{I}}}}) + \frac{3 \ln \lambda}{4}\,(\kern 0.06em \boldsymbol{p} \boldsymbol{\cdot} \boldsymbol{d}) (\kern 0.06em \boldsymbol{p} \times {{\boldsymbol{\mathsf{I}}}}) + \frac{3}{4}\,\boldsymbol{p}(\boldsymbol{d} \times \boldsymbol{p}) \right], \end{gather}$$
$$\begin{gather}{{\boldsymbol{\mathsf{M}}}}_{\boldsymbol{\varOmega} \boldsymbol{L}} \sim \frac{1}{8 {\rm \pi}\eta |_{\boldsymbol{x}_c} a^3} \left[ \frac{3 \lambda^2}{2}\,\boldsymbol{p} \boldsymbol{p} + 3 \ln \lambda\,({{\boldsymbol{\mathsf{I}}}} - \boldsymbol{p} \boldsymbol{p}) \right]. \end{gather}$$
In this limit, we obtain the mobility matrix for an asymptotically slender spheroid in a constant viscosity gradient. Calculating the leading-order translational and rotational velocities with external force 
$\boldsymbol {F}_{ext} = m \boldsymbol {g}$ and torque 
$\boldsymbol {L}_{ext} = \boldsymbol {0}$, our results coincide exactly with the sedimenting velocities of slender filaments in viscosity gradients found by Kamal & Lauga (Reference Kamal and Lauga2023).
Recent work by Anand & Narsimhan (Reference Anand and Narsimhan2024) also explored the dynamics of sedimenting passive spheroids in viscosity gradients numerically. The authors of that work constructed a dimensionless mobility matrix, and following their approach, we rescale so that the dimensionless torque-translational velocity tensor is
\begin{equation} \tilde{{{\boldsymbol{\mathsf{M}}}}}_{\boldsymbol{U} \boldsymbol{L}} = \frac{6 {\rm \pi}\eta_{\infty} a^2}{\lambda^{4/3}}\,{{\boldsymbol{\mathsf{M}}}}_{\boldsymbol{U} \boldsymbol{L}} = \varepsilon \lambda^{{-}2/3} [ \tilde{\varLambda}_1 ( \boldsymbol{d} \times {{\boldsymbol{\mathsf{I}}}}) + \tilde{\varLambda}_2 (\kern 0.06em \boldsymbol{p} \boldsymbol{\cdot} \boldsymbol{d}) (\kern 0.06em \boldsymbol{p} \times {{\boldsymbol{\mathsf{I}}}}) + \tilde{\varLambda}_3 \boldsymbol{p}(\boldsymbol{d} \times \boldsymbol{p}) ], \end{equation}where
\begin{equation} \tilde{\varLambda}_i = \frac{\varLambda_i}{\lambda^{2/3}}, \quad i = 1,2,3. \end{equation}
We then compare dimensionless coefficients 
$\tilde {\varLambda }_i$ with the corresponding numerical results by Anand & Narsimhan (Reference Anand and Narsimhan2024) for different aspect ratios, as shown in figure 2. We find very good agreement, with only minor discrepancies for 
$\tilde {\varLambda }_2$ and 
$\tilde {\varLambda }_3$ at larger aspect ratios. As another point of comparison, we also calculate the corresponding values of 
$\tilde {\varLambda }_i$ from Kamal & Lauga (Reference Kamal and Lauga2023), and as shown in figure 2, when the aspect ratio is very large, our analytical results align closely.
Figure 2. A plot of mobility coefficients 
$\tilde {\varLambda }_i$ as functions of aspect ratio 
$\lambda$. Solid lines represent the present work, and triangles are those found by Anand & Narsimhan (Reference Anand and Narsimhan2024). Also shown are the data for a sphere from Datt & Elfring (Reference Datt and Elfring2019) (circles) and for an asymptotically slender spheroid from the resistive force theory (RFT) of Kamal & Lauga (Reference Kamal and Lauga2023) (squares).
5. Active spheroids
Microswimmers are often considered to be neutrally buoyant; we do the same here, hence we assume that there is no externally applied force or torque, 
${{\boldsymbol{\mathsf{F}}}}_{ext}={{\boldsymbol{\mathsf{0}}}}$, on an the active spheroid swimming in a viscosity gradient.
At leading order in 
$\varepsilon$, we have an active spheroid swimming in a homogeneous Newtonian fluid of viscosity 
$\eta _{\infty }$. The swim speed is well known (Keller & Wu Reference Keller and Wu1977; Theers et al. Reference Theers, Westphal, Gompper and Winkler2016; Pöhnl et al. Reference Pöhnl, Popescu and Uspal2020; van Gogh et al. Reference van Gogh, Demir, Palaniappan and Pak2022):
\begin{equation} {{\boldsymbol{\mathsf{U}}}}_0 = \hat{{{\boldsymbol{\mathsf{R}}}}}^{{-}1}_{{{\boldsymbol{\mathsf{F}}}} {{\boldsymbol{\mathsf{U}}}}} \boldsymbol{\cdot} {{\boldsymbol{\mathsf{F}}}}_s = \begin{pmatrix} \dfrac{2e - (1 - e^2) \mathcal{L}_e}{2e^3}\,B_1 \boldsymbol{p} \\ \boldsymbol{0} \end{pmatrix}. \end{equation}The corresponding flow field is given in Appendix A.
At first order, the translational and rotational velocities
\begin{equation} \varepsilon{{\boldsymbol{\mathsf{U}}}}_1 = \hat{{{\boldsymbol{\mathsf{R}}}}}^{{-}1}_{{{\boldsymbol{\mathsf{F}}}} {{\boldsymbol{\mathsf{U}}}}} \boldsymbol{\cdot} {{\boldsymbol{\mathsf{F}}}}_{NN} \end{equation}
are obtained using (3.6), with 
$\dot {\boldsymbol {\gamma }}_0$ calculated from the flow field solutions of a two-mode active spheroidal squirmer in Appendix A. Combining (5.1) with (5.2), 
${{\boldsymbol{\mathsf{U}}}} = {{\boldsymbol{\mathsf{U}}}}_0 + \varepsilon {{\boldsymbol{\mathsf{U}}}}_1$, we obtain expressions valid up to 
$O(\varepsilon )$ for the translational and rotational velocities of a prolate spheroidal squirmer:
$$\begin{gather} \boldsymbol{U} = \boldsymbol{U}_0 - \frac{a B_2}{5}\,(\mathcal{X}^{U} \boldsymbol{I} - \mathcal{Y}^{U} 3 \boldsymbol{pp}) \boldsymbol{\cdot} {\boldsymbol{\nabla}} \left(\frac{\eta }{ \eta_{\infty}}\right), \end{gather}$$
$$\begin{gather}\boldsymbol{\varOmega} ={-} \frac{1}{2}\,\mathcal{X}^{\varOmega} \boldsymbol{U}_{0} \times {\boldsymbol{\nabla}} \left(\frac{\eta }{ \eta_{\infty}}\right), \end{gather}$$where the coefficients
$$\begin{gather} \mathcal{X}^U = \frac{5 [{-}6e + 4 e^3 + 3 (1 - e^2) \mathcal{L}_e ] [ - 6e + 10 e^3 + 3 (1 - e^2)^2 \mathcal{L}_e]}{24 e^5 [ 6e - (3 - e^2)\mathcal{L}_e]}, \end{gather}$$
$$\begin{gather}\mathcal{Y}^U = \frac{5 [{-}6e + 4 e^3 + 3 (1 - e^2) \mathcal{L}_e] [ - 18e + 6 e^3 + (9 - 6e^2 + 5e^4) \mathcal{L}_e]}{72 e^5 [ 6e - (3 - e^2) \mathcal{L}_e]}, \end{gather}$$
$$\begin{gather}\mathcal{X}^{\varOmega} = \frac{(1 - e^2) [{-}2 e + (1 + e^2) \mathcal{L}_e ]}{(2 - e^2) [ 2e - (1 - e^2) \mathcal{L}_e ]}, \end{gather}$$
are monotonically decreasing functions of the eccentricity. In the spherical limit 
$e\rightarrow 0$, we have 
$\mathcal {X}^U=1$, 
$\mathcal {Y}^U=1$ and 
$\mathcal {X}^{\varOmega }=1$, and we recover exactly the dynamics for spheres found by Datt & Elfring (Reference Datt and Elfring2019). Conversely, in the slender limit 
$e\rightarrow 1$, we have 
$\mathcal {X}^U=0$, 
$\mathcal {Y}^U=5/9$ and 
$\mathcal {X}^{\varOmega }=0$, meaning that infinitely slender squirmers do not reorient in viscosity gradients, but there is still a change in their translational velocity due to the interaction of the dipolar flow with the spatial variations in viscosity. Generally (for 
$e \lessapprox 0.9988$), the speed change is greater for spheroids than spheres when aligned with the gradient.
In general, the behaviour of spheroidal squirmers is qualitatively similar to spherical squirmers as they navigate through constant viscosity gradients (Datt & Elfring Reference Datt and Elfring2019): all swimmers display negative viscotaxis by reorienting to swim down viscosity gradients, except that the impact of the gradient is diminished with increasing slenderness. The mechanistic reason for this change is straightforward: the viscosity difference across a slimmer body is reduced, which leads to slower reorientation, and in the slender limit viscotaxis ceases. Examining (5.4) more closely, we see that the angular velocity scales as the viscosity difference across the particle, 
${\varOmega }/{(U_0/a)} \sim b\,| \boldsymbol{\nabla} ({\eta }/{\eta _{\infty }})|\,({\ln \lambda }/{\lambda }) = {\varepsilon \ln \lambda }/{\lambda ^2}$. Thus for weak viscosity gradients 
$\varepsilon \ll 1$ and slender particles 
$\lambda \gg 1$, the rotation rate tends to zero. But one can certainly envision a scenario where for sufficiently sharp gradients, 
${\varepsilon \ln \lambda }/{\lambda ^2} = O(1)$. Our formulas are not strictly valid in this setting, but formulas for spheres in weak gradients (Datt & Elfring Reference Datt and Elfring2019) reproduce accurately the reorientation found crossing sharp gradients (Gong et al. Reference Gong, Shaik and Elfring2023).
In figure 3(a), we compare trajectories of spherical squirmers and spheroidal squirmers (
$e = 0.5$) for all three types of swimmers (
$\beta = \pm 2$ for pullers and pushers, and 
$\varepsilon = 0.1$). Spheroidal pushers still exhibit the greatest range of movement, traversing both horizontally across the gradient and vertically along it, whereas pullers cover the least distance. As expected, figure 3(a) shows that spheroidal squirmers take longer to reorient than spherical squirmers. In figure 3(b), we show the effect on a neutral squirmer as the eccentricity increases, making the spheroid more elliptical in shape, illustrating that the effect on the dynamics becomes dramatic for increasingly slender swimmers.
Figure 3. (a) Trajectories of spheroidal (
$e = 0.5$) and spherical squirmers with an initial orientation 
$\boldsymbol {p}$ orthogonal to the viscosity gradient 
$\boldsymbol{\nabla} \eta$ from 
$t=0$ to 
$t = 100a/B_1$. (b) Trajectories of neutral spheroidal squirmers of different eccentricities swimming at an initial orientation 
$\boldsymbol {p}$ orthogonal to the viscosity gradient 
$\boldsymbol{\nabla} \eta$ from 
$t=0$ to 
$t = 250a/B_1$. All squirmers eventually swim down the viscosity gradient.
We also plot, in figure 4, the trajectories of squirmers swimming in a radially varying viscosity field,
\begin{equation} \boldsymbol{\nabla} (\eta/\eta_{\infty}) = \varepsilon \boldsymbol{e}_r /a, \end{equation}as shown by Datt & Elfring (Reference Datt and Elfring2019) for spheres. Here, the assumption is that (5.3) and (5.4) still hold as a local approximation of dynamics even in radial viscosity gradients because at the particle length scale, the distinctions between the two types of gradients should be minimal. In this viscosity field, the dynamics of all three types of spheroidal squirmers again closely resembles that of spherical squirmers except that the reorientation dynamics is slowed as the squirmers become more slender. In particular, as with spheres, pushers and neutral swimmers have a stable orbit about the viscosity minimum, and as the particle becomes more slender, the radius of that orbit expands correspondingly.
Figure 4. Planar trajectories of three types of spheroidal swimmers: (a) neutral swimmers, (b) pushers, and(c) pullers, from 
$t=0$ to 
$t = 4000a/B_1$. The initial position of each swimmer is 
$x/a = 1$, 
$y/a = 1$, indicated by a red dot, with the swimmers initially pointing in the positive 
$x$-axis direction. These swimmers are placed in a radial viscosity gradient, where the viscosity increases radially outwards from the original point. The dynamics of the spheroidal squirmers qualitatively resembles that of spherical swimmers, except that the reorientation is slowed so orbits have a larger radius.
6. Disturbance viscosity effects
Up to this point we have assumed that spatial variations in viscosity are prescribed and not disturbed by the presence of the particle. However, because variations in the viscosity generally arise from variations in an underlying field that affects the viscosity, such as temperature, salt or nutrient concentration, we should take into account the effect of boundary conditions on the surface of the particle for that underlying field. For example, in an otherwise linear salt concentration field, the presence of a particle may disrupt the field (and thus the coupled viscosity field) due to salt impermeability, or in an otherwise linear temperature field, the particle may disrupt the field due to differences in thermal conductivity between the fluid and the particle. Although these disturbances diminish with distance from the particle, the disturbance does have a leading-order effect on the dynamics of the active particle (Shaik & Elfring Reference Shaik and Elfring2021).
Here, we determine the dynamics of a prolate spheroid swimmer in an otherwise constant viscosity gradient while considering the disturbance viscosity caused by a no-flux condition on the boundary of the particle, following the work of Shaik & Elfring (Reference Shaik and Elfring2021) for spheres. We write the total viscosity field as the superposition of an ambient viscosity field (denoted as 
$\eta _0$) and a disturbance viscosity field (indicated by a prime),
\begin{equation} \eta = \eta_0 + \eta', \end{equation}where the disturbance viscosity diminishes in the far-field region:
\begin{equation} \eta' \rightarrow 0 \quad \text{as} \ |\boldsymbol{r}| \rightarrow \infty. \end{equation}The transport of a scalar such as temperature or salt concentration is governed by an advection–diffusion equation. When the scalar variations are weak, the changes in viscosity are directly proportional to the changes in the underlying scalar field, hence viscosity transport is governed by a similar advection–diffusion equation. For microswimmers moving slowly in a highly diffusive scalar such as temperature or salt concentration, advection is usually small. In this limit, the distribution of viscosity satisfies Laplace's equation. As the ambient viscosity field is linear, the disturbance viscosity must also satisfy Laplace's equation:
\begin{equation} \nabla^2 \eta = \nabla^2 \eta' = 0. \end{equation}The disturbance viscosity is also determined by the boundary conditions present on the particle's surface. Here, we consider that the surface is impermeable to nutrient or salt concentration, or insulating to the temperature. In this scenario, the particle surface maintains a no-flux condition for viscosity, where
\begin{equation} \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\nabla} \eta = 0 \quad \text{on} \ S_p. \end{equation}
The detailed disturbance viscosity field is given in Appendix C (where we also give solutions with an alternative boundary condition 
$\eta (\boldsymbol {x}\in S_p)={\rm const}$). Here, we explain only the effect of disturbance viscosity on the dynamics of the active spheroid.
The impact of the total viscosity field (both ambient and disturbance viscosities) on the swimming velocity of a particle with a no-flux condition is, to leading order,
$$\begin{gather} \boldsymbol{U}_1 ={-} \frac{13 a B_2}{60}\,(\mathcal{X}^{U,nf} {{\boldsymbol{\mathsf{I}}}} - \mathcal{Y}^{U,nf} 3 \boldsymbol{pp}) \boldsymbol{\cdot} \boldsymbol{\nabla} \left(\frac{\eta_0 }{ \eta_{\infty}}\right) , \end{gather}$$
$$\begin{gather}\boldsymbol{\varOmega}_1 ={-} \frac{5}{8}\,\mathcal{X}^{\varOmega,nf} \boldsymbol{U}_{0} \times {\boldsymbol{\nabla}} \left(\frac{\eta_0 }{ \eta_{\infty}}\right), \end{gather}$$where
\begin{align} \mathcal{X}^{U,nf} & = 5 [
4e^2(63 - 117e^2 + 52 e^4 ) + 12 e (1 - e^2)^2 ({-}21 +
2e^2) \mathcal{L}_e \nonumber\\ & \quad - 9 ({-}7 + e^2) (1
- e^2)^2 \mathcal{L}_e^2 - 6e (1 - e^2)^2 \mathcal{L}_e^3 ]
\nonumber\\ & \quad \times \{ 13 e^2 [6e +
({-}3+e^2)\mathcal{L}_e][{-}2e+4e^3-({-}1+e^2)\mathcal{L}_e]
\}^{{-}1},
\end{align}
\begin{align}\mathcal{Y}^{U,nf} & = 10 [ 8e^4(48
- 77e^2 + 32 e^4 ) - 4 e^3 ( 75 - 129e^2 + 58 e^4)
\mathcal{L}_e \nonumber\\ & \quad + 2e^4 ( 33 - 56e^2 +
23e^4 ) \mathcal{L}_e^2 - e (1 - e^2)^2 ({-}33 + 37e^2)
\mathcal{L}_e^3 - 3(1 - e^2)^3 \mathcal{L}_e^4 ] \nonumber\\
&\quad \times \{ 39 e [6e +
({-}3+e^2)\mathcal{L}_e][{-}2e+4e^3-({-}1+e^2)\mathcal{L}_e][2e+({-}1+e^2)\mathcal{L}_e]
\}^{{-}1},
\end{align}
\begin{align} \mathcal{X}^{\varOmega,nf} & = \frac{2 (1 - e^2) [ 4 e^2 (7 - 8e^2) + 4e({-}7 + 7e^2 + 2e^4) \mathcal{L}_e - ({-}7 + 6e^2 + e^4) \mathcal{L}_e^2]}{5(2 - e^2) [{-}2e + 4e^3 +(1 - e^2) \mathcal{L}_e ] [2e + ({-}1+e^2)\mathcal{L}_e]}, \end{align}
are monotonically decreasing functions of the eccentricity. In the spherical limit 
$e\rightarrow 0$, we have 
$\mathcal {X}^{U,nf}=1$, 
$\mathcal {Y}^{U,nf}=1$ and 
$\mathcal {X}^{\varOmega,nf}=1$, and we recover exactly the dynamics for spheres found by Shaik & Elfring (Reference Shaik and Elfring2021). Conversely, in the slender limit 
$e\rightarrow 1$, we have 
$\mathcal {X}^{U,nf}=0$, 
$\mathcal {Y}^{U,nf}=20/39$ and 
$\mathcal {X}^{\varOmega,nf}=0$.
We see that the disturbance viscosity does not alter the fundamental physics of a spheroidal particle governed in comparison to effects of the ambient viscosity alone. It primarily increases the rate at which the particle rotates to align against the viscosity gradient. It also enhances the effects of the ambient viscosity field on various swimmer types: pushers speed up, pullers slow down, while neutral swimmers maintain consistent speeds relative to those in a homogeneous fluid.
Finally, we note that any underlying field that changes fluid viscosity will also affect the fluid density, be it temperature or salinity or the concentration of a nutrient or a polymeric additive. However, changes in viscosity can be significant without meaningful changes in density, as shown in experiments using polymers (Coppola & Kantsler Reference Coppola and Kantsler2021; Stehnach et al. Reference Stehnach, Waisbord, Walkama and Guasto2021) or glucose (Esparza López et al. Reference Esparza López, Gonzalez-Gutierrez, Solorio-Ordaz, Lauga and Zenit2021).
7. Conclusion
In this paper, we analysed the hydrodynamics of prolate spheroids, both passive and active, in constant viscosity gradients. For passive spheroids, we determined the mobility tensor that governs the dynamics of a spheroid under an external force and torque in viscosity gradients. Our analytical expression agrees with, and generalizes, previous results for spheres (Datt & Elfring Reference Datt and Elfring2019) and asymptotically slender bodies (Kamal & Lauga Reference Kamal and Lauga2023). We also derived formulas for the dynamics of active spheroids in constant viscosity gradients. These results generalize previous results for active spherical squirmers (Datt & Elfring Reference Datt and Elfring2019; Shaik & Elfring Reference Shaik and Elfring2021) to include the effects of particle shape. In general, the behaviour of spheroidal squirmers is qualitatively similar to spherical squirmers as they navigate through constant viscosity gradients. All swimmers display negative viscotaxis by reorienting to swim down viscosity gradients, except that the impact of the gradient is diminished with increasing slenderness. The viscosity difference across their body is reduced for slimmer swimmers, which leads to slower reorientation, and in the slender limit, viscotaxis ceases. The implications of this may seem limited, but it actually raises interesting new possibilities. For example, consider a swimmer that consists of a slim ‘tail’ that produces thrust but is too slender to drive reorientation in a viscosity gradient, coupled with a large spherical passive ‘head’ that strongly interacts with a viscosity gradient. In such a case, the reorientation of the swimmer would be dictated entirely by the head, and using our results for the mobility of a passive sphere (4.14) driven by a propulsive force from the tail, 
$\boldsymbol {F}_s$, we obtain 
$\boldsymbol {\varOmega } \propto \boldsymbol {F}_s\times \boldsymbol {\nabla }\eta$, indicating that such a swimmer would display positive viscotaxis by reorienting to swim up viscosity gradients in a fashion analogous to what was originally proposed by Liebchen et al. (Reference Liebchen, Monderkamp, ten Hagen and Löwen2018). Extending this idea further, one can see that geometry and activity can be tailored to control or eliminate viscotaxis. These results enrich the current understanding of how particle shape impacts viscotaxis, and the insights gleaned from this study may have implications not only for understanding the complex dynamics of natural microswimmers, but also for guiding the design and manipulation of synthetic active particles in complex fluidic systems. One can envision designing active particles that navigate gradients in temperature, or a particular chemical species (that affects fluid viscosity) acting as autonomous sensors. Our study quantifies this mechanism only for simple shapes and weak gradients, and further experimental or computational studies would be needed to really probe the limits of this mechanism.
Funding
This work was supported by the Natural Sciences and Engineering Research Council of Canada (RGPIN-2020-04850) and by a UBC Killam Accelerator Research Fellowship to G.J.E.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Spheroids in Stokes flow
Here, we give solutions to the Stokes equations for a passive and active spheroid in a Newtonian fluid with constant viscosity.
A.1. Spheroidal multipoles
Before proceeding to the solution of a passive prolate spheroid, we first introduce the spheroidal multipole solutions (Chwang & Wu Reference Chwang and Wu1975) that are used.
The Green's function 
$\boldsymbol{\mathsf{G}}$ of the Stokes equations and derivatives is given by, in component form,
\begin{equation} \left.\begin{array}{ll}
G_{ij}=\dfrac{\delta_{ij}}{r} +\dfrac{x_ix_j}{r^3}, &
\text{Stokeslet}, \\
G^d_{ijk}=G_{ij,k}=\dfrac{\delta_{jk}x_i+\delta_{ik}x_j-\delta_{ij}x_k}{r^3}
-3\,\dfrac{x_ix_j}{r^5}, & \text{dipole},\\
G^D_{ij}=G_{ij,ll}=2\,\dfrac{\delta_{ij}}{r^3}
-6\,\dfrac{x_ix_j}{r^5}, & \text{potential doublet},\\
G^R_{ijk}=\dfrac{1}{2}\,(G_{ij,k}-G_{ik,j})=\dfrac{\delta_{ik}x_j-\delta_{ij}x_k}{r^3},
& \text{rotlet},\\ G^{S}_{ijk}=\dfrac{1}{2}\,
(G_{ij,k}+G_{ik,j})=\dfrac{\delta_{kj}x_i}{r^3}-3\,\dfrac{x_ix_jx_k}{r^5},
& \text{stresslet},\\
G^Q_{ijk}=G_{ij,llk}={-}6\,\dfrac{\delta_{jk}x_i+\delta_{ik}x_j+\delta_{ij}x_k}{r^5}
+30\, \dfrac{x_ix_jx_k}{r^7}, & \text{potential
quadrupole}. \end{array}\right\}
\end{equation} Spheroidal multipoles are a weighted distribution of the above multipoles between the foci 
$\xi =-c$ and 
$c$, where 
$c = ae$, used to represent flows around spheroidal particles:
$$\begin{gather} Q_{ij}=\int^{c}_{{-}c} G_{ij}(\boldsymbol{x}-\xi \boldsymbol{p}) \,{\rm d}\xi, \end{gather}$$
$$\begin{gather}Q^D_{ij}=\int^{c}_{{-}c} (c^2-\xi^2)\,G^D_{ij}(\boldsymbol{x}-\xi \boldsymbol{p}) \,{\rm d}\xi , \end{gather}$$
$$\begin{gather}Q^R_{ijk}=\int^{c}_{{-}c} (c^2-\xi^2)\,G^R_{ijk}(\boldsymbol{x}-\xi \boldsymbol{p}) \,{\rm d}\xi, \end{gather}$$
$$\begin{gather}Q^{S}_{ijk}=\int^{c}_{{-}c}(c^2-\xi^2)\,G^{S}_{ijk}(\boldsymbol{x}-\xi \boldsymbol{p}) \,{\rm d}\xi , \end{gather}$$
$$\begin{gather}Q^Q_{ijk}=\int^{c}_{{-}c} (c^2-\xi^2)^2\,G^Q_{ijk}(\boldsymbol{x}-\xi \boldsymbol{p}) \,{\rm d}\xi. \end{gather}$$Explicit expressions for spheroidal multipoles are taken from Einarsson et al. (Reference Einarsson, Candelier, Lundell, Angilella and Mehlig2015) and Abtahi & Elfring (Reference Abtahi and Elfring2019):
\begin{align} Q_{ij}&=\delta_{ij} I^0_1+x_ix_j I^0_3 -(x_ip_j+x_jp_i) I^1_3+ p_ip_jI^2_3, \end{align}
\begin{align} Q^D_{ij}&=2\delta_{ij}J^0_3+6 [{-}x_ix_jJ^0_5+(x_ip_j+x_jp_i)J_5^1-p_ip_jJ^2_5], \end{align}
\begin{align}
Q^R_{ijk}&=(\delta_{ik}x_j-\delta_{ij}x_k)J^0_3+(\delta_{ij}p_k-\delta_{ik}p_j)J^1_3,
\end{align}
\begin{align}Q^{S}_{ijk}&=\delta_{jk}x_iJ^0_3-\delta_{jk}p_iJ^1_3 \nonumber\\
&\quad+3 [{-}x_ix_jx_kJ^0_5+(x_ix_kp_j+x_jx_kp_i+x_ix_jp_k)J^1_5\nonumber\\
&\quad-(x_kp_ip_j+x_ip_jp_k+x_jp_ip_k)J^2_5+p_ip_jp_kJ^3_5],
\end{align}
\begin{align} Q^Q_{ijk}&=6 [-(\delta_{jk}x_i+\delta_{ik}x_j+\delta_{ij}x_k)K^0_5+(\delta_{jk}p_i+\delta_{ik}p_j+\delta_{ij}p_k)K^1_5 ]\nonumber\\ &\quad+30 [ x_ix_jx_kK^0_7-(x_ix_kp_j+x_jx_kp_i+x_ix_jp_k)K^1_7 \nonumber\\ &\quad+(x_kp_ip_j+x_ip_jp_k+x_jp_ip_k)K^2_7-p_ip_jp_kK^3_7 ], \end{align}
\begin{align} Q_{ij,k} & = (- \delta_{ij} x_k + \delta_{ik} x_j + \delta_{jk} x_i) I^0_3 + (\delta_{ij} p_k - \delta_{ik} p_j - \delta_{jk} p_i) I^1_3 \nonumber\\ & \quad + 3 [ - x_i x_j x_k I^0_5 + (x_ix_kp_j+x_jx_kp_i+x_ix_jp_k) I^1_5 ], \end{align}
\begin{align} Q^D_{ij,k} & = 6 [ - ( \delta_{ij} x_k + \delta_{ik} x_j + \delta_{jk} x_i) J^0_5 + ( \delta_{ij} p_k + \delta_{ik} p_j + \delta_{jk} p_i) J^1_5 ] \nonumber\\ & \quad + 30 [ x_i x_j x_k J^0_7 - (x_ix_kp_j+x_jx_kp_i+x_ix_jp_k) J^1_7 \nonumber\\ & \quad + (x_kp_ip_j+x_ip_jp_k+x_jp_ip_k) J^1_7 - p_ip_jp_k J^3_7 ], \end{align}
\begin{align} Q^{R}_{ijk,m}&=(\delta_{ik}\delta_{jm}-\delta_{ij}\delta_{km})J^0_3 \nonumber\\ & \quad+3(\delta_{ik}x_j-\delta_{ij}x_k)(\kern 0.06em p_m J^1_5-x_m J^0_5)\nonumber\\ & \quad+3(\delta_{ij}p_k-\delta_{ik}p_j)(p_m J^2_5-x_m J^1_5), \end{align}
\begin{align} Q^{S}_{ijk,m}&=\delta_{jk}\delta_{im}J^0_3+3\delta_{jk}x_i (p_mJ^1_5-x_mJ^0_5) -3\delta_{jk}p_i(p_mJ^2_5-x_mJ^1_5) \nonumber\\ &\quad+3 [ -(\delta_{im}x_jx_k+\delta_{jm}x_ix_k+\delta_{km}x_ix_j)J^0_5-5x_ix_jx_k(p_mJ^1_7-x_mJ^0_7) \nonumber\\ &\quad+(\delta_{im}x_kp_j+\delta_{km}x_ip_j + \delta_{jm} x_kp_i+\delta_{km} x_jp_i + \delta_{im} x_jp_k+\delta_{jm}x_ip_k)J^1_5 \nonumber\\ &\quad+5(x_ix_kp_j+x_jx_kp_i+x_ix_jp_k)(p_mJ^2_7-x_mJ^1_7) \nonumber\\ &\quad-(\delta_{km}p_ip_j+\delta_{im}p_jp_k+\delta_{jm}p_ip_k)J^2_5 +5p_ip_jp_k(p_mJ^4_7-x_mJ^3_7) \nonumber\\ &\quad-5(x_kp_ip_j+x_ip_jp_k+x_jp_ip_k)(p_mJ^3_7-x_mJ^2_7)], \end{align}
\begin{align} Q^Q_{ijk,m}&=6[-(\delta_{jk}\delta_{im}+\delta_{ik}\delta_{jm}+\delta_{ij}\delta_{km})K^0_5-5(\delta_{jk}x_i+\delta_{ik}x_j+\delta_{ij}x_k)(p_mK^1_7-x_mK^0_7)\nonumber\\ &\quad+5(\delta_{jk}p_i+\delta_{ik}p_j+\delta_{ij}p_k)(p_mK^2_7-x_mK^1_7) ]\nonumber\\ &\quad+30 [ (\delta_{im}x_jx_k+\delta_{jm}x_ix_k+\delta_{km}x_ix_j)K^0_7+7x_ix_jx_k(p_mK^1_9-x_mK^0_9)\nonumber\\ &\quad-(\delta_{im}x_kp_j+\delta_{km}x_ip_j+\delta_{jm}x_kp_i+\delta_{km}x_jp_i+\delta_{im}x_jp_k+\delta_{jm}x_ip_k)K^1_7\nonumber\\ &\quad-7(x_ix_kp_j+x_jx_kp_i+x_ix_jp_k)(p_mK^2_9-x_mK^1_9)\nonumber\\ &\quad+(\delta_{km}p_ip_j+\delta_{im}p_jp_k+\delta_{jm}p_ip_k)K^2_7 -p_ip_jp_kK^3_7\nonumber\\ &\quad+7(x_kp_ip_j+x_ip_jp_k+x_jp_ip_k)(p_mK^3_9-x_mK^2_9) ], \end{align}where
$$\begin{gather} I^n_m =\int^{c}_{{-}c} {\rm d}\xi\,\frac{\xi^n}{|\boldsymbol{x}-\xi \boldsymbol{p}|^m}, \end{gather}$$
$$\begin{gather}J^n_m =c^2 I^n_m -I^{n+2}_m, \end{gather}$$
$$\begin{gather}K^n_m=c^2 J^n_m -J^{n+2}_m=c^4 I^n_m - 2 c^2I^{n+2}_m +I^{n+4}_m. \end{gather}$$ The integrals 
$I^n_{m}$ satisfy the relationship
\begin{equation} \frac{\partial}{\partial x_i} I^n_{m} = mp_i I^{n+1}_{m+2} - mx_i I^{n}_{m+2}. \end{equation}
To simplify integration, one may employ an auxiliary coordinate system 
$(x',y',z')$, with 
$x'$ aligned with 
$\boldsymbol {p}$ such that
\begin{equation} I^n_m=\int^c_{{-}c} {\rm d}\xi\,\frac{\xi^n}{[(x'-\xi)^2+(y')^2+(z')^2]^{m/2}}=\int^c_{{-}c} {\rm d}\xi\, \frac{\xi^n}{[(x'-\xi)^2+R^2]^{m/2}}, \end{equation}where on the surface of the particle we have
\begin{equation} \left.\begin{array}{c} R_1=\sqrt{(x'+c)^2+R^2},\\[3pt] R_2=\sqrt{(x'-c)^2+R^2},\\[3pt] R=\sqrt{(1-e^2) (a^2-x'^2)}. \end{array}\right\} \end{equation}The integrals also satisfy the relationship
\begin{equation} I^n_{m} = x' I^{n-1}_m+\frac{(n-1) I^{n-2}_{m-2}}{m-2}-\frac{c^{n-1} (({-}1)^{n} R_1^{2-m}+R_2^{2-m})}{m-2}. \end{equation}
Integrals 
$J^n_m$ and 
$K^n_m$ can be calculated easily from (A18) and (A19).
A.2. A passive prolate spheroid
A.2.1. Rigid-body translation
The flow field due to a prolate spheroid translating with velocity 
$\hat {\boldsymbol {U}}$ in a quiescent fluid is
\begin{equation} \hat{u}_{i}= (Q_{ij} + \alpha_1 Q^D_{ij}) [\mathcal{A}^U p_j p_m + \mathcal{B}^U (\delta_{jm}-p_jp_m) ] \hat{U}_{m}, \end{equation}where
\begin{equation}
\left.\begin{array}{c} \alpha_1=\dfrac{1-e^2}{4e^2},\\
\mathcal{A}^U=\dfrac{e^2}{ - 2e + (1 +
e^2)\mathcal{L}_e},\\ \mathcal{B}^U=\dfrac{2e^2}{ 2e +
(3e^2 - 1 )\mathcal{L}_e}. \end{array}\right\}
\end{equation}The strain-rate tensor can be written as
\begin{equation} \hat{\dot{\boldsymbol{\gamma}}} = 2 \hat{\boldsymbol{E}}_{\boldsymbol{U}} \boldsymbol{\cdot} \hat{\boldsymbol{U}}, \end{equation}where
\begin{equation} \hat{E}_{U_{ikm}} =\tfrac{1}{2} (Q^T_{ijk} + \alpha_1 Q^{DT}_{ijk}) [\mathcal{A}^U p_j p_m + \mathcal{B}^U (\delta_{jm}-p_jp_m) ], \end{equation}
and 
$Q^{T}_{ijk}=Q_{ij,k}+Q_{kj,i}$. Here, 
$Q^{DT}_{ijk}$ is defined similarly to 
$Q^T_{ijk}$.
A.2.2. Rigid-body rotation
The flow field due to a prolate spheroid rotating with angular velocity 
$\hat {\boldsymbol {\varOmega }}$ in another quiescent fluid is
\begin{align} \hat{u}_{i}&=\{-\epsilon_{jkl} Q^R_{ijk} [ \mathcal{A}^{\varOmega} p_l p_s +\mathcal{B}^{\varOmega} (\delta_{ls}-p_lp_s)] \nonumber\\ &\quad +( Q^{S}_{ijk}+ \alpha_2 Q^Q_{ijk} ) \mathcal{C}^{\varOmega} (\epsilon_{jsm}p_kp_m+\epsilon_{ksm}p_jp_m) \} \hat{\varOmega}_{s}, \end{align}where
\begin{equation} \left.\begin{array}{c} \alpha_2=\dfrac{1-e^2}{8e^2},\\ \mathcal{A}^{\varOmega}=\dfrac{1 - e^2}{-4 e + 2 (1 - e^2) \mathcal{L}_e},\\ \mathcal{B}^{\varOmega}=\dfrac{2 - e^2}{4 e - 2 (1 + e^2) \mathcal{L}_e},\\ \mathcal{C}^{\varOmega}=\dfrac{e^2}{ 4e - 2 (1 + e^2) \mathcal{L}_e}. \end{array}\right\} \end{equation}The strain-rate tensor can be written as
\begin{equation} \hat{\dot{\boldsymbol{\gamma}}} = 2 \hat{\boldsymbol{E}}_{\boldsymbol{\varOmega}} \boldsymbol{\cdot} \hat{\boldsymbol{\varOmega}}, \end{equation}where
\begin{align} \hat{E}_{\varOmega_{ims}} &=\tfrac{1}{2} \{ -\epsilon_{jkl} Q^{RT}_{ijkm} [\mathcal{A}^{\varOmega} p_l p_s +\mathcal{B}^{\varOmega} (\delta_{ls}-p_lp_s) ]\nonumber\\ &\quad +( Q^{ST}_{ijkm}+ \alpha Q^{QT}_{ijkm} ) \mathcal{C}^{\varOmega} (\epsilon_{jsm}p_kp_m+\epsilon_{ksm}p_jp_m)\}, \end{align}
and 
$Q^{RT}_{ijkm}=Q^{R}_{ijk,m}+Q^{R}_{mjk,i}$. Here, 
$Q^{ST}_{ijkm}$ and 
$Q^{QT}_{ijkm}$ are defined similarly to 
$Q^{RT}_{ijkm}$.
Finally, the tensor 
$\hat{{{\boldsymbol{\mathsf{E}}}}}_{{{\boldsymbol{\mathsf{U}}}}}$ used in the integral (3.3) is simply defined as
\begin{equation} \hat{{{\boldsymbol{\mathsf{E}}}}}_{{{\boldsymbol{\mathsf{U}}}}} = \begin{pmatrix} \hat{\boldsymbol{E}}_{\boldsymbol{U}}\\ \hat{\boldsymbol{E}}_{\boldsymbol{\varOmega}} \end{pmatrix}. \end{equation}A.2.3. Mobility tensor for a prolate spheroid in a Newtonian fluid
The mobility tensor
\begin{equation} \hat{{{\boldsymbol{\mathsf{M}}}}}_{{{\boldsymbol{\mathsf{U}}}} {{\boldsymbol{\mathsf{F}}}}} = \hat{{{\boldsymbol{\mathsf{R}}}}}^{{-}1}_{{{\boldsymbol{\mathsf{F}}}} {{\boldsymbol{\mathsf{U}}}}} =\begin{pmatrix} \hat{{{\boldsymbol{\mathsf{M}}}}}_{\boldsymbol{U} \boldsymbol{F}} & \hat{{{\boldsymbol{\mathsf{M}}}}}_{\boldsymbol{U} \boldsymbol{L}} \\ \hat{{{\boldsymbol{\mathsf{M}}}}}_{\boldsymbol{\varOmega} \boldsymbol{F}} & \hat{{{\boldsymbol{\mathsf{M}}}}}_{\boldsymbol{\varOmega} \boldsymbol{L}} \end{pmatrix} \end{equation}
couples force and torque to rigid-body translation and rotation for a body in Stokes flows. For a prolate spheroid in a Newtonian fluid with constant viscosity 
$\eta _{\infty }$, there is no torque-translation (or force-rotation) coupling. Specifically, the terms are (Kim & Karilla Reference Kim and Karilla1991)
\begin{equation} \left.\begin{array}{c} \displaystyle\hat{{{\boldsymbol{\mathsf{M}}}}}_{\boldsymbol{U} \boldsymbol{F}} = \dfrac{1}{6 {\rm \pi}\eta_{\infty} a} \left[ \dfrac{1}{\mathcal{X}^{A}}\,\boldsymbol{pp} + \dfrac{1}{\mathcal{Y}^{A}}\,({{\boldsymbol{\mathsf{I}}}} - \boldsymbol{pp}) \right],\\ \displaystyle\hat{{{\boldsymbol{\mathsf{M}}}}}_{\boldsymbol{\varOmega} \boldsymbol{L}} = \dfrac{1}{8 {\rm \pi}\eta_{\infty} a^3} \left[ \dfrac{1}{\mathcal{X}^{C}}\,\boldsymbol{pp} + \dfrac{1}{\mathcal{Y}^{C}}\,({{\boldsymbol{\mathsf{I}}}} - \boldsymbol{pp}) \right],\\ \hat{{{\boldsymbol{\mathsf{M}}}}}_{\boldsymbol{U} \boldsymbol{L}} = \hat{{{\boldsymbol{\mathsf{M}}}}}_{\boldsymbol{\varOmega} \boldsymbol{F}} = {{\boldsymbol{\mathsf{0}}}}, \end{array}\right\} \end{equation} where 
$\mathcal {X}^{A}$, 
$\mathcal {Y}^{A}$, 
$\mathcal {X}^{C}$ and 
$\mathcal {Y}^{C}$ are functions of eccentricity 
$e$:
\begin{equation} \left.\begin{array}{c} \mathcal{X}^A=\dfrac{8 e^3}{3[ - 2e + (1 + e^2)\mathcal{L}_e]},\\ \mathcal{Y}^A=\dfrac{16 e^3}{3[ 2e + (3 e^2 - 1)\mathcal{L}_e]},\\ \mathcal{X}^C=\dfrac{4e^3(1 - e^2)}{3 [ 2e - ( 1 - e^2)\mathcal{L}_e ]},\\ \mathcal{Y}^C=\dfrac{4e^3(2 - e^2)}{ 3 [{-}2e + ( 1 + e^2)\mathcal{L}_e ]}. \end{array}\right\} \end{equation}A.2.4. Extra stress tensor
In (4.4), we defined the tensor
\begin{equation} \boldsymbol{{{\boldsymbol{\mathsf{R}}}}}_{NN} =\int_{\mathcal{V}}2(\eta (\boldsymbol{x}) -\eta_{\infty})\hat{{{\boldsymbol{\mathsf{E}}}}}_{{{\boldsymbol{\mathsf{U}}}}}:\hat{{{\boldsymbol{\mathsf{E}}}}}_{{{\boldsymbol{\mathsf{U}}}}}\, {\rm d}V. \end{equation}Writing
\begin{equation} \boldsymbol{{{\boldsymbol{\mathsf{R}}}}}_{NN} = \begin{pmatrix} {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{F} \boldsymbol{U}} & {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{F} \boldsymbol{\varOmega}} \\ {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{L} \boldsymbol{U}} & {{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{L} \boldsymbol{\varOmega}} \end{pmatrix}, \end{equation}we have
\begin{equation} \left.\begin{array}{c} \displaystyle{{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{F} \boldsymbol{U}} = \int_{\mathcal{V}} 2 (\eta (\boldsymbol{x}) - \eta_{\infty}) \hat{\boldsymbol{E}}_{\boldsymbol{U}} : \hat{\boldsymbol{E}}_{\boldsymbol{U}} \, {\rm d} V,\\ \displaystyle{{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{F} \boldsymbol{\varOmega}} = \int_{\mathcal{V}} 2 (\eta (\boldsymbol{x}) - \eta_{\infty}) \hat{\boldsymbol{E}}_{\boldsymbol{U}} : \hat{\boldsymbol{E}}_{\boldsymbol{\varOmega}} \, {\rm d} V,\\ \displaystyle{{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{L} \boldsymbol{U}} = \int_{\mathcal{V}} 2 (\eta (\boldsymbol{x}) - \eta_{\infty}) \hat{\boldsymbol{E}}_{\boldsymbol{\varOmega}} : \hat{\boldsymbol{E}}_{\boldsymbol{U}} \, {\rm d} V,\\ \displaystyle{{\boldsymbol{\mathsf{R}}}}_{\boldsymbol{L} \boldsymbol{\varOmega}} = \int_{\mathcal{V}} 2 (\eta (\boldsymbol{x}) - \eta_{\infty}) \hat{\boldsymbol{E}}_{\boldsymbol{\varOmega}} : \hat{\boldsymbol{E}}_{\boldsymbol{\varOmega}} \, {\rm d} V, \end{array}\right\} \end{equation} and 
$\boldsymbol{\mathsf{R}}_{\boldsymbol {F} \boldsymbol {\varOmega }} = \boldsymbol{\mathsf{R}}_{\boldsymbol {L} \boldsymbol {U}}^\textrm {T}$.
A.3. An active prolate spheroid in Stokes flow
The flow field 
$\boldsymbol {u}_0$ around a two-mode spheroidal squirmer swimming in a Newtonian fluid with constant viscosity can be written in terms of a stream function 
$\psi _0$ (Keller & Wu Reference Keller and Wu1977; Theers et al. Reference Theers, Westphal, Gompper and Winkler2016; van Gogh et al. Reference van Gogh, Demir, Palaniappan and Pak2022) using a particle-aligned coordinate system 
$O'\zeta _1 \zeta _2 \phi$ (see Appendix B for details), as
\begin{equation} \boldsymbol{u}_0 = \frac{1}{h_{\zeta_2} h_{\phi}}\,\frac{\partial \psi_0}{\partial \zeta_2}\,\boldsymbol{e}_{\zeta_1} - \frac{1}{h_{\zeta_1} h_{\phi}}\,\frac{\partial \psi_0}{\partial \zeta_1}\,\boldsymbol{e}_{\zeta_2}, \end{equation}where
\begin{align} \psi_0 & = C_1\,H_2 (\zeta_1)\,G_2 (\zeta_2) + C_2 \zeta_1 ( 1 - \zeta_2^2) \nonumber\\ & \quad + C_3\,H_3 (\zeta_1)\,G_3 (\zeta_2) + C_4 \zeta_2 ( 1 - \zeta_2^2) + \tfrac{1}{2} U_0 c^2 (\zeta_1^2 - 1)(1 - \zeta_2^2). \end{align}
Here, 
$H_n(x)$ and 
$G_n(x)$ are Gegenbauer functions of the first and second order of degree 
$-1/2$ (Theers et al. Reference Theers, Westphal, Gompper and Winkler2016). The coefficients 
$C_n$ are
\begin{equation} \left.\begin{array}{c} C_1 = 2c^2\,\dfrac{U_0 (\tilde{\zeta}_1^2 + 1) - 2 B_1 \tilde{\zeta}_1^2}{- \tilde{\zeta}_1 + (1 + \tilde{\zeta}_1^2 ) \coth^{{-}1} \tilde{\zeta}_1},\\ C_2 = c^2\,\dfrac{B_1 \tilde{\zeta}_1 [\tilde{\zeta}_1 - (\tilde{\zeta}_1^2 - 1) \coth^{{-}1} \tilde{\zeta}_1] - U_0}{- \tilde{\zeta}_1 + (1 + \tilde{\zeta}_1^2) \coth^{{-}1} \tilde{\zeta}_1},\\ C_3 = c^2\,\dfrac{ 4 B_2 \tilde{\zeta}_1}{3 \tilde{\zeta}_1 + (1 - 3\tilde{\zeta}_1^2) \coth^{{-}1} \tilde{\zeta}_1}, \\ C_4 = c^2\,\dfrac{ B_2 \tilde{\zeta}_1 [2/3 - \tilde{\zeta}_1^2 + \tilde{\zeta}_1(\tilde{\zeta}_1^2 - 1) \coth^{{-}1} \tilde{\zeta}_1]}{3 \tilde{\zeta}_1 + (1 - 3\tilde{\zeta}_1^2) \coth^{{-}1} \tilde{\zeta}_1} \end{array}\right\} \end{equation} and 
$\tilde {\zeta }_1 = 1 / e$.
Appendix B. Coordinate transformation
We choose an arbitrary point 
$O$ and construct a lab-frame Cartesian coordinate system with unit vectors 
$\boldsymbol {e}_i$ (
$i = 1, 2, 3$) and position vector 
$\boldsymbol {x} = x\boldsymbol {e}_1 + y\boldsymbol {e}_2 + z\boldsymbol {e}_3$. The centre of the particle can be expressed as 
$\boldsymbol {x}_c = x_c\boldsymbol {e}_1 + y_c\boldsymbol {e}_2 + z_c\boldsymbol {e}_3$. Without loss of generality, we can always adjust the axes to make sure that the ambient viscosity varies only in the 
$\boldsymbol {e}_1$ direction. However, the volume integrations in the reciprocal theorem are difficult to evaluate analytically in the lab-frame coordinate system 
$Oxyz$. To solve this problem, we use a particle-aligned Cartesian coordinate system 
$O'XYZ$ and the corresponding spheroidal coordinate system 
$O'\zeta _1 \zeta _2 \phi$, where 
$O'$ is the centre of the spheroid at 
$\boldsymbol {x}_c$ (figure 5). The Cartesian coordinate axes are determined by the viscosity gradient direction 
$\boldsymbol {d} = {{\boldsymbol{\nabla} } \eta }/{|{\boldsymbol{\nabla} } \eta |}$ and the swimming direction 
$\boldsymbol {p}$. The unit vectors are
\begin{equation} \left.\begin{array}{c} \boldsymbol{e}_X = \dfrac{(\boldsymbol{d} \times \boldsymbol{p}) \times \boldsymbol{p}}{|(\boldsymbol{d} \times \boldsymbol{p}) \times \boldsymbol{p}|},\\ \boldsymbol{e}_Y = \dfrac{\boldsymbol{d} \times \boldsymbol{p}}{|\boldsymbol{d} \times \boldsymbol{p}|},\\ \boldsymbol{e}_Z = \boldsymbol{p}. \end{array}\right\} \end{equation}
Figure 5. The particle-aligned Cartesian coordinate system 
$O'XYZ$ and prolate spheroid coordinate system 
$O'\zeta _1 \zeta _2 \phi$.
When 
$\boldsymbol {d}$ is parallel or anti-parallel to 
$\boldsymbol {p}$, we can simply set 
$\boldsymbol {e}_X = \boldsymbol {e}_1$, 
$\boldsymbol {e}_Y = \boldsymbol {e}_2$ and 
$\boldsymbol {e}_Z = \boldsymbol {e}_3$ without loss of generality. The position vector in this coordinate system is 
$\boldsymbol {r} = \boldsymbol {x} - \boldsymbol {x}_c = X \boldsymbol {e}_X + Y \boldsymbol {e}_Y + Z \boldsymbol {e}_Z$. We then write the viscosity field as
\begin{equation} \eta = \eta_{\infty} + \varepsilon\,\frac{\eta_{\infty}}{a}\,(x_c + \boldsymbol{r} \boldsymbol{\cdot} \boldsymbol{e}_1), \end{equation}
where 
$\boldsymbol {e}_1$ is obtained by inverting (B1).
In the particle-aligned Cartesian coordinate system, the surface of a spheroid satisfies
\begin{equation} \frac{Z^2}{a^2} + \frac{X^2 + Y^2}{b^2} = 1. \end{equation}
Cartesian coordinates 
$(X, Y, Z)$ can be written in terms of 
$(\zeta _1, \zeta _2, \phi )$ as
\begin{equation} \left.\begin{array}{c} X = c \sqrt{\zeta_1^2 - 1} \sqrt{ 1 - \zeta_2^2} \cos \phi, \\ Y = c \sqrt{\zeta_1^2 - 1} \sqrt{ 1 - \zeta_2^2} \sin \phi, \\ Z = c \zeta_1 \zeta_2, \end{array}\right\} \end{equation} where 
$1 \leq \zeta _1 < \infty$, 
$-1 \leq \zeta _2 \leq 1$ and 
$0 \leq \phi < 2{\rm \pi}$. Here, 
$c = \sqrt {a^2 - b^2}$ is half of the focal length. The unit normal vector and the tangent vector to the particle are 
$\boldsymbol {n} = \boldsymbol {e}_{\zeta _1}$ and 
$\boldsymbol {s} = -\boldsymbol {e}_{\zeta _2}$. Scale factors are
\begin{equation} \left.\begin{array}{c} h_{\zeta_1} = c \sqrt{\dfrac{\zeta_1^2 - \zeta_2^2}{\zeta_1^2 - 1}},\\ h_{\zeta_2} = c \sqrt{\dfrac{\zeta_1^2 - \zeta_2^2}{1 - \zeta_2^2}},\\ h_{\phi} = c \sqrt{\zeta_1^2 - 1} \sqrt{ 1 - \zeta_2^2}. \end{array}\right\} \end{equation}Appendix C. Disturbance viscosity field
The general solution of (6.3) in spheroidal coordinates satisfies
\begin{equation} \eta ' = \sum_{k = 0}^{ \infty} \sum_{m=k}^{\infty} [A_{k,m} \cos (m \phi) + B_{k,m} \sin (m \phi)]\,P_k^m (\zeta_2)\,Q_k^m (\zeta_1), \end{equation}
where 
$A_{k,m}$, 
$B_{k,m}$ are the constant coefficients, while 
$P_k^m$ and 
$Q_k^m$ are the associated Legendre polynomials of the first and second kind, respectively, with 
$k$ the degree and 
$m$ the order. Mathematical expressions for 
$P_k^m$ and 
$Q_k^m$ can be found in Abramowitz & Stegun (Reference Abramowitz and Stegun1964). Below, we determine the coefficients for first a no-flux boundary and then a constant viscosity boundary condition.
C.1. No flux
Supposing that the ambient viscosity field is aligned with 
$\boldsymbol {e}_1$, we can write the no-flux constraint in particle-aligned coordinates as
\begin{equation} \left.\frac{\partial \eta '}{\partial \zeta_1} \right|_{\zeta_1 = \tilde{\zeta}_1} ={-} \varepsilon \eta_{\infty} p_1 e \zeta_2 + \frac{\varepsilon \eta_{\infty} \sqrt{1 - p_1^2}}{e \sqrt{1 - e^2}}\,\sqrt{1 - \zeta_2^2} \cos \phi, \end{equation}
where 
$p_1 = \boldsymbol {p} \boldsymbol {\cdot } \boldsymbol {e}_1$. The expression for the disturbance viscosity field is
\begin{equation} \eta ' = A_{1,0}\,P_1^0 (\zeta_2)\,Q_1^0 (\zeta_1) + A_{1,1}\,P_1^1 (\zeta_2)\,Q_1^1 (\zeta_1) \cos (\phi), \end{equation}where
$$\begin{gather} A_{1,0} = \varepsilon \eta_{\infty} p_1\,\frac{2e (1 - e^2)}{2e - (1 - e^2) \mathcal{L}_e}, \end{gather}$$
$$\begin{gather}A_{1,1} = \varepsilon \eta_{\infty} \sqrt{1 - p_1^2}\,\frac{2e (1 - e^2)}{2e - 4 e^3 - (1 - e^2) \mathcal{L}_e}. \end{gather}$$The changes in the translational and rotational velocity due to the disturbance viscosity are
$$\begin{gather} \boldsymbol{U}_1' ={-} \frac{a B_2}{60}\,(\mathcal{X}^{U',nf} \boldsymbol{I} - \mathcal{Y}^{U',nf} 3 \boldsymbol{pp}) \boldsymbol{\cdot} {\boldsymbol{\nabla}} \left(\frac{\eta_0 }{ \eta_{\infty}}\right), \end{gather}$$
$$\begin{gather}\boldsymbol{\varOmega}_1' ={-} \frac{1}{8}\,\mathcal{X}^{\varOmega',nf} \boldsymbol{U}_{0} \times {\boldsymbol{\nabla}} \left(\frac{\eta_0 }{ \eta_{\infty}}\right), \end{gather}$$where
\begin{align} \mathcal{X}^{U',nf} & = 5 ({-}1 + e^2) [ 8e^3({-}9 -33e^2 + 32 e^4 ) - 4 e^2 ( - 27 - 39e^2 + 62e^4) \mathcal{L}_e\nonumber\\ & \quad + 6e ({-}9 + 7e^2 -3e^4 + 5e^6) \mathcal{L}_e^2 - 3 ({-}3 + 9e^2 - 13e^4 +7 e^6) \mathcal{L}_e^3 ]\nonumber\\ & \quad \times \{ 2 e^5 [6e + ({-}3+e^2)\mathcal{L}_e][{-}2e+4e^3-({-}1+e^2)\mathcal{L}_e] \}^{{-}1}, \end{align}
\begin{align} \mathcal{Y}^{U',nf} & = 5 ({-}1 + e^2) [ 16e^4({-}27 + 54e^2 - 102 e^4 + 64 e^6 ) \nonumber\\ & \quad - 16 e^3 ( - 54 + 117e^2 -135 e^4 +74e^6 ) \mathcal{L}_e \nonumber\\ & \quad + 8e^2 ({-}81 + 189e^2 - 150e^4 + 31e^6 + 13 e^8) \mathcal{L}_e^2 \nonumber\\ & \quad - 4e ({-}54 + 135e^2 - 99e^4 + e^6 + 17 e^8) \mathcal{L}_e^3 - 3(3-4e^2+e^4)^2 \mathcal{L}_e^4 ]\nonumber\\ & \quad \times \{ 6 e^5 [6e + ({-}3+e^2)\mathcal{L}_e][{-}2e+4e^3-({-}1+e^2)\mathcal{L}_e][2e+({-}1+e^2)\mathcal{L}_e] \}^{{-}1}, \end{align}
\begin{align} \mathcal{X}^{\varOmega',nf} & = \frac{2 (1 - e^2) [ 4 e^2 (5 - 4e^2) + 20e({-}1+e^2) \mathcal{L}_e + (5 - 6e + e^4) \mathcal{L}_e^2]}{(2 - e^2) [{-}2e + 4e^3 + (1 - e^2) \mathcal{L}_e ] [2e + ({-}1+e^2)\mathcal{L}_e]}. \end{align}Adding (C6) and (C7) to (5.3) and (5.4), one obtains (6.5) and (6.6).
C.2. Constant viscosity
The constant viscosity constraint 
$\eta (x \in S_p) = \eta _p = const$ in spheroidal coordinates is
\begin{equation} \eta' |_{\zeta_1 = \tilde{\zeta}_1} = \eta_c - \varepsilon \eta_{\infty} p_1 \zeta_2 + \varepsilon \eta_{\infty} \sqrt{(1 - p_1^2)(1 - e^2)(1 - \zeta_2^2)} \cos \phi, \end{equation}
where 
$\eta _c = \eta _p - \eta _{\infty } - \varepsilon ({\eta _{\infty }}/{a}) x_c$ is a constant, while the other term varies on the surface of the spheroid. The disturbance viscosity field satisfying this constraint is
\begin{equation} \eta ' = \acute{A}_{0,0}\,P_0^0 (\zeta_2)\,Q_0^0 (\zeta_1) + \acute{A}_{1,0}\,P_1^0 (\zeta_2)\,Q_1^0 (\zeta_1) + \acute{A}_{1,1}\,P_1^1 (\zeta_2)\,Q_1^1 (\zeta_1) \cos (\phi), \end{equation}where
$$\begin{gather} \acute{A}_{0,0} = \frac{2 \eta_c}{\mathcal{L}_e}, \end{gather}$$
$$\begin{gather}\acute{A}_{1,0} = \varepsilon \eta_{\infty} p_1\,\frac{2e}{ 2e - \mathcal{L}_e}, \end{gather}$$
$$\begin{gather}\acute{A}_{1,1} = \varepsilon \eta_{\infty} \sqrt{1 - p_1^2}\,\frac{2e (1 - e^2)}{2e - (1 - e^2) \mathcal{L}_e}. \end{gather}$$The changes in the translational and rotational velocity due to the disturbance viscosity are then
$$\begin{gather} \boldsymbol{U}_1' = \frac{\eta_c}{12 \varepsilon \eta_{\infty}}\,\mathcal{Z}^{U'} \boldsymbol{U}_0 + \frac{a B_2}{30}\,(\mathcal{X}^{U',c} {{\boldsymbol{\mathsf{I}}}} - \mathcal{Y}^{U',c} 3 \boldsymbol{pp}) \boldsymbol{\cdot} {\boldsymbol{\nabla}} \left(\frac{\eta_0 }{ \eta_{\infty}}\right), \end{gather}$$
$$\begin{gather}\boldsymbol{\varOmega}_1' = \frac{1}{4}\,\mathcal{X}^{\varOmega',c} \boldsymbol{U}_{N} \times {\boldsymbol{\nabla}} \left(\frac{\eta_0 }{ \eta_{\infty}}\right), \end{gather}$$where
\begin{align} \mathcal{Z}^{U'} & =
\frac{6 (1 - e^2) (2 e - \mathcal{L}_e)^2}{e^2
\mathcal{L}_e [2e - (1 - e^2 ) \mathcal{L}_e]},
\end{align}
\begin{align} \mathcal{X}^{U',c} & = 5 ({-}1 + e^2) [ 8e^3({-}9 -33e^2 +
32 e^4 ) - 4 e^2 ( - 27 - 39e^2 + 62e^4)
\mathcal{L}_e\nonumber\\ & \quad + 6e ({-}9 + 7e^2 -3e^4 +
5e^6) \mathcal{L}_e^2 - 3 ({-}3 + 9e^2 - 13e^4 +7 e^6)
\mathcal{L}_e^3 ]\nonumber\\ & \quad \times \{ 4 e^5
[6e + ({-}3+e^2)\mathcal{L}_e][2e+({-}1+e^2)\mathcal{L}_e]
\}^{{-}1},
\end{align}
\begin{align} \mathcal{Y}^{U',c} & = 5 [ 16e^4(27 - 27e^2 - 33 e^4 + 32 e^6 ) - 8 e^3 ( 108 - 153e^2 -9 e^4 +62e^6) \mathcal{L}_e\nonumber\\ & \quad + 4e^2 ( 162 - 297e^2 + 147e^4 - 23e^6 + 15 e^8) \mathcal{L}_e^2 \nonumber\\ & \quad - 2e (108 - 243e^2 + 189e^4 - 77 e^6 + 23 e^8) \mathcal{L}_e^3 + 3(3-4e^2+e^4)^2 \mathcal{L}_e^4 ]\nonumber\\ & \quad \times \{ 12 e^5(2e -\mathcal{L}_e) [6e + ({-}3+e^2)\mathcal{L}_e][2e+({-}1+e^2)\mathcal{L}_e] \}^{{-}1}, \end{align}
\begin{align} \mathcal{X}^{\varOmega',c} & = \frac{ (1 - e^2) [ 4 e^2 (5 - 4e^2) + 20e({-}1+e^2) \mathcal{L}_e + (5 - 6e + e^4) \mathcal{L}_e^2]}{(2 - e^2) [2e + ({-}1+e^2)\mathcal{L}_e]^2}. \end{align}Adding (C16) and (C17) to (5.3) and (5.4), we obtain the combined effects of the ambient and disturbance viscosities on the particle's translational and rotational velocities:
$$\begin{gather} \boldsymbol{U}_1 = \frac{\eta_c}{12 \varepsilon \eta_{\infty}}\,\mathcal{Z}^{U} \boldsymbol{U}_0 - \frac{a B_2}{6}\,(\mathcal{X}^{U,c} \boldsymbol{I} - \mathcal{Y}^{U,c} 3 \boldsymbol{pp}) \boldsymbol{\cdot} \boldsymbol{\nabla} \left(\frac{\eta_0 }{ \eta_{\infty}}\right), \end{gather}$$
$$\begin{gather}\boldsymbol{\varOmega}_1 ={-} \frac{1}{4}\,\mathcal{X}^{\varOmega,c} \boldsymbol{U}_{0} \times \boldsymbol{\nabla} \left(\frac{\eta_0 }{ \eta_{\infty}}\right), \end{gather}$$where
\begin{align} \mathcal{Z}^{U} & = \frac{6 (1 - e^2) (2 e - \mathcal{L}_e)^2}{e^2 \mathcal{L}_e [2e - (1 - e^2 ) \mathcal{L}_e]}, \end{align}
\begin{align} \mathcal{X}^{U,c} & = [{-}4e^2(45 - 75e^2 + 32 e^4 ) + 12 e (15 - 28e^2 + 13e^4) \mathcal{L}_e \nonumber\\ & \quad - 9 (1 - e^2)^2 (5 + e^2) \mathcal{L}_e^2 + 6e (1 - e^2)^2 \mathcal{L}_e^3] \nonumber\\ & \quad \times \{ 2 e^2 [6e + ({-}3+e^2)\mathcal{L}_e][2e + ({-}1+e^2)\mathcal{L}_e]\}^{{-}1},\nonumber\\ \mathcal{Y}^{U,c} & = [{-}4e^4({-}39 + 32e^2 ) - ( 168 e^3 - 156 e^5) \mathcal{L}_e \nonumber\\ & \quad - e^2 ({-}63 + 54e^2 + 5e^4 ) \mathcal{L}_e^2 - e (15 - 16e^2 + e^4) \mathcal{L}_e^3 + 3(1 - e^2)^2 \mathcal{L}_e^4 ]\nonumber\\ & \quad \times \{ 3 e [6e + ({-}3+e^2)\mathcal{L}_e](2e - \mathcal{L}_e)[2e+({-}1+e^2)\mathcal{L}_e] \}^{{-}1}, \end{align}
\begin{align} \mathcal{X}^{\varOmega,c} & = \frac{ (1 - e^2) [ 4 e^2 (- 7 + 4e^2) - 4e({-}7 + 5e^2) \mathcal{L}_e + ({-}7 + 6e^2 + e^4) \mathcal{L}_e^2]}{(2 - e^2) [2e - (1-e^2)\mathcal{L}_e]^2}. \end{align}Compared to the no-flux condition, the disturbance viscosity here introduces a more complex influence on the swimming dynamics of a spheroidal particle. However, the particles will still generally display viscophobic dynamics.
