1. Introduction
Swimming microorganisms are widely found in nature and are important in a diverse array of biological processes (Vogel Reference Vogel1996; Lauga Reference Lauga2016; Wadhwa & Berg Reference Wadhwa and Berg2022; Ishikawa & Pedley Reference Ishikawa and Pedley2023a,Reference Ishikawa and Pedleyb). An understanding of the dynamics of microswimmers and other active particles (biological or synthetic) in homogeneous Newtonian fluids is now reasonably well developed (Brennen & Winet Reference Brennen and Winet1977; Lauga & Powers Reference Lauga and Powers2009; Yeomans, Pushkin & Shum Reference Yeomans, Pushkin and Shum2014; Elgeti, Winkler & Gompper Reference Elgeti, Winkler and Gompper2015; Ishikawa Reference Ishikawa2024). However, natural environments can often be quite complex and inhomogeneous (Bechinger et al. Reference Bechinger, Di Leonardo, Löwen, Reichhardt, Volpe and Volpe2016; Martínez-Calvo, Trenado-Yuste & Datta Reference Martínez-Calvo, Trenado-Yuste and Datta2023).
Inhomogeneity of fluid environments can arise due to spatial variations of different physical quantities, including light (Jékely Reference Jékely2009), heat (Bahat et al. Reference Bahat, Tur-Kaspa, Gakamsky, Giojalas, Breitbart and Eisenbach2003), chemical concentration and nutrients (Berg & Brown Reference Berg and Brown1972; Berg Reference Berg2004), and can lead to directed motion, which is known as taxis. Natural bodies of water, such as lakes, ponds and oceans, often have gradients in temperature or salinity, and these inhomogeneities can also lead to a stratification of the mechanical properties of the fluid such as density and viscosity (Stocker Reference Stocker2012). In response to gradients in viscosity, bacteria such as Leptospira and Spiraplasma have been observed to perform positive viscotaxis by swimming up the 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), while Escherichia coli and Chlamydomonas reinhardtii perform negative viscotaxis by swimming down the gradients (Sherman, Timkina & Glagolev Reference Sherman, Timkina and Glagolev1982; Coppola & Kantsler Reference Coppola and Kantsler2021; Stehnach et al. Reference Stehnach, Waisbord, Walkama and Guasto2021). A gradient in the viscosity of the intestinal mucosal barrier is thought to control the spatial organization of intestinal microbiota (Swidsinski et al. Reference Swidsinski, Sydora, Doerffel, Loening-Baucke, Vaneechoutte, Lupicki, Scholze, Lochs and Dieleman2007).
While there is an obvious benefit associated with moving to an environment where the energetic penalty of motion is lower, recent work has shown that such directed motion can arise as an immediate consequence of interaction with an environment with inhomogeneous mechanical properties. Inhomogeneity can break rotational symmetry, and this leads naturally to reorientation in order to conserve angular momentum. In pioneering work, Liebchen et al. (Reference Liebchen, Monderkamp, ten Hagen and Löwen2018) showed that active particles, modelled as linked spheres moving with a constant propulsive force, in weak viscosity gradients exhibited positive viscotaxis due to the mismatch of viscous drag on the spheres. Further investigations incorporated the impact of viscosity changes on propulsion using the spherical (Datt & Elfring Reference Datt and Elfring2019; Shaik & Elfring Reference Shaik and Elfring2021; Gong, Shaik & Elfring Reference Gong, Shaik and Elfring2023) and spheroidal (Gong, Shaik & Elfring Reference Gong, Shaik and Elfring2024) squirmer models, which simulates a swimming gait by a surface slip velocity. These studies found that the interaction between the particle's active slip and spatial variations of viscosity tends to dominate the dynamics and typically results in negative viscotaxis, although the effect of the gradient decreases with increased slenderness.
Changes in translational and rotational dynamics that arise as a consequence of viscosity gradients are now relatively well understood, for both passive (Kamal & Lauga Reference Kamal and Lauga2023; Anand & Narsimhan Reference Anand and Narsimhan2024) and active (Gong et al. Reference Gong, Shaik and Elfring2024) particles, but here we focus on the power expended by microswimmers and any associated changes in efficiency that are caused by the inhomogeneous environment. Even though some microswimmers are capable of achieving higher speeds in viscosity gradients, we still need to quantify the resultant changes in power expenditure in order to determine whether such motion is in fact more efficient. We assess the efficiency of microswimmers using the Froude efficiency, which is defined as the ratio of the power expended to simply drag a body of the same shape to the power expended by propulsion (Lighthill Reference Lighthill1952). This efficiency has been used broadly to characterize different microswimmers in homogeneous Newtonian fluids (Stone & Samuel Reference Stone and Samuel1996; Chattopadhyay et al. Reference Chattopadhyay, Moldovan, Yeung and Wu2006; Ishimoto & Gaffney Reference Ishimoto and Gaffney2014), although more recently, other efficiency measures have been proposed (Childress Reference Childress2012; Nasouri, Vilfan & Golestanian Reference Nasouri, Vilfan and Golestanian2021; Daddi-Moussa-Ider, Golestanian & Vilfan Reference Daddi-Moussa-Ider, Golestanian and Vilfan2023), to guarantee a measure that can never exceed unity.
Several studies have investigated the efficiency of microswimmers in homogeneous non-Newtonian environments. One study showed that pusher-type spherical squirmers always expend more power and swim less efficiently than puller-type squirmers in second-order fluids (De Corato, Greco & Maffettone Reference De Corato, Greco and Maffettone2015). In fluids with shear-thinning properties, all types of two-mode spherical squirmers swim more slowly, expending less power with higher efficiency (Nganguia, Pietrzyk & Pak Reference Nganguia, Pietrzyk and Pak2017), but by introducing a third squirming mode, it is possible to design a microswimmer that can swim faster and more efficiently in shear-thinning fluids. Specifically for mechanically inhomogeneous fluids, a swimming sheet was shown to expend more power and move with reduced efficiency in a viscosity stratified fluid in comparison to a homogeneous Newtonian fluid (Dandekar & Ardekani Reference Dandekar and Ardekani2020).
In this work, we aim to explore how viscosity gradients affect the swimming efficiency of spheroidal squirmers in order to understand the role of swimming gait and shape on power expended and efficiency in an inhomogeneous fluid. In what follows, we derive analytical formulas for the energy expended by these microswimmers and the associated efficiency assuming weak constant gradients in the viscosity.
2. Swimmers in viscosity gradients
2.1. The squirmer model
Microswimmers are modelled as squirmers in this research. The spherical squirmer model is a classical hydrodynamic model for the motion of self-propelling particles, such as protozoa or volvocine green algae, wherein the complex motions of ciliated surfaces are represented as a tangential slip velocity on the surface of a spherical body (Lighthill Reference Lighthill1952; Blake Reference Blake1971). The slip velocity is commonly written in terms of Legendre polynomials:
\begin{equation} \boldsymbol{u}^s ={-} \sum_{n=1}^{\infty} \frac{2 B_n}{n(n+1)}\,P_n'(\kern0.7pt\boldsymbol{p} \boldsymbol{\cdot} \boldsymbol{n} )\,\boldsymbol{p} \boldsymbol{\cdot}(\boldsymbol{\mathsf{I}} - \boldsymbol{n} \boldsymbol{n} ), \end{equation}
where 
$\boldsymbol {p}$ is the swimming direction of the particle, 
$P_n$ is the Legendre polynomial of degree 
$n$, and 
$\boldsymbol {n}$ is a unit normal to the squirmer surface. The coefficients 
$B_n$, often called squirming modes, can be related to Stokes flow singularity solutions (Pak & Lauga Reference Pak and Lauga2014). The swimming velocity of a particle is determined by the 
$B_1$ mode (here, we assume 
$B_1 \geqslant 0$), while the 
$B_2$ mode sets the magnitude of the force dipole, the slowest decaying contribution to the far-field flow.
Since many organisms are elongated in shape, such as Paramecium caudatum or Opalina, the squirmer model was extended to accommodate spheroidal geometries by Keller & Wu (Reference Keller and Wu1977), and the streamlines predicted by their spheroidal model indeed closely matched the experimental streak photographs of Paramecium caudatum. The original model by Keller & Wu (Reference Keller and Wu1977) included only the first squirming mode, and latter studies added a force-dipole mode to better represent other types of microswimmers (Gaffney et al. Reference Gaffney, Gadêlha, Smith, Blake and Kirkman-Brown2011; Ishimoto & Gaffney Reference Ishimoto and Gaffney2014; Theers et al. Reference Theers, Westphal, Gompper and Winkler2016). We use the two-mode spheroidal squirmer here to represent prolate spheroidal active swimmers (Theers et al. Reference Theers, Westphal, Gompper and Winkler2016; van Gogh et al. Reference van Gogh, Demir, Palaniappan and Pak2022; Gong et al. Reference Gong, Shaik and Elfring2024). The slip velocity on the surface of a prolate spheroidal squirmer can be expressed as
\begin{equation} \boldsymbol{u}^s ={-}B_1 (\boldsymbol{s} \boldsymbol{\cdot} \boldsymbol{p})\,\boldsymbol{s} - \frac{B_2}{a}\,(\boldsymbol{r} \boldsymbol{\cdot} \boldsymbol{p} ) (\boldsymbol{s} \boldsymbol{\cdot} \boldsymbol{p}) \boldsymbol{s}, \end{equation}
where 
$\boldsymbol {r}$ is a vector from the centre of the particle to a point on the surface, while 
$\boldsymbol {s}$ is the unit tangent at that point (see figure 1). This form is equivalent to the formula in (2.1) for a sphere with only two modes. Recent research points out that the swimming speed and stresslet of such a squirmer are influenced by more than just the 
$B_1$ and 
$B_2$ modes (Pöhnl, Popescu & Uspal Reference Pöhnl, Popescu and Uspal2020); however, these additional modes significantly affect the outcome only when the particle is notably slender. For most practical cases, the two-mode prolate squirmer model (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; Gong et al. Reference Gong, Shaik and Elfring2024) suffices to capture swimming dynamics. Thus identifying the coefficients 
$B_1$ and 
$B_2$ completely determines the swimming mechanism here.
Figure 1. Schematic of a spheroidal squirmer swimming in a constant viscosity gradient 
$\boldsymbol {\nabla } \eta$. Here, 
$\boldsymbol {p}$ represents the swimming direction of the swimmer.
It is common to define the squirming ratio 
$\beta = B_2/B_1$, in which case the sign of 
$\beta$ determines whether the propulsion of the swimmer originates from its front or rear. Organisms like Escherichia coli, which generate propulsion from their rear, are classified as pushers and are characterized by 
$\beta < 0$. Conversely, organisms that draw fluid towards them from the front, such as Chlamydomonas reinhardtii, are known as pullers and are characterized by 
$\beta > 0$. There are also swimmers whose propulsion mechanisms are not predominantly directed from the front or the rear; one example is Volvox carteri, which has flagella evenly distributed across its spherical surface. For these organisms, referred to as neutral swimmers, 
$\beta = 0$.
2.2. Newtonian fluid with viscosity gradients
We endeavour to find the power expended by an active particle swimming in a Newtonian fluid with a spatially varying viscosity, and thereafter determine the efficiency of that motion. The active particle is modelled as a prolate spheroid. The spheroid is an axisymmetric ellipsoid with two equal-length equatorial semi-axes 
$b$ and one longer polar semi-axis 
$a$ (
$b \leqslant a$). The eccentricity
\begin{equation} e = \sqrt{1 - \frac{b^2}{a^2}} \end{equation}
represents the slenderness of the particle. When 
$e = 0$ (
$b = a$), the particle is a sphere with radius 
$a$. The swimming direction 
$\boldsymbol {p}$ of the particle is the direction of its major axis.
Viscosity variations in fluids are caused by similar variations in temperature, salinity, nutrient or polymer concentrations that often occur over length scales much larger than the swimming microorganisms. For instance, the viscosity variations in microfluidic channels occur over 
$10^2\unicode{x2013}10^3$ 
$\mathrm{\mu }$m, while those in oceans, lakes and ponds occur over several metres to kilometres. Generally, it is reasonable to assume that viscosity varies linearly at the scale of small (
$\mathrm {\mu }$m-sized) organisms considered here, and indeed such local to particle viscosity variations have been shown to largely capture changes in particle dynamics (Shaik & Elfring Reference Shaik and Elfring2021). Moreover, one can also achieve close to linear viscosity profiles in microfluidic channels (Stehnach et al. Reference Stehnach, Waisbord, Walkama and Guasto2021). Hence we consider a linearly varying viscosity field 
$\eta (\boldsymbol {x})$ varying over a large macroscopic length scale 
$L$ such that
\begin{equation} \boldsymbol{\nabla} \eta = \frac{\eta_{\infty}}{L}\,\boldsymbol{d}, \end{equation}
where 
$\eta _{\infty }$ is the viscosity at an arbitrary point near the particle, and 
$\boldsymbol {d}$ is the unit direction of the viscosity gradient. Defining the dimensionless parameter 
$\epsilon = a/L \ll 1$ that characterizes the size of the particle vs the length scale of the gradient, we can write 
$\boldsymbol {\nabla } \eta = \epsilon ({\eta _{\infty }}/{a}) \boldsymbol {d}$.
The flow around the particle in a Newtonian fluid with a spatially varying viscosity at low Reynolds number is governed by the continuity equation and Cauchy's equation of motion:
\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 } = - p \boldsymbol{\mathsf{I}} + \eta (\boldsymbol {x})\,\dot {\boldsymbol {\gamma }}$ is the stress tensor. Here, 
$p$ is the pressure, and 
$\dot {\boldsymbol {\gamma }} = \boldsymbol {\nabla } \boldsymbol {u} + (\boldsymbol {\nabla } \boldsymbol {u})^{\rm T}$ is twice the strain-rate tensor. Defining an extra stress
\begin{equation} \boldsymbol{\tau}_{NN} = (\eta (\boldsymbol{x}) - \eta_{\infty}) \dot{\boldsymbol{\gamma}}, \end{equation}we can alternatively write the stress tensor as
\begin{equation} \boldsymbol{\sigma} ={-} p \boldsymbol{\mathsf{I}} + \eta_{\infty} \dot{\boldsymbol{\gamma}} + \boldsymbol{\tau}_{NN}. \end{equation}The boundary conditions on the flow field are that far away from the particle, the disturbance flow generated by the microswimmer should decay, which means that
\begin{equation} \boldsymbol{u} \rightarrow \boldsymbol{0} \quad \text{as} \ |\boldsymbol{r}| \rightarrow \infty. \end{equation}
While on the particle surface (denoted by 
$S_p$), the fluid velocity 
$\boldsymbol {u}$ should satisfy no-slip conditions
\begin{equation} \boldsymbol{u} ( \boldsymbol{x} \in S_p) = \boldsymbol{U} + \boldsymbol{\varOmega} \times \boldsymbol{r} + \boldsymbol{u}^s. \end{equation}
The slip velocity 
$\boldsymbol {u}^s$ is prescribed by (2.2), while the translational and rotational velocities, 
$\boldsymbol {U}$ and 
$\boldsymbol {\varOmega }$, are determined by satisfying the additional constraints that since the microswimmers are considered as neutrally buoyant, there is no external force or torque acting on the particle, and in the absence of particle inertia, the hydrodynamic force and torque on the particle must both vanish:
\begin{gather} \boldsymbol{F} = \int_{S_p} \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\sigma} \, \text{d} S = \boldsymbol{0}, \end{gather}
\begin{gather}\boldsymbol{L} = \int_{S_p} \boldsymbol{r} \times (\boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\sigma}) \, \text{d} S = \boldsymbol{0}. \end{gather} The power 
$\mathcal {P}$ expended by swimmer in its motion through the fluid is written as
\begin{equation} \mathcal{P} ={-}\int_{S_p} \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\sigma} \boldsymbol{\cdot} \boldsymbol{u} \, \text{d} S ={-}\int_{S_p} \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\sigma} \boldsymbol{\cdot} \boldsymbol{u}^s \, \text{d} S, \end{equation}
where the latter form comes from substitution of the boundary conditions (2.10) and application of the force- and torque-free conditions. Alternatively, an application of the divergence theorem yields the power expended by the swimmer in terms of the viscous dissipation in the fluid volume 
$\mathcal {V}$ in which the swimmer is immersed:
\begin{equation} \mathcal{P} = \int_{\mathcal{V}} \boldsymbol{\sigma}:\boldsymbol{\nabla} \boldsymbol{u} \, \text{d} V = \frac{1}{2}\int_{\mathcal{V}} \eta(\boldsymbol{x})\, \dot{\boldsymbol{\gamma}}:\dot{\boldsymbol{\gamma}} \, \text{d} V. \end{equation} Lighthill (Reference Lighthill1952) defined the efficiency 
$\mathcal {E}$ of a particle swimming at low Reynolds number as the ratio of the power 
$\hat {\mathcal {P}}$ required to simply drag the particle (as a rigid body) at the same velocity as it swims, to the power expended by swimming 
$\mathcal {P}$, in other words,
\begin{equation} \mathcal{E} = \frac{\hat{\mathcal{P}}}{\mathcal{P}}. \end{equation} The towing power 
$\hat {\mathcal {P}}$ in (2.15) is given by
\begin{equation} \hat{\mathcal{P}} ={-}\hat{\boldsymbol{F}} \boldsymbol{\cdot} \boldsymbol{U} - \hat{\boldsymbol{L}} \boldsymbol{\cdot} \boldsymbol{\varOmega}, \end{equation}
where the hydrodynamic force 
$\hat {\boldsymbol {F}}$ and torque 
$\hat {\boldsymbol {L}}$ due to undergoing the rigid-body motion with velocities 
$\boldsymbol {U}$ and 
$\boldsymbol {\varOmega }$ can be written as
\begin{equation} \begin{pmatrix} \hat{\boldsymbol{F}} \\ \hat{\boldsymbol{L}} \end{pmatrix}={-}\begin{pmatrix} \boldsymbol{\mathsf{R}}_{\boldsymbol{F U}} & \boldsymbol{\mathsf{R}}_{\boldsymbol{F \varOmega}} \\ \boldsymbol{\mathsf{R}}_{\boldsymbol{L U}} & \boldsymbol{\mathsf{R}}_{\boldsymbol{L \varOmega}} \end{pmatrix}\boldsymbol{\cdot} \begin{pmatrix} \boldsymbol{U} \\ \boldsymbol{\varOmega} \end{pmatrix}, \end{equation}
where 
$(\boldsymbol{\mathsf{R}}_{\boldsymbol {F U}}, \boldsymbol{\mathsf{R}}_{\boldsymbol {F \varOmega }}, \boldsymbol{\mathsf{R}}_{\boldsymbol {L U}}, \boldsymbol{\mathsf{R}}_{\boldsymbol {L \varOmega }})$ are rigid-body resistance tensors dependent on the eccentricity 
$e$, the orientation vector 
$\boldsymbol {p}$, and the viscosity gradient 
$\boldsymbol {\nabla } \eta$. The detailed expressions are given in Appendix B.
3. Asymptotic analysis
Since the gradient in viscosity is considered weak, as characterized by 
$\epsilon \ll 1$, we expand all flow terms 
$\{ \boldsymbol {u}, \boldsymbol {\sigma }, \boldsymbol {\tau }_{NN}, \dot {\boldsymbol {\gamma }}, \mathcal {P}, \hat {\mathcal {P}}, \mathcal {E} \}$ assuming a regular perturbation expansion in 
$\epsilon$, for example, 
$\boldsymbol {u} = \boldsymbol {u}_0 + \epsilon \boldsymbol {u}_1 + O(\epsilon ^2)$.
The spatially varying viscosity field can be written in terms of the centre of the particle 
$\boldsymbol {x}_c$ as
\begin{equation} \eta(\boldsymbol{x}) = \eta(\boldsymbol{x}_c) +(\boldsymbol{x}-\boldsymbol{x}_c) \boldsymbol{\cdot} \boldsymbol{\nabla} \eta = \eta(\boldsymbol{x}_c) +\epsilon\, \frac{\eta_{\infty}}{a}\,\boldsymbol{d} \boldsymbol{\cdot} (\boldsymbol{x} - \boldsymbol{x}_c). \end{equation}
The centre of the particle need not necessarily be in the plane 
$\mathcal {X}_\infty$ where 
$\eta (\boldsymbol {x}) = \eta _\infty$, but we can always define 
$\eta _{\infty }$ near the particle such that the distance 
$| \boldsymbol {d} \boldsymbol {\cdot } (\boldsymbol {x} -\boldsymbol {x}_c) |/a \ll O(1/\epsilon )$ for 
$\boldsymbol {x} \in \mathcal {X}_{\infty }$, hence 
$\eta (\boldsymbol {x}_c) = \eta _\infty + O(\epsilon )$.
3.1. Homogeneous fluids
At leading order, 
$\epsilon \rightarrow 0$, the viscosity is homogeneous, so we simply have an active particle swimming in a Newtonian fluid with a constant viscosity. The corresponding flow field is well known, for both spherical and spheroidal active particles, and given in Appendix A. According to (2.14), the leading-order value of power dissipation has the form
\begin{equation} \mathcal{P}_0 = \frac{1}{2}\int_{\mathcal{V}} \eta_\infty \dot{\boldsymbol{\gamma}}_0:\dot{\boldsymbol{\gamma}}_0 \,\text{d} S. \end{equation}Evaluating the integral, we obtain
\begin{equation} \mathcal{P}_0 = {\rm \pi}a \eta_\infty (\mathcal{A} + \mathcal{B} \beta^2)B_1^2, \end{equation}
where the terms 
$\mathcal {A}$ and 
$\mathcal {B}$ are functions of only the eccentricity and have the form
\begin{gather} \mathcal{A} = \frac{2(1-e^2)[{-}2e + (1 + e^2) \mathcal{L}_e] }{e^3}, \end{gather}
\begin{gather}\mathcal{B} \,{=}\,
\frac{-8e^2(45\,{-}\, 51e^2\,{+}\, 8e^4)\,{+}\, 24e(15 \,{-}\,
22e^2 \,{+}\, 7e^4)\mathcal{L}_e\,{+}\, 6({-}15\,{+}\, 27e^2\,{-}\, 13e^4\,{+}\, e^6)\mathcal{L}_e^2}{3e^5[6e + ({-}3 +
e^2)\mathcal{L}_e]}, \end{gather}
\begin{gather}\mathcal{L}_e = \ln ((1+e)/(1-e)). \end{gather} Keller & Wu (Reference Keller and Wu1977) derived the above result for a neutral spheroidal squirmer, but here we have added the contribution of the second mode to obtain a formula valid for pushers and pullers. When 
$e \rightarrow 0$, 
$\mathcal {A}= 16/3$ and 
$\mathcal {B} = 8/3$, and we obtain the result for a spherical squirmer derived by Lighthill (Reference Lighthill1952) and Blake (Reference Blake1971):
\begin{equation} \mathcal{P}_{0,{sphere}} = \tfrac{8}{3} {\rm \pi}(2 + \beta ^2) \eta_{\infty} a B_1^2 = 6 {\rm \pi}a \eta_\infty (2 + \beta ^2) U_{0,{sphere}}^2, \end{equation}
where 
$U_{0,{sphere}} = 2B_1/3$ is the swimming speed of a spherical squirmer in Newtonian fluid with uniform viscosity. Here, 
$\mathcal {A}$ and 
$\mathcal {B}$ are monotonically decreasing functions of eccentricity, meaning that power diminishes as a spheroidal squirmer becomes more slender, and when 
$e \rightarrow 1$, 
$\mathcal {A} =\mathcal {B} = 0$, meaning that an infinitely slender squirmer requires no power to move (with speed equal to 
$B_1$) through a homogeneous Newtonian fluid.
The efficiency at leading order is
\begin{equation} \mathcal{E}_0 = \frac{\hat{\mathcal{P}}_0}{\mathcal{P}_0}= \frac{\boldsymbol{U}_0 \boldsymbol{\cdot} \boldsymbol{\mathsf{R}}_{\boldsymbol{F U}}^{(0)} \boldsymbol{\cdot} \boldsymbol{U}_0 + \boldsymbol{\varOmega}_0 \boldsymbol{\cdot} \boldsymbol{\mathsf{R}}_{\boldsymbol{L \varOmega}}^{(0)} \boldsymbol{\cdot} \boldsymbol{\varOmega}_0 }{\mathcal{P}_0}. \end{equation}Using classic results of the rigid-body resistance (see Appendix B) and the velocity of a spheroidal squirmer in a homogeneous Newtonian fluid (Keller & Wu Reference Keller and Wu1977),
\begin{gather} \boldsymbol{U}_0 = \frac{2 e - (1 - e^2) \mathcal{L}_e}{2 e^3}\,B_1\boldsymbol{p}, \end{gather}
\begin{gather}\boldsymbol{\varOmega}_0 = \boldsymbol{0}, \end{gather}we obtain the efficiency
\begin{equation} \mathcal{E}_0 = \frac{\mathcal{F}}{\mathcal{A} + \mathcal{B} \beta^2}, \end{equation}where
\begin{equation} \mathcal{F} = \frac{4[2e + ({-}1 + e^2) \mathcal{L}_e]^2}{e^3[{-}2e + (1 + e^2) \mathcal{L}_e]} \end{equation}
is a monotonically decreasing function of slenderness 
$e$. When 
$e\rightarrow 0$, 
$\mathcal {F} = 8/3$ and we obtain the classic efficiency of a two-mode spherical squirmer,
\begin{equation} \mathcal{E}_{0,{sphere}}=\frac{1}{2 + \beta^2}. \end{equation}
In the slender limit, 
$\mathcal {E}_0(e\rightarrow 1) = 3/\beta ^2$, meaning that efficiency is finite given the presence of the wasteful 
$B_2$ mode.
3.2. Effect of gradients
The effect of viscosity variations on the power dissipation in (2.14) are captured at 
$O(\epsilon )$ by
\begin{equation} \epsilon \, \mathcal{P}_1 = \frac{1}{2} \int_{\mathcal{V}} ( \eta(\boldsymbol{x}) - \eta_{\infty}) \dot{\boldsymbol{\gamma}}_0 : \dot{\boldsymbol{\gamma}}_0 \, \text{d} V. \end{equation}
Terms involving 
$\dot {\boldsymbol {\gamma }}_1$ are identically zero due to the force- and torque-free condition (De Corato et al. Reference De Corato, Greco and Maffettone2015; Nganguia et al. Reference Nganguia, Pietrzyk and Pak2017), so power at this order involves only integrating the leading-order homogeneous flow field. Similarly, the particle dynamics in viscosity gradients can be written as an integral of the leading-order flow field, although against a different kernel (Gong et al. Reference Gong, Shaik and Elfring2024), hence as the efficiency depends on power and particle dynamics, it can also be found from the leading-order homogeneous flow field. With the above formula, we obtain the power dissipation, valid up to 
$O(\epsilon )$:
\begin{equation} \mathcal{P} ={\rm \pi} (\mathcal{A} + \mathcal{B} \beta^2)\,\eta(\boldsymbol{x}_c)\,a B_1^2 + {\rm \pi}a^2 \mathcal{C} B_1^2 \beta( \boldsymbol{p} \boldsymbol{\cdot} \boldsymbol{\nabla} \eta), \end{equation}where
\begin{equation} \mathcal{C} = \frac{4({-}1 + e^2)[6e + ({-}3 + e^2) \mathcal{L}_e]}{e^5} \end{equation}
is a monotonically decreasing function of eccentricity. Taking the limit 
$e \rightarrow 0$, 
$\mathcal {C}= 32/15$ and we obtain for a sphere that
\begin{equation} \mathcal{P}_{{sphere}} = \frac{8}{3}\,{\rm \pi} (2 + \beta^2)\,\eta(\boldsymbol{x}_c)\,a B_1^2 + \frac{32}{15}\,{\rm \pi} a^2 B_1^2\beta( \boldsymbol{p} \boldsymbol{\cdot} \boldsymbol{\nabla}\eta). \end{equation} What we see is that pushers (
$\beta <0$) require less power to swim up viscosity gradients, and more to swim down viscosity gradients, in comparison to a homogeneous fluid, while the opposite is true for pullers (
$\beta >0$). Neutral squirmers (
$\beta =0$) do not see a change in power. The reason for this change in power expended can be understood in terms of the change in speed experienced by the squirmers. Previous research (Datt & Elfring Reference Datt and Elfring2019; Gong et al. Reference Gong, Shaik and Elfring2024) showed that pushers generate more thrust and thus swim faster and dissipate more energy when moving down viscosity gradients; conversely, they generate less thrust, swim slower and dissipate less energy when moving up viscosity gradients. The opposite is true for pullers.
A spheroidal squirmer always dissipates less power (for the same 
$B_1$ and 
$B_2$ values) as 
$\mathcal {A}$, 
$\mathcal {B}$ and 
$\mathcal {C}$ are all monotonically decreasing functions of the slenderness 
$e$, that all go to zero when 
$e \rightarrow 1$.
The correction to the efficiency, at 
$O(\epsilon )$, can be written as
\begin{equation} \mathcal{E}_1 = \frac{\hat{\mathcal{P}}_1}{\mathcal{P}_0}-\mathcal{E}_0\,\frac{\mathcal{P}_1}{\mathcal{P}_0}. \end{equation}
From Appendix B, we know that 
$\boldsymbol{\mathsf{R}}_{\boldsymbol {F U}}$, 
$\boldsymbol{\mathsf{R}}_{\boldsymbol {L \varOmega }}$ are 
$O(1)$ at the leading order, while 
$\boldsymbol{\mathsf{R}}_{\boldsymbol {F \varOmega }}, \boldsymbol{\mathsf{R}}_{\boldsymbol {L U}}$ are 
$O(\epsilon )$ as a result:
\begin{equation} \hat{\mathcal{P}}_1 = 2 \boldsymbol{U}_1 \boldsymbol{\cdot} \boldsymbol{\mathsf{R}}_{\boldsymbol{F U}}^{(0)} \boldsymbol{\cdot} \boldsymbol{U}_0 + \boldsymbol{U}_0 \boldsymbol{\cdot} \boldsymbol{\mathsf{R}}_{\boldsymbol{F U}}^{(1)} \boldsymbol{\cdot} \boldsymbol{U}_0. \end{equation}
The correction to particle velocities due to viscosity variations, 
$\boldsymbol {U}_1$ and 
$\boldsymbol {\varOmega }_1$, are from Gong et al. (Reference Gong, Shaik and Elfring2024):
\begin{gather} \boldsymbol{U}_1 ={-}\frac{B_2}{5}\,( \mathcal{X}^U \boldsymbol{I} - \mathcal{Y}^U 3 \boldsymbol{pp} ) \boldsymbol{\cdot} \boldsymbol{d}, \end{gather}
\begin{gather}\boldsymbol{\varOmega}_1 ={-}\frac{1}{2a}\,\mathcal{X}^{\varOmega} \boldsymbol{U}_0 \times \boldsymbol{d}, \end{gather}where
\begin{gather} \mathcal{X}^U = \frac{5 [{-}6e + 4e^3 + 3( 1 - e^2 ) \mathcal{L}_e ] [{-}6e + 10e^3 + 3 ( 1- e^2 )^2 \mathcal{L}_e ]}{24e^5 [ 6e - ( 3 - e^2) \mathcal{L}_e ]}, \end{gather}
\begin{gather}\mathcal{Y}^U = \frac{5 [{-}6e + 4e^3 + 3 ( 1 - e^2 ) \mathcal{L}_e ] [{-}18e + 6e^3 + ( 9 - 6e^2 + 5e^4 ) \mathcal{L}_e ]}{72e^5 [ 6e - ( 3 - e^2 ) ]}, \end{gather}
\begin{gather}\mathcal{X}^{\varOmega} = \frac{( 1 - e^2 ) [{-}2e + ( 1 + e^2 ) \mathcal{L}_e ]}{(2 - e^2 ) [2e - ( 1 - e^2 ) \mathcal{L}_e ] }. \end{gather}
However, note that 
$\boldsymbol {\varOmega }_1$ is not needed for 
$\hat {\mathcal {P}}_1$. Substituting 
$\hat {\mathcal {P}}_1$ and 
$\mathcal {P}_1$ from (3.15) into (3.18), we can calculate the efficiency valid to 
$O(\epsilon )$ as
\begin{equation} \mathcal{E} = \frac{\mathcal{F}}{\mathcal{A} + \mathcal{B} \beta^2} + \epsilon\, \frac{\mathcal{G} + \mathcal{H} \beta^2}{(\mathcal{A} + \mathcal{B} \beta^2)^2}\, \beta(\kern0.7pt\boldsymbol{p} \boldsymbol{\cdot} \boldsymbol{d}),\end{equation}
where 
$\mathcal {G}$ (for 
$e < 0.964$) and 
$\mathcal {H}$ are monotonically decreasing function of slenderness:
\begin{align} \mathcal{G} & = 8(1\,{-}\,
e^2)[2e\,{+}\, ({-}1\,{+}\, e^2) \mathcal{L}_e][8e^3(45\,{+}\,
3e^2\,{+}\, 2e^4)\,{-}\, 4e^2(135\,{-}\, 84e^2\,{+}\, 5e^4\,{+}\, 4e^6)\mathcal{L}_e \nonumber\\ & \quad + 2e(135 -177e^2
+ 51e^4 + 5e^6 + 2e^8)\mathcal{L}_e^2 -
3(1-e^2)^2(15+e^4)\mathcal{L}_e^3] \nonumber\\ & \quad
\times \{3e^8 [6e + ({-}3 + e^2) \mathcal{L}_e][{-}2e + (1
+ e^2) \mathcal{L}_e] \}^{{-}1}, \end{align}
\begin{align} \mathcal{H} & = 4[2e +
({-}1 + e^2) \mathcal{L}_e] \nonumber\\ & \quad
\times [8e^2(45 \,{-}\, 51e^2\,{+}\, 8e^4)\,{-}\,
24e(15\,{-}\, 22e^2\,{+}\, 7e^4)\mathcal{L}_e\nonumber\\ &\quad -6({-}15 \,{+}\, 27e^2\,{-}\, 13e^4 \,{+}\, e^6)\mathcal{L}_e^2] \nonumber\\ &
\quad \times [4e^2({-}9 \,{+}\, 3e^2\,{+}\,
2e^4) \,{-}\, 4e({-}9 + 6e^2 + e^6)\mathcal{L}_e + 3({-}3
+3e^2\,{-}\, e^4 + e^6)\mathcal{L}_e^2] \nonumber\\& \quad
\times \{9e^{10} [6e + ({-}3 + e^2) \mathcal{L}_e]^2 [{-}2e
+ (1 + e^2) \mathcal{L}_e] \}^{{-}1}.
\end{align}
For spheres, 
$e\rightarrow 0$, the efficiency simplifies to
\begin{equation} \mathcal{E}_{{sphere}} = \frac{1}{2 + \beta^2} + \epsilon\,\frac{8 + 6 \beta^2}{5(2+ \beta^2)^2}\,\beta(\kern0.7pt\boldsymbol{p} \boldsymbol{\cdot} \boldsymbol{d}).\end{equation} The essential message of these formulas is that, similar to speed and power, the correction to efficiency due to viscosity gradients is proportional to 
$\beta (\boldsymbol {p}\boldsymbol {\cdot }\boldsymbol {d})$, meaning that, for example, when swimming down viscosity gradients, pushers are faster and expend more power, but are more efficient, while pullers are slower, expend less power and are less efficient. These trends are illustrated in figure 2, which shows the variation of relative power and efficiency as a function of 
$\beta$ for various eccentricities 
$e$. The opposite is true when swimming up viscosity gradients; however, we know that swimming up viscosity gradients is dynamically unstable, and all squirmers tend to perform negative viscotaxis.
Figure 2. Variation of (a) relative power expenditure and (b) swimming efficiency of a swimmer moving down viscosity gradients with squirming ratio 
$\beta$ and eccentricity 
$e$. Power expenditure is evaluated by choosing the viscosity at the centre of the swimmer as 
$\eta _{\infty }$ (
$\eta ( \boldsymbol {x}_c ) = \eta _{\infty }$).
In deriving these formulas, we have assumed that the background viscosity remains fixed even in the presence of the particle. However, variations in viscosity generally arise from variations in an underlying field that affects the viscosity, such as temperature, salt or nutrient concentration, and taking into account the effect of boundary conditions on the surface of the particle in relation to the underlying field will then lead to changes in the viscosity. In Appendix C, we determine the effects of imposing a ‘no-flux’ condition for viscosity at the surface of the particle, and the associated effects on power expended and efficiency. Although there are quantitative differences in the parameters, the qualitative picture described above remains unchanged.
We have also neglected the higher-order (
$n>2$) squirming modes. The higher-order modes do not affect the velocity of a spherical squirmer in a homogeneous fluid (Lighthill Reference Lighthill1952; Blake Reference Blake1971) or in viscosity gradients (Datt & Elfring Reference Datt and Elfring2019; Shaik & Elfring Reference Shaik and Elfring2021). But they do affect the velocity of a spheroidal squirmer in homogeneous fluids (Pöhnl et al. Reference Pöhnl, Popescu and Uspal2020), and we expect the effect to persist even in viscosity gradients. However, higher-order modes affect power dissipation and swimming efficiency in all cases, regardless of the presence of viscosity gradients or swimmer shape. In Appendix D, we give general formulas for power dissipation and efficiency for arbitrary spherical squirmers, and from those equations, we see that the effect of higher-order modes generally diminishes rapidly with the mode number 
$n$, and so may often be reasonably neglected. However, the presence of additional modes can substantially change the picture depending on the magnitude of the higher modes. For example, a three-mode pusher swimming down the viscosity gradients can dissipate more or less power and swim more or less efficiently as compared to in homogeneous fluids. We expect these broad insights on the effect of higher-order modes for a spherical squirmer to translate qualitatively to a spheroidal squirmer.
4. Conclusion
In this work, we derived analytical formulas for the power expended by spheroidal squirmers swimming in linearly varying viscosity fields, and the associated efficiency of that motion. We found that pushers are faster and more efficient when moving down gradients, but slower and less efficient moving up viscosity gradients, and the opposite is true for pullers. We also evaluated the effect of shape on power expenditure and efficiency when swimming in viscosity gradients, and found that the change in both due to gradients diminishes monotonically with increasing slenderness.
Swimming down viscosity gradients is favourable for pushers because they generate thrust from the rear, where the viscosity is highest, which leads to faster swimming and thus a greater level of power expenditure; however, the energetic cost of this boost is lower (relative to rigid body motion) in comparison to a homogeneous fluid, thus the increase in efficiency. While both pushers and pullers can minimize or maximize efficiency depending on orientation relative to the gradient, both display negative viscotaxis, which means that pushers dynamically tend to the most efficient orientation while pullers tend to the least. However, this conclusion is sensitive to geometry. For example, if we construct a pusher made from a very thin ‘tail’ that generates thrust, and a large head that bears the majority of the drag, then this sort of swimmer would display positive viscotaxis (Gong et al. Reference Gong, Shaik and Elfring2024), and be efficiency minimizing. This shows how shape and gait play a critical role in driving dynamics in inhomogeneous environments, and how both efficiency minimizing and maximizing dynamical states can be stable depending on the specific shape and gait. This raises questions about whether organisms might adjust their swimming gait in light of this fact. Finally, we note that while a puller is less efficient moving down a viscosity gradient in comparison to a homogeneous fluid of the same viscosity, it will still be energetically preferable to move to the region of lower viscosity, even if the instantaneous dynamics is slower and less efficient.
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. An active prolate spheroid in Stokes flow
The flow field 
$\boldsymbol {u}_0$ of an active spheroid swimming in a Newtonian fluid with uniform viscosity can be expressed 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) as
\begin{equation} \boldsymbol{u}_0 = \frac{1}{c^2 \sqrt{\zeta_1^2 - \zeta_2^2}} \left(\frac{1}{\sqrt{\zeta_1^2 - 1}}\,\frac{\partial \psi_0}{\partial \zeta_2}\,\boldsymbol{e}_{\zeta_1} - \frac{1}{ \sqrt{1 - \zeta_2^2 }}\,\frac{\partial \psi_0}{\partial \zeta_1}\,\boldsymbol{e}_{\zeta_2}\right), \end{equation}
using prolate spheroid coordinates 
$O'\zeta _1 \zeta _2 \phi$, where 
$O'$ is the centre of the spheroid at 
$\boldsymbol {x}_c$. The spheroidal coordinate system can be related to a particle-aligned Cartesian coordinate system 
$O'X_1X_2X_3$ as follows:
\begin{align} \left.\begin{gathered} X_1 = c \sqrt{\zeta_1^2 - 1} \sqrt{ 1 - \zeta_2^2} \cos \phi, \\ X_2 = c \sqrt{\zeta_1^2 - 1} \sqrt{ 1 - \zeta_2^2} \sin \phi, \\ X_3 = c \zeta_1 \zeta_2, \end{gathered}\right\} \end{align}
where 
$1 \leqslant \zeta _1 < \infty$, 
$-1 \leqslant \zeta _2 \leqslant 1$ and 
$0 \leqslant \phi < 2{\rm \pi}$. Here, 
$c = \sqrt {a^2 - b^2}$ is half of the focal length. We can define 
$\widetilde {\zeta _1} = 1 / e$, with 
$\zeta _1 > \tilde {\zeta }_1$ corresponding to the fluid domain exterior to the surface 
$\zeta _1 = \tilde {\zeta }_1$ of the particle.
The stream function 
$\psi _0$ has the solution
\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{gathered} C_1 = 2c^2\,\frac{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\,\frac{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\,\frac{ 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\,\frac{ 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{gathered}\right\} \end{equation}
and 
$U_0 = |\boldsymbol {U}_0|$ is given in (3.9).
Appendix B. Resistance tensor for a passive prolate spheroid in viscosity gradient
Here, we derive the resistance tensor for a spheroidal particle undergoing rigid-body motion in a fluid with a linearly varying viscosity field. We follow steps similar to those shown in a previous paper deriving the mobility of a spheroidal particle in a viscosity gradient (Gong et al. Reference Gong, Shaik and Elfring2024), and likewise, we use compact six-dimensional vectors for velocities 
$\boldsymbol{\mathsf{U}} = (\boldsymbol {U}, \boldsymbol {\varOmega })^\textrm {T}$ and forces and torques 
$\hat {\boldsymbol{\mathsf{F}}} = ( \hat {\boldsymbol {F}}, \hat {\boldsymbol {L}})^\textrm {T}$ to simplify formulas. Following the notation of this paper, the hat notation corresponds to values associated with rigid-body motion.
The resistance tensor to leading order, 
$\boldsymbol{\mathsf{R}}^{(0)}_{\boldsymbol{\mathsf{F}} \boldsymbol{\mathsf{U}}}$, corresponding to a homogeneous fluid with uniform viscosity 
$\eta _\infty$, is well known and given by Kim & Karilla (Reference Kim and Karilla1991). The hydrodynamic force and torque in a viscosity gradient can be written as
\begin{equation} \hat{\boldsymbol{\mathsf{F}}} ={-}\boldsymbol{\mathsf{R}}^{(0)}_{\boldsymbol{\mathsf{F}} \boldsymbol{\mathsf{U}}} \boldsymbol{\cdot} \boldsymbol{\mathsf{U}} + \hat{\boldsymbol{\mathsf{F}}}_{NN}. \end{equation}
The ‘extra’ hydrodynamic force and torque corresponding to the extra stress 
$\hat {\boldsymbol {\tau }}_{NN}$ from the gradient can be shown, by the reciprocal theorem, to be
\begin{equation} \hat{\boldsymbol{\mathsf{F}}}_{NN} ={-} \int_{\mathcal{V}} \hat{\boldsymbol{\tau}}_{NN} : \hat{\boldsymbol{\mathsf{E}}}_{\boldsymbol{\mathsf{U}}}^{(0)} \, \text{d} V, \end{equation}
where the operator 
$\hat {\boldsymbol{\mathsf{E}}}_{\boldsymbol{\mathsf{U}}}^{(0)}$ corresponds to the rate of strain tensor in a homogeneous fluid,
\begin{equation} \hat{\dot{\boldsymbol{\gamma}}}_0 = 2 \hat{\boldsymbol{\mathsf{E}}}_{\boldsymbol{\mathsf{U}}}^{(0)} \boldsymbol{\cdot} \boldsymbol{\mathsf{U}}, \end{equation}
and may be calculated with the flow field due to a spheroid undergoing rigid-body motion in a homogeneous Newtonian fluid, given by Kim & Karilla (Reference Kim and Karilla1991). Substituting 
$\hat {\boldsymbol {\tau }}_{NN} = (\eta (\boldsymbol {x}) - \eta _{\infty }) \widehat {\dot {\boldsymbol {\gamma }}}$ into the equation above, and expanding terms in powers of 
$\epsilon$, we obtain
\begin{equation} \hat{\boldsymbol{\mathsf{F}}}_{NN} ={-}\int_{\mathcal{V}} 2 (\eta (\boldsymbol{x}) - \eta_{\infty}) \hat{\boldsymbol{\mathsf{E}}}_{\boldsymbol{\mathsf{U}}}^{(0)} : \hat{\boldsymbol{\mathsf{E}}}_{\boldsymbol{\mathsf{U}}}^{(0)}\boldsymbol{\cdot} \boldsymbol{\mathsf{U}}\, \text{d} V + O (\epsilon^2), \end{equation}hence we can identify
\begin{equation} \epsilon \boldsymbol{\mathsf{R}}^{(1)}_{\boldsymbol{\mathsf{F}} \boldsymbol{\mathsf{U}}} = \int_{\mathcal{V}} 2 (\eta (\boldsymbol{x}) - \eta_{\infty}) \hat{\boldsymbol{\mathsf{E}}}_{\boldsymbol{\mathsf{U}}}^{(0)} : \hat{\boldsymbol{\mathsf{E}}}_{\boldsymbol{\mathsf{U}}}^{(0)}\, \text{d} V. \end{equation} Combining terms, the resistance tensor accurate to 
$O(\epsilon )$ is 
$\boldsymbol{\mathsf{R}}_{\boldsymbol{\mathsf{FU}}} = \boldsymbol{\mathsf{R}}^{(0)}_{\boldsymbol{\mathsf{F}} \boldsymbol{\mathsf{U}}} + \epsilon \boldsymbol{\mathsf{R}}^{(1)}_{\boldsymbol{\mathsf{F}} \boldsymbol{\mathsf{U}}}$, where
\begin{equation} \boldsymbol{\boldsymbol{\mathsf{R}}}_{\boldsymbol{\mathsf{FU}}} = \begin{pmatrix} \boldsymbol{\mathsf{R}}_{\boldsymbol{F} \boldsymbol{U}} \quad \boldsymbol{\mathsf{R}}_{\boldsymbol{F} \boldsymbol{\varOmega}} \\ \boldsymbol{\mathsf{R}}_{\boldsymbol{L} \boldsymbol{U}} \quad \boldsymbol{\mathsf{R}}_{\boldsymbol{L} \boldsymbol{\varOmega}} \end{pmatrix}, \end{equation}
and 
$\boldsymbol{\mathsf{R}}_{\boldsymbol {F} \boldsymbol {\varOmega }} = \boldsymbol{\mathsf{R}}_{\boldsymbol {L} \boldsymbol {U}}^\textrm {T}$. The resistance depends on the eccentricity 
$e$, the orientation vector 
$\boldsymbol {p}$, and the viscosity gradient 
$\boldsymbol {\nabla } \eta$. Detailed expressions for all couplings in (B6) are
\begin{gather} \boldsymbol{\mathsf{R}}_{\boldsymbol{FU}} = 6 {\rm \pi}\,\eta (\boldsymbol{x}_c)\,a [ \mathcal{X}^{A} \boldsymbol{p} \boldsymbol{p} + \mathcal{Y}^{A} (\boldsymbol{\mathsf{I}} - \boldsymbol{p} \boldsymbol{p})], \end{gather}
\begin{gather}\boldsymbol{\mathsf{R}}_{\boldsymbol{L\varOmega}} = 8 {\rm \pi}\,\eta (\boldsymbol{x}_c)\,a^3 [ \mathcal{X}^{C} \boldsymbol{p} \boldsymbol{p} + \mathcal{Y}^{C} (\boldsymbol{\mathsf{I}} - \boldsymbol{p} \boldsymbol{p})], \end{gather}
\begin{gather}\boldsymbol{\mathsf{R}}_{\boldsymbol{F\varOmega}} = \boldsymbol{\mathsf{R}}_{\boldsymbol{L U}}^{{\top}} = \epsilon 6 {\rm \pi}\eta_{\infty} a^2 [\lambda_1 ( \boldsymbol{d} \times \boldsymbol{\mathsf{I}}) + \lambda_2 ( \boldsymbol{p} \boldsymbol{\cdot} \boldsymbol{d}) ( \boldsymbol{p} \times \boldsymbol{\mathsf{I}}) + \lambda_3 \boldsymbol{p} (\boldsymbol{d} \times \boldsymbol{p})], \end{gather}where
\begin{gather} \lambda_1 = \frac{16 e^3 (1-e^2)}{9[ - 2e + (1 - 3 e^2) \mathcal{L}_e]}, \end{gather}
\begin{gather}\lambda_2 = \frac{8[ 6e^4(1 - e^2) + e^3 ({-}3 + 4e^2 + 3e^4) \mathcal{L}_e]}{9 [{-}2e + (1 + e^2) \mathcal{L}_e][ - 2e + (1 - 3 e^2) \mathcal{L}_e]}, \end{gather}
\begin{gather}\lambda_3 = \frac{4[e^5(36 - 28e^2) + e^4({-}36 + 40e^2 + 4e^4)\mathcal{L}_e + e^3(9 - 13e^2 - e^4 + 5e^6)\mathcal{L}_e^2]}{9 [{-}2e + (1 + e^2) \mathcal{L}_e]^2[ - 2e + (1 - 3 e^2) \mathcal{L}_e]}, \end{gather}
\begin{gather}\mathcal{X}^A =\frac{8 e^3}{3[ - 2e + (1 + e^2)\mathcal{L}_e]}, \end{gather}
\begin{gather}\mathcal{Y}^A =\frac{16 e^3}{3[ 2e + (3 e^2 - 1)\mathcal{L}_e]}, \end{gather}
\begin{gather}\mathcal{X}^C =\frac{4e^3(1 - e^2)}{3 [ 2e - ( 1 - e^2)\mathcal{L}_e ]}, \end{gather}
\begin{gather}\mathcal{Y}^C =\frac{4e^3(2 - e^2)}{ 3 [{-}2e + ( 1 + e^2)\mathcal{L}_e ]} \end{gather}
are functions of eccentricity 
$e$.
Appendix C. Disturbance viscosity effects
Variations in viscosity generally arise from variations in an underlying field that affects the viscosity, such as temperature, salt or nutrient concentration. Taking into account the effect of boundary conditions on the surface of the particle in relation to the underlying field will then lead to changes in the viscosity. For example, in an otherwise linear salt concentration field, the presence of a particle will disrupt the field due to salt impermeability. 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 power dissipation and swimming efficiency, factoring in the disturbance to viscosity due to the application of a no-flux condition at the particle's surface. The total viscosity field can be represented as the superposition of an ambient viscosity field (denoted as 
$\eta _0$) and a disturbance viscosity field (denoted by a prime):
\begin{equation} \eta = \eta_0 + \eta'. \end{equation} The boundary conditions on the disturbance viscosity field 
$\eta '$ are that the disturbance viscosity should diminish in the far-field region:
\begin{equation} \eta' \rightarrow 0 \quad \text{as} \ |\boldsymbol{r}| \rightarrow \infty. \end{equation}The transport of scalar fields, such as temperature, salt or nutrient concentrations, can be characterized using an advection–diffusion equation. When variations in these scalar fields are minimal, changes in viscosity can be modelled with a similar advection–diffusion equation. For slow-moving microswimmers in environments with highly diffusive scalar fields, such as temperature or salt concentration, advection is minimal, resulting in the viscosity distribution following the Laplace equation. Since the ambient viscosity field is linear, the disturbance viscosity field must also conform to the Laplace equation,
\begin{equation} \nabla^2 \eta = \nabla^2 \eta' = 0. \end{equation}The disturbance viscosity field is also governed by the boundary conditions on the particle's surface. We consider the surface to be impermeable to salt or insulated against temperature changes, and this leads to an associated ‘no-flux’ condition for viscosity at the particle's boundary,
\begin{equation} \boldsymbol{n} \boldsymbol{\cdot}\boldsymbol{\nabla} \eta = 0 \quad \text{on} \ S_p. \end{equation} Supposing that the ambient viscosity field is along 
$\boldsymbol {e}_1$, the mathematical expression of the disturbance viscosity field 
$\eta '$ in a particle-aligned coordinate 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} = \epsilon \eta_{\infty} p_1\,\frac{2e (1 - e^2)}{ [ 2e - (1 - e^2) \mathcal{L}_e]}, \end{gather}
\begin{gather}A_{1,1} = \epsilon \eta_{\infty} \sqrt{1 - p_1^2}\,\frac{2e (1 - e^2)}{ [ 2e - 4 e^3 - (1 - e^2) \mathcal{L}_e]}. \end{gather}
Here, 
$p_1 = \boldsymbol {p} \boldsymbol {\cdot } \boldsymbol {e}_1$, and 
$P_k^m$ and 
$Q_k^m$ denote the associated Legendre polynomials of the first and second kinds, with degree 
$k$ and order 
$m$. Their detailed mathematical formulations can be found in Abramowitz & Stegun (Reference Abramowitz and Stegun1964).
The power dissipation of a prolate spheroidal squirmer with a no-flux condition is, to first order,
\begin{equation} \mathcal{P} ={\rm \pi} (\mathcal{A} + \mathcal{B} \beta^2)\,\eta(\boldsymbol{x}_c)\,a B_1^2 + {\rm \pi}a^2 \mathcal{C}^{nf} B_1^2 \beta( \boldsymbol{p} \boldsymbol{\cdot} \boldsymbol{\nabla} \eta), \end{equation}where
\begin{equation} \mathcal{C}^{nf} = \frac{8({-}1 + e^2)[6e + ({-}3 + e^2) \mathcal{L}_e]}{e^2[2e + ({-}1 + e^2) \mathcal{L}_e]} \end{equation}
is a monotonically decreasing function of eccentricity. For a sphere, 
$e \rightarrow 0$, 
$\mathcal {C}^{nf} = 16/5$, and we have
\begin{equation} \mathcal{P}_{{sphere}}^{nf} = \tfrac{8}{3} {\rm \pi}(2 + \beta^2)\,\eta(\boldsymbol{x}_c)\, a B_1^2 + \tfrac{16}{5} {\rm \pi}a^2 B_1^2\beta( \boldsymbol{p} \boldsymbol{\cdot} \boldsymbol{\nabla}\eta). \end{equation}
The corresponding swimming efficiency up to 
$O(\epsilon )$ is
\begin{equation} \mathcal{E} = \frac{\mathcal{F}}{\mathcal{A} + \mathcal{B} \beta^2} + \epsilon\, \frac{\mathcal{G}^{nf} + \mathcal{H}^{nf} \beta^2}{(\mathcal{A} + \mathcal{B} \beta^2)^2}\, \beta(\kern0.7pt\boldsymbol{p} \boldsymbol{\cdot} \boldsymbol{d}), \end{equation}where
\begin{align} \mathcal{G}^{nf} & = 8(1-e^2)[8e^3(45+87e^2-32e^4)+4e^2({-}135 - 36e^2 + 19e^4 + 32e^6)\mathcal{L}_e \nonumber\\ & \quad + (270e\! -\! 426e^3\! +\! 402e^5\! -\! 182e^7)\mathcal{L}_e^2 + ({-}45\! +\! 186e^2 \! -\! 212e^4\! +\! 54e^6\! +\! 17e^8)\mathcal{L}_e^3 \nonumber\\ & \quad - 12e (1-e^2)^2 (1+e^2) \mathcal{L}_e^4 ]\times \{3e^5 [6e + ({-}3 + e^2) \mathcal{L}_e][{-}2e + (1 + e^2) \mathcal{L}_e] \}^{{-}1}, \end{align}
\begin{align} \mathcal{H}^{nf} & = 8 [{-}4e^2(45\! -\! 51e^2\! +\! 8e^4)\! +\! 12e(15\! -\! 22e^2\! +\! 7e^4)\mathcal{L}_e\! -\! 3(15\! -\! 27e^2\! +\! 13e^4\! -\! e^6)\mathcal{L}_e^2] \nonumber\\ & \quad\times [4e^2(63\! -\! 87e^2\! +\! 32e^4)\! -\! 4e(63\! -\! 96e^2\! +\! 37e^4)\mathcal{L}_e\! +\! (63 \! -\! 81e^2\! +\! e^4\! +\! 17e^6)\mathcal{L}_e^2\nonumber \\ & \quad -12e(1-e^2)^2 \mathcal{L}_e^3] \times \{9e^{7} [6e + ({-}3 + e^2) \mathcal{L}_e]^2 [{-}2e + (1 + e^2) \mathcal{L}_e] \}^{{-}1}. \end{align}
Here, 
$\mathcal {G}^{nf}$ (for 
$e < 0.954$) and 
$\mathcal {H}^{nf}$ are monotonically decreasing functions of slenderness, and they both vanish when 
$e \rightarrow 1$. For a sphere, we obtain
\begin{equation} \mathcal{E}_{{sphere}}^{nf} = \frac{1}{2 + \beta^2} + \epsilon\,\frac{14 + 13 \beta^2}{10(2+ \beta^2)^2}\,\beta(\kern0.7pt\boldsymbol{p} \boldsymbol{\cdot} \boldsymbol{d}). \end{equation}We see that taking into account the disturbance viscosity due to a no-flux boundary condition does not change the qualitative picture when compared to the effects of ambient viscosity alone.
Appendix D. Higher-order modes
A spherical squirmer with an arbitrary number of squirming modes 
$B_n$ expends power
\begin{align} \mathcal{P}_{{sphere}} &= \left( \frac{16}{3}\,B_1^2 + \sum_{n = 2}^{\infty} \frac{16}{n( n+ 1 )}\,B_n^2 \right) {\rm \pi}\,\eta ( \boldsymbol{x}_c )\,a\nonumber\\ &\quad + \left( \frac{32}{15}\,B_1 B_2 + \sum_{n = 2}^{\infty} \frac{32}{( n+ 1 ) (2n+1 )}\,B_n B_{n+1} \right) {\rm \pi}a^2 ( \boldsymbol{p} \boldsymbol{\cdot} \boldsymbol{\nabla} \eta ), \end{align}and swims with efficiency
\begin{align} \mathcal{E}_{{sphere}} &= \frac{B_1^2}{2 B_1^2 + \displaystyle\sum_{n = 2}^{\infty} \dfrac{6}{n (n + 1 )}\,B_n^2} \nonumber\\ &\quad + \epsilon ( \boldsymbol{p} \boldsymbol{\cdot} \boldsymbol{d} )\,\frac{ \dfrac{8}{5}\,B_1^3 B_2 + \displaystyle\sum_{n = 2}^{\infty} \left( \dfrac{36}{5n ( n+ 2 ) }\,B_1 B_2 B_n^2 - \dfrac{12}{( n+ 1 ) ( 2n+1 )}\,B_1^2 B_n B_{n+1} \right) }{\left( 2 B_1^2 + \displaystyle\sum_{n = 2}^{\infty} \dfrac{6}{n (n + 1 )}\,B_n^2 \right)^2}. \end{align}
