1. Introduction
Fluids with inhomogeneous viscosity fields are ubiquitous around us. For example, certain biological fluids like mucus and extracellular microbial polymers are mixtures of fluids with different viscosities (Berg 2004), and therefore, exhibit variable viscosity, either with (Esparza López et al. Reference Esparza López, Gonzalez-Gutierrez, Solorio-Ordaz, Lauga and Zenit2021) or without sharp viscosity gradients (Du et al. Reference Du, Keener, Guy and Fogelson2012). Similarly, gradients in temperature, salinity or concentration may induce spatial variation in viscosity, most commonly observed in marine ecosystems (Arrigo et al. Reference Arrigo, Robinson, Worthen, Dunbar, DiTullio, VanWoert and Lizotte1999). Finally, suspensions of particles in Newtonian fluids (both active and passive) may be treated at the continuum level as fluids with viscosity varying with local volume fraction (Hatwalne et al. Reference Hatwalne, Ramaswamy, Rao and Simha2004; Rafaï, Jibuti & Peyla Reference Rafaï, Jibuti and Peyla2010).
In this paper we examine an idealized problem of a single spheroid sedimenting in a spatially varying viscosity field. We discuss the dynamics that are observed, and how they differ from other situations studied in the literature. By now, it is well known that in Stokes flow, a spheroid in gravity does not change its orientation due to the particle symmetry and the reversibility of the Stokes equations. If the orientation starts out neither parallel or perpendicular to the gravity direction, the particle will move in a straight diagonal line, the direction of which is determined by the resistances parallel and perpendicular to the particle's orientation vector (figure 1a). These dynamics will change only when symmetry breaking is present in the system. One way in which symmetry breaking occurs is if fluid inertia is present (Cox Reference Cox1965; Khayat & Cox Reference Khayat and Cox1989; Auguste, Magnaudet & Fabre Reference Auguste, Magnaudet and Fabre2013), or if the suspending fluid has normal stresses due to the presence of polymers (Kim Reference Kim1986; Galdi Reference Galdi2000; Galdi et al. Reference Galdi, Vaidya, Pokorný, Joseph and Feng2011). For example, small fluid inertia generates a torque that orients the spheroid's longest axis perpendicular to the external force – the so-called ‘broad side on’ configuration (Dabade, Marath & Subramanian Reference Dabade, Marath and Subramanian2015). Conversely, fluid viscoelasticity orients the spheroid such that its longest axis is along the force direction, i.e. an ‘edge wise’ configuration (Kim Reference Kim1986; Dabade et al. Reference Dabade, Marath and Subramanian2015). These effects markedly change the particle trajectory as well as the sedimentation speed (figure 1), since the particle's drag coefficient is a function of orientation and is minimized when the longest axis is along the force direction.
Figure 1. Illustration of spheroid orientation and trajectory during sedimentation in (a) Stokes flow (zero Reynolds number), (b) fluid with finite inertia and (c) polymeric fluid with normal stresses (large elasticity number). This paper investigates the behaviour when viscosity stratification is present, i.e. case (d).
Another way in which symmetry breaking could occur is if there is a stratified fluid, i.e. variations in density, viscosity or other fluid properties that alter the force and torque on the particle (More & Ardekani Reference More and Ardekani2022). This area of research is relatively modern, and most of the efforts have examined the effect of density stratification on particle dynamics including both solid particles (Doostmohammadi, Dabiri & Ardekani Reference Doostmohammadi, Dabiri and Ardekani2014; Ardekani, Doostmohammadi & Desai Reference Ardekani, Doostmohammadi and Desai2017) and liquid droplets (Mandel et al. Reference Mandel, Zhou, Waldrop, Theillard, Kleckner and Khatri2020). When density increases along the gravity direction, it is found that the drag on a sphere is enhanced as confirmed by theory (Mehaddi, Candelier & Mehlig Reference Mehaddi, Candelier and Mehlig2018), experiments (Lofquist & Purtell Reference Lofquist and Purtell1984; Yick et al. Reference Yick, Torres, Peacock and Stocker2009) and simulations (Hanazaki, Konishi & Okamura Reference Hanazaki, Konishi and Okamura2009; More et al. Reference More, Ardekani, Brandt and Ardekani2021). The buoyancy force also leads to continuous deceleration and absence of a terminal velocity (Doostmohammadi et al. Reference Doostmohammadi, Dabiri and Ardekani2014). For anisotropic particles like spheroids, there has been some research to understand their settling behaviour in density stratified fluids. Using a reciprocal theorem based approach, Varanasi & Subramanian (Reference Varanasi and Subramanian2022) showed that the hydrostatic torque due to buoyancy originating from density stratification tends to rotate the particle in a broad side on configuration (similar to inertia), which had earlier been also shown by Dandekar, Shaik & Ardekani (Reference Dandekar, Shaik and Ardekani2020). In the cited papers, it was assumed that the fluid density is not altered by the presence of the particle, and gives rise to a so-called ‘hydrostatic torque’. However, the particle itself can alter the density field, and this additional effect can modify the particle torque (More et al. Reference More, Ardekani, Brandt and Ardekani2021; Varanasi & Subramanian Reference Varanasi and Subramanian2022). For example, density is often linked to a scalar field like temperature, which depends on a convection–diffusion equation. Depending on the Péclet number, the density around the particle may or may not be coupled with the fluid flow. In the low-Péclet-number limit, this additional torque is opposite to the hydrostatic torque (More et al. Reference More, Ardekani, Brandt and Ardekani2021; Varanasi & Subramanian Reference Varanasi and Subramanian2022).
Despite the advances in understanding microhydrodynamics of particles in density stratified fluids, there is a relative lack of literature examining viscosity-stratified fluids, even though there is recent evidence suggesting that these effects would be more important than those due to variations in density in a variety of applications (Jacquemin et al. Reference Jacquemin, Husson, Padua and Majer2006; Dandekar & Ardekani Reference Dandekar and Ardekani2020). For example, viscosity gradients are present in the swimming of micro-organisms, and it is of much interest to biologists to understand how organisms move in such complex environments (Hatwalne et al. Reference Hatwalne, Ramaswamy, Rao and Simha2004; Sokolov & Aranson Reference Sokolov and Aranson2009; Rafaï et al. Reference Rafaï, Jibuti and Peyla2010; Sengupta, Carrara & Stocker Reference Sengupta, Carrara and Stocker2017; Liebchen et al. Reference Liebchen, Monderkamp, ten Hagen and Löwen2018), as well as roboticists who design microrobots in such fluids (Nelson, Kaliakatsos & Abbott Reference Nelson, Kaliakatsos and Abbott2010; Palagi et al. Reference Palagi, Jager, Mazzolai and Beccai2013; Kim et al. Reference Kim, Lee, Lee, Nelson, Zhang and Choi2016; Li et al. Reference Li, de Ávila, Gao, Zhang and Wang2017a; Zhuang, Park & Sitti Reference Zhuang, Park and Sitti2017; Asghar et al. Reference Asghar, Javid, Waqas, Ghaffari and Khan2020). Some questions that arise are: how does a spatially varying fluid viscosity affect the common swimming speed, propulsion and efficiency (Gagnon & Arratia Reference Gagnon and Arratia2016)? Do microswimmers orient themselves in preferable positions in response to the viscosity gradients (Takabe et al. Reference Takabe, Tahara, Islam, Affroze, Kudo and Nakamura2017)? The common approach is to leverage a prototypical swimmer model (for squirmers, see Datt & Elfring Reference Datt and Elfring2019 and Shaik & Elfring Reference Shaik and Elfring2021; for a swimming sheet, see Dandekar & Ardekani Reference Dandekar and Ardekani2020 and Eastham & Shoele Reference Eastham and Shoele2020; for Purcell's swimmer, see Qin & Pak Reference Qin and Pak2023; for cilia, see Palagi et al. Reference Palagi, Jager, Mazzolai and Beccai2013 and Asghar et al. Reference Asghar, Javid, Waqas, Ghaffari and Khan2020) and then couple it to the Stokes flow field with a variable viscosity.
Currently, work has been performed on the motion of a single sphere in a viscosity varying fluid (Datt & Elfring Reference Datt and Elfring2019), but the effect of particle shape has yet to be considered. We note that the authors in the cited paper found that viscosity gradients give rise to force/rotation and torque/translation coupling for the sphere's motion, which would otherwise not exist if the viscosity gradient were absent. This type of coupling is likely to give rise to unique rotational dynamics for orientable particles, which we investigate in this paper. We note that similar to the case of density stratified flows, the motion of the particle may induce changes in viscosity in the fluid around the particle that could lead to novel hydrodynamic phenomena. For the case of spheres, as shown by Shaik & Elfring (Reference Shaik and Elfring2021), the changes in viscosity induced by particle motion do not lead to qualitative changes in the particle dynamics. Thus, in our paper, we do not account for changes in viscosity induced by the motion of the particle and only consider a spatially varying, time invariant viscosity field. However, later in § 5.4, we expound the probable strategy to incorporate the changes in viscosity due to particle motion, in the low-Péclet-number (
$Pe$) limit.
With this motivation in mind, in this paper we examine a problem of a single spheroid sedimenting in a Newtonian fluid with a spatially varying viscosity field. The viscosity field varies linearly in space, and its gradient points in an arbitrary direction with respect to the direction of sedimentation (external force). Section 2 outlines the particle geometry and equations of motion. Section 3 numerically solves for the particle's rigid body motion in the limit of weak viscosity gradient using the reciprocal theorem. Section 4 uses the principles of symmetry to obtain a general expression for the particle mobility tensor relating the particle's rigid body motion with the force and torque on the spheroid. The force/translation and torque/rotation relationships are unaltered due to the presence of a viscosity gradient, but the viscosity gradient gives rise to new force/rotation and torque/translation coupling terms that depend on three undetermined coefficients. We determine the values of these coefficients numerically, and thus, are able to solve the rigid body problem for an arbitrary set of forcing, viscosity gradient direction and particle geometry. Section 5 discusses some illustrative examples, wherein the orientations and trajectories of settling spheroids are analysed for different directions of the viscosity gradient. We find that depending on the viscosity gradient direction, particle shape (prolate versus oblate spheroid) and particle aspect ratio, the spheroid can take on different steady orientation angles. The section concludes on how to extend the analysis to more complicated situations, followed by § 6 which summarizes all results.
We note that although this work primarily focuses on passive particles in viscosity-stratified fluids, the results here will likely be important in a variety of contexts beyond this work. For example, scientists are interested in quantifying the swimming of particles in viscosity varying fluids, and the mobility relationships developed here can be used for such applications. Furthermore, understanding the rotation behaviour and velocity field from a single, orientable particle can help understand their far-field hydrodynamic interactions in a dilute suspension, which is important in understanding concentration instabilities that arise in fibrous suspensions (Koch & Shaqfeh Reference Koch and Shaqfeh1989, Reference Koch and Shaqfeh1991; Nicolai et al. Reference Nicolai, Herzhaft, Hinch, Oger and Guazzelli1998; Herzhaft & Guazzelli Reference Herzhaft and Guazzelli1999; Butler & Shaqfeh Reference Butler and Shaqfeh2002; Kuusela, Lahtinen & Ala-Nissila Reference Kuusela, Lahtinen and Ala-Nissila2003; Shin, Koch & Subramanian Reference Shin, Koch and Subramanian2006, Reference Shin, Koch and Subramanian2009; Vishnampet & Saintillan Reference Vishnampet and Saintillan2012). The force/rotation and torque/translation couplings shown in our paper have also been very recently shown in sedimentation of spheroids with non-uniform density, where the centre of mass of the particle does not match its gravitational centre (Nissanka, Ma & Burton Reference Nissanka, Ma and Burton2023) and, therefore, we expect that the theory developed in this paper may also aid in the understanding of the sedimentation of such mass polar spheroids as well. We will not comment on this point further, noting that the work acts as a stepping stone for these more complicated problems when viscosity gradients are present.
2. Problem statement
2.1. Problem geometry
The schematic of our system is shown in figures 2 and 3. We consider a torque-free spheroid under an external force 
$\boldsymbol {F}$ in a Newtonian fluid with a constant viscosity gradient 
$\boldsymbol {\nabla } \eta$. The force is in the positive 3-direction. The viscosity gradient 
$\boldsymbol {\nabla } \eta$ can be co-linear with the force (figure 2, where 
$\boldsymbol {\nabla } \eta$ is in the 
$\pm$3-direction) or perpendicular to the force (figure 3, where 
$\boldsymbol {\nabla } \eta$ is in the 
$+$1-direction). The spheroid has three semi-major axes of lengths (
$a, b, c$), with 
$a\neq b = c$. The initial centre of mass of the spheroid is 
$(x_{01},x_{02},x_{03})=(0,0,0)$.
Figure 2. Schematic of (a) a prolate and (b) oblate spheroid falling under an external force acting in the 3-direction. The viscosity gradient is along the 3-direction (parallel or anti-parallel). The particle's orientation vector 
$\boldsymbol {p}$ makes a polar angle 
$\alpha \in [0,{\rm \pi} ]$ with respect to the sedimentation direction.
Figure 3. Schematic of (a) a prolate and (b) oblate spheroid falling under an external force 
$\boldsymbol {F}$ acting in the 
$3$-direction, while the viscosity varies spatially in the 
$1$-direction. The particle's orientation 
$\boldsymbol {p}$ makes a polar angle 
$\alpha \in [0,{\rm \pi} ]$ with respect to the 3-direction, and makes an azimuthal angle 
$\phi \in [0,2{\rm \pi} )$ in the 1-2 plane.
We define the spheroid's orientation vector 
$\boldsymbol {p}$ as the direction along its unequal axis (i.e. the 
$a$ axis). Two different cases arise. A prolate spheroid has 
$\boldsymbol {p}$ along its longest axis, while an oblate spheroid has 
$\boldsymbol {p}$ along its shortest axis. Another way to parameterize the particle shape is through an aspect ratio parameter 
$A_R$ and equivalent radius 
$R$. Here, 
$A_R$ is the ratio 
$a/b$, while 
$R$ is the radius of an equivalent sphere with the same volume.
\begin{equation} A_R =\frac{a}{b}, \quad R =(abc)^{1/3}. \end{equation}The two systems of parameterization are connected by the following relationship:
\begin{equation} a=R A_{R}^{2 / 3}, \quad b=c= R A_R^{-1/3}. \end{equation} Evidently, a prolate particle has its aspect ratio parameter 
$A_R > 1$, while an oblate particle has its aspect ratio parameter 
$A_R < 1$. Figures 2 and 3 describe the polar and azimuthal angles 
$\alpha \in [0,{\rm \pi} ]$ and 
$\phi \in [0,2{\rm \pi} )$ for the particle orientation. The next section discusses the equations of motion and the rheology of the fluid.
2.2. Equations of motion and fluid rheology
The fluid surrounding the particle is incompressible and Newtonian. The fluid also has negligible inertia – in other words, the Reynolds number based on the particle's largest length scale 
$Re=(\rho _f U L_{max})/\eta _0 \approx 0$. Here, 
$\rho _f$ is the density of the fluid surrounding the particle, 
$U$ is the translation speed of the particle, 
$L_{max} =\max (a,b)$ is the largest axis of the particle and 
$\eta _0$ is the fluid's viscosity at the origin if the particle were absent.
When these conditions hold, the momentum and mass balance equations in the fluid are given as
\begin{equation} \frac{\partial \sigma_{ij}}{\partial x_j}=0, \quad \frac{\partial v_i}{\partial x_i}=0, \end{equation}
where 
$\sigma _{ij}$ is the stress tensor and 
$v_i$ is the velocity field. Einstein summation convention is assumed, i.e. repeated indices are summed. The stress tensor takes the form
\begin{equation} \sigma_{ij} =-p\delta_{ij}+\eta(\boldsymbol{x})\dot{\gamma}_{ij}, \end{equation}
where 
$p$ is the pressure, 
$\dot {\gamma }_{ij} = {\partial v_j}/{\partial x_i}+{\partial v_i}/{\partial x_j}$ is twice the strain rate tensor and 
$\eta$ is the viscosity of the medium. In this problem, the viscosity is independent of the strain rate but exhibits a spatial dependence. The viscosity field is
\begin{equation} \eta(\boldsymbol{x}) =\eta_{0}\left( 1 +\frac{\beta}{R} \hat{\boldsymbol{d}} \boldsymbol{\cdot} \boldsymbol{x} \right)\!. \end{equation}
Here 
$\eta _0$ is the viscosity at the origin and 
$\boldsymbol {\nabla } \eta = ({\eta _0}/{R}) \beta \hat {\boldsymbol {d}}$ is a constant viscosity gradient with dimensionless magnitude 
$\beta$ and unit direction 
$\hat {\boldsymbol {d}}$.
The goal of the problem is to solve (2.3), (2.4) and (2.5) for the stress and velocity around the particle. The equations have to be solved with the boundary conditions
$$\begin{gather} v_i \rightarrow 0, \quad |x_i| \rightarrow \infty, \end{gather}$$
$$\begin{gather}v_i=U_i+\epsilon_{ijk}\omega_j (x_k-x_k^{cm}), \quad x_i \in S_p, \end{gather}$$
where 
$(U_i,\omega _i)$ are the rigid body velocities of the particle, 
$S_p$ is the particle surface, 
$x_k^{cm}$ is the centre of mass and 
$\epsilon _{ijk}$ is the Levi-Civita symbol. An additional constraint is that the particle's external force and torque are specified. These are
\begin{equation} F_i =-\int_{S_p}\sigma_{ij}n_j\,{\rm d}S, \quad T_i =-\int_{S_p}\epsilon_{ikl}(x_k-x_k^{cm})\sigma_{lj}n_j\,{\rm d}S, \end{equation}
where 
$n_i$ is the outward-pointing vector on the particle surface. For this problem, 
$T_i =0$.
In this problem, we specify the viscosity field to have a constant gradient, while for other problems, the viscosity field is often found by solving a scalar quantity like temperature or concentration that is a solution to a convection–diffusion equation around the particle. For such problems in the limit of a small Péclet number (one-way coupling), the results will be very similar to the problem formulated here, albeit with minor quantitative differences. A more detailed discussion will be provided at the end of the paper (§ 5.4).
2.3. Non-dimensionalization, dimensionless numbers and perturbation expansion
Unless otherwise noted, all quantities from here on out will be written in non-dimensional form. Lengths will be scaled by the average particle size 
$R$, forces by its magnitude 
$F$ and viscosities by its value at the origin 
$\eta _0$. Velocities will be scaled by the Newtonian sedimentation velocity 
$U ={F}/{6{\rm \pi} \eta _0 R}$, times by 
$t_c=R/U$, strain rates and rotational velocities by 
$\dot {\gamma }_c =1/t_c$, stresses and pressures by 
$\eta _0\dot {\gamma }_c$, and torques (if present) by 
$FR$.
The dynamics of the spheroid will depend on the following dimensionless quantities – the particle aspect ratio parameter 
$A_R$, the particle orientation 
$\boldsymbol {p}$ (characterized by angles 
$\alpha$ and 
$\phi$) and the non-dimensional viscosity gradient 
$\boldsymbol {\nabla } \eta$ (characterized by magnitude 
$\beta$ and direction 
$\hat {\boldsymbol {d}}$):
\begin{equation} A_R = \frac{b}{a}; \quad \boldsymbol{p} = [\sin\alpha \cos \phi, \sin \alpha \sin \phi, \cos \alpha]; \quad \boldsymbol{\nabla} \eta = \beta \hat{\boldsymbol{d}}. \end{equation}In dimensionless form, the viscosity of the fluids in figures 2 and 3 are
\begin{equation} \eta =1\pm \beta x_3, \quad \eta = 1+\beta x_1, \end{equation}
where the first equation, namely 
$\eta =1 \pm \beta x_3$, corresponds to the situation where the viscosity gradient is parallel (
$+3$-direction) or anti-parallel (
$-3$-direction) to the external force, and the second equation, namely 
$\eta =1 + \beta x_1$, corresponds to the case where the viscosity gradient is perpendicular to the external force. For a general viscosity gradient 
$\boldsymbol {\nabla } \eta$, the particle motion will be a superposition of the solutions for the two cases listed above. We will examine particle dynamics in the limit of a small viscosity gradient:
\begin{equation} Re \ll \beta \ll 1. \end{equation}
The above condition indicates that one can neglect fluid inertia and perform a regular perturbation expansion in 
$\beta$. We will solve for the rigid body motion up to 
$O(\beta )$, both numerically and semi-analytically using symmetry arguments listed in the next sections.
3. Numerical solution to particle dynamics
3.1. Reciprocal theorem
We determine the rigid body motion of the spheroid by performing a perturbation expansion in the non-dimensional viscosity gradient 
$\beta \ll 1$. We perturb the dependent variables as
\begin{align} {\{v_i, p, \sigma_{ij},\dot{\gamma}_{ij},F_i, U_i, \omega_i \}}&= \{v_{i}^{(0)}, p^{(0)}, \sigma_{ij}^{(0)},\dot{\gamma}_{ij}^{(0)},F_i^{(0)},U_i^{(0)}, \omega_i^{(0)}\}\nonumber\\ &\quad +\beta\{u_{i}^{(1)}, p^{(1)}, \tau_{ij}^{(1)},\dot{\gamma}_{ij}^{(1)},F_i^{(1)},U_i^{(1)}, \omega_i^{(1)}\}+\cdots, \end{align}
and solve for the momentum and mass balances (2.3)–(2.7) at each order in 
$\beta$. At leading order, the spheroid sediments in a zero-Reynolds-number fluid with a constant, non-dimensional viscosity 
$\eta =1$ and a non-dimensional external force 
$F = 1$:
\begin{equation} \frac{\partial ^2 v_i^{(0)}}{\partial x_k \partial x_k}-\frac{\partial p^{(0)}}{\partial x_i} =0, \quad \frac{\partial v_i^{(0)}}{\partial x_i} =0,\quad F_i =\delta_{i3} ,\quad T_i =0. \end{equation}The solution to the above problem is given in many classical texts (for example, see Kim & Karilla Reference Kim and Karilla2005). The velocity field is presented in Appendix A, while the rigid body motion satisfies the classical resistance relationship
\begin{equation} \left(\begin{array}{@{}cc@{}}
\boldsymbol{R}^{FU} & \boldsymbol{R}^{F\omega} \\
\boldsymbol{R}^{TU} & \boldsymbol{R}^{T\omega} \\
\end{array}\right) {\cdot} \left(\begin{array}{@{}c@{}}
\boldsymbol{U}^{(0)} \\ \omega^{(0)} \\
\end{array}\right)=\left(\begin{array}{@{}c@{}} \boldsymbol{F} \\
\boldsymbol{T} \\ \end{array}\right)\!.
\end{equation}
Here 
$(\boldsymbol {R}^{FU},\boldsymbol {R}^{F\omega },\boldsymbol {R}^{TU},\boldsymbol {R}^{T\omega })$ are the resistance tensors for a spheroid, which are given in Appendix B. The external force and torque are given in (3.2).
At the next order of approximation 
$O(\beta )$, the momentum and mass balance equations become Stokes flow with an extra fluid body force 
$b_i$:
\begin{equation} \frac{\partial ^2v_i^{(1)}}{\partial x_k\partial x_k} -\frac{\partial p^{(1)}}{\partial x_i} +b_i =0,\quad \frac{\partial v_i^{(1)}}{\partial x_i} =0. \end{equation}Here the body force is due to the spatially varying viscosity field:
\begin{equation} b_i =\frac{\partial\tau_{ij}^{ex}}{\partial x_j}, \quad \tau_{ij}^{ex} = (\hat{\boldsymbol{d}} \boldsymbol{\cdot} \boldsymbol{x})\dot{\gamma}_{ij}^{(0)}. \end{equation}
In (3.5), 
$\tau _{ij}^{ex}$ is the extra stress tensor and 
$\dot {\gamma }^{(0)}_{ij}$ is twice the rate of strain tensor from the leading-order velocity field.
We employ the reciprocal theorem to solve for the translational and rotational velocity for the 
$O(\beta )$ problem. This theorem has a storied history in the Stokes flow community, as it is often used to solve for the rigid body motion of particles in Stokes flow with a fluid body force. The derivation is stated in Appendix C and we present the main results below. In brief, the translational and rotational velocities follow a resistance relationship similar to (3.3), except the forces and torques are replaced by an effective viscosity-stratified force and torque:
\begin{equation} \left(\begin{array}{@{}cc@{}}
\boldsymbol{R}^{FU} & \boldsymbol{R}^{F\omega} \\
\boldsymbol{R}^{TU} & \boldsymbol{R}^{T\omega} \\
\end{array}\right) {\cdot} \left(\begin{array}{@{}c@{}}
\boldsymbol{U}^{(1)} \\ \omega^{(1)} \\
\end{array}\right)=\left(\begin{array}{@{}c@{}}
\boldsymbol{F}^{vs} \\ \boldsymbol{T}^{vs} \\
\end{array}\right)\!.
\end{equation}The viscosity-stratified force and torque are given as
\begin{equation} F_i^{vs} =-\int_{V_{out}}\frac{\partial v_{ki}^{trans}}{\partial x_j}\tau_{kj}^{ex} \,{\rm d}V, \quad T_i^{vs}=-\int_{V_{out}}\frac{\partial v_{ki}^{rot}}{\partial x_j}\tau_{kj}^{ex} \,{\rm d}V, \end{equation}
where the integrals are evaluated over the volume 
$V_{out}$ outside the particle. The quantities 
$v_{ki}^{trans}$ and 
$v_{ki}^{rot}$ are the Stokes flow velocity fields around a spheroid in the 
$k$ direction due to unit translation or unit rotation in the 
$i$ direction. These quantities are derived from the same velocity fields listed in Appendix A.
3.2. Numerical implementation
The volume integrals in (3.7) are difficult to evaluate analytically. A custom-made MATLAB code was written to calculate the spheroid's rigid body motion. This code is similar to the approach used in our prior papers to investigate the motion of ellipsoids in weakly viscoelastic fluids (Wang, Tai & Narsimhan Reference Wang, Tai and Narsimhan2020), except that here the extra stress is modified to account for the viscosity gradient. First, we transform from the laboratory frame to the particle frame of reference where the origin is at the particle's centre of mass and the Cartesian coordinate axes align with the particle's principle axes. We then evaluate the volume integrals in (3.7) for the viscosity-stratified force and torque, using an elliptical coordinate system and performing Gaussian quadrature via Legendre polynomials. A mesh convergence study showed that 
$60$ elements in the radial direction, 
$80$ elements in the polar direction and 
$100$ elements in the azimuthal direction yielded a sufficiently accurate mesh for our analysis. We then solve the matrix equations (3.3) and (3.6) for the rigid body motions at 
$O(1)$ and 
$O(\beta )$, and transform back to the laboratory frame. The particle's centre of mass and orientation are evolved by solving the rigid body dynamics
\begin{equation} \frac{{\rm d}\,\boldsymbol{x}^{cm}}{{\rm d}t} = \boldsymbol{U},\quad \frac{{\rm d} \boldsymbol{p}}{{\rm d}t} = \boldsymbol{\omega} \times \boldsymbol{p}. \end{equation}
We use a forward Euler scheme with 
$\Delta t = 0.01$. More details are found in our prior publications (Wang et al. Reference Wang, Tai and Narsimhan2020; Anand & Narsimhan Reference Anand and Narsimhan2023).
3.2.1. Verification of code
For the case of a sphere sedimenting in a linear, imposed viscosity gradient, we refer to the work by Datt & Elfring (Reference Datt and Elfring2019). Specifically, Eqs.(7,8) in Datt & Elfring (Reference Datt and Elfring2019) are the resistance relationships for the external force 
$\boldsymbol {F}$ and torque 
$\boldsymbol {T}$ on a sphere of radius 
$a$ in a fluid with a constant viscosity gradient 
$\boldsymbol {\nabla } \eta$, with translational velocity 
$\boldsymbol {U}$ and rotational velocity 
$\boldsymbol {\omega }$. For convenience, these equations are reproduced in dimensional form here:
$$\begin{gather} \boldsymbol{F}= 6 {\rm \pi}a \eta_{0} \boldsymbol{U} - 2 {\rm \pi}a^{3} \boldsymbol{\nabla} \eta \times \boldsymbol{\omega}, \end{gather}$$
$$\begin{gather}\boldsymbol{T}= 2 {\rm \pi}a^{3} \boldsymbol{\nabla} \eta \times \boldsymbol{U} + 8 {\rm \pi}\eta_{0} a^{3} \boldsymbol{\omega}. \end{gather}$$We employ the above equations ((3.9) and (3.10)) to validate our theory for two different cases namely, spatial variation of viscosity in y direction (i) and spatial variation of viscosity in x direction (ii), as illustrated below.
(i) Spatial variation in the
$y$ direction: For a torque-free $(\boldsymbol {T} = 0)$ sphere sedimenting in the $x$ direction where the dimensional viscosity gradient is along the $y$ direction $\boldsymbol {\nabla }\eta = \beta ({\eta _0}/{a}) \hat {\boldsymbol {y}}$, the above equations give us
(3.11a,b)$$\begin{gather} U_x = \frac{F_x}{{\rm \pi} a \eta_0 (6-0.5\beta^2)}, \quad U_y = U_z = 0, \end{gather}$$(3.12a,b)$$\begin{gather}\omega_x = \omega_y = 0, \quad \omega_z = 0.25\beta \frac{F_x}{{\rm \pi} \eta_0 a^2(6-0.5\beta^2)}. \end{gather}$$The analytical result for
$\omega _z$ in (3.12) is compared against the results of the numerical simulation and the comparison is shown in figure 4(a) showing an accurate match.(ii) Spatial variation in the
$x$ direction: Similarly, for a torque-free sphere sedimenting in the $x$ direction where the dimensional viscosity gradient is along the $x$ direction $\boldsymbol {\nabla } \eta = \beta ({\eta _0}/{a}) \hat {\boldsymbol {x}}$, the above equations give
(3.13a,b)$$\begin{gather} U_x = \frac{F_x}{6{\rm \pi} a \eta_0}, \quad U_y = U_z = 0, \end{gather}$$(3.14)We compare the results of the simulation against (3.13) in figure 4(b) where a good match is seen. The relative error in the simulation with respect to the simulation results is <1 %.$$\begin{gather}\omega_x = \omega_y = \omega_z = 0. \end{gather}$$
Figure 4. Code validation for a sphere sedimenting in a fluid with a prescribed viscosity gradient in the (a) 
$y$ direction and (b) 
$x$ direction. For all the cases, the external force is a unit vector acting in the 
$x$ direction, while the external torque is 
$\boldsymbol {T}=0$. The radius and fluid viscosity are 
$a = 1$ and 
$\eta _0 = 1$, respectively. The results of the theory are from Datt & Elfring (Reference Datt and Elfring2019), expanded in § 3.2.1. Results are shown for (a) 
$\boldsymbol {F} \propto \hat {\boldsymbol {x}}$, 
$\boldsymbol {\nabla } \eta \propto \beta \hat {\boldsymbol {y}}$ and (b) 
$\boldsymbol {F} \propto \hat {\boldsymbol {x}}$, 
$\boldsymbol {\nabla } \eta \propto \beta \hat {\boldsymbol {x}}$.
4. Semi-analytical theory
4.1. Introduction and motivation
The simulations described in the previous section solve the rigid body motion of the particle, but are computationally intensive. At each time step, one has to evaluate six volume integrals in (3.7) to obtain the viscosity-stratified force and torque. Furthermore, a new time sweep has to be performed if one examines a different viscosity gradient direction and magnitude.
An alternative approach to obtain the same dynamics is to develop a semi-analytical theory based on the symmetry of the problem. Such a theory will give the general form of the particle's motion in terms of three undetermined constants, which in turn can be found by performing simulations at three specific configurations. The result of this analysis is that one can cheaply obtain the particle's motion for an arbitrary set of particle orientations, forcing and viscosity gradients.
What we are doing is essentially finding the general form of the mobility tensor when a viscosity gradient is present. Thus, the analysis below will not only give general information about the force/rotation coupling of these orientable particles, but can also give results for the case when a torque is applied; for example, the torque/translation coupling. A description is given below.
4.2. General form of mobility tensor
The governing momentum and continuity equations (2.3)–(2.7) are linear in the external force and torque 
$(\boldsymbol {F}, \boldsymbol {T})$. Thus, the translational and rotational velocities 
$(\boldsymbol {U},\boldsymbol {\omega })$ are also linear in these quantities and obey the following relationship:
\begin{equation} \left(\begin{array}{@{}c@{}}
\boldsymbol{U} \\ \omega
\end{array}\right)=\left(\begin{array}{@{}cc@{}} \boldsymbol{A} &
\boldsymbol{B} \\ \boldsymbol{B}^{T} & \boldsymbol{D}
\end{array}\right) {\cdot} \left(\begin{array}{@{}c@{}}
\boldsymbol{F} \\ \boldsymbol{T} \end{array}\right)\!.
\end{equation}
Here, 
$(\boldsymbol {A},\boldsymbol {B},\boldsymbol {D})$ are mobility tensors that are non-dimensionalized by 
$({U}/{F},{U}/{FR}, {U}/{FR^2}) = ({1}/{6{\rm \pi} \eta _0 R}, {1}/{6{\rm \pi} \eta _0 R^2}, {1}/{6{\rm \pi} \eta _0 R^3})$, respectively. In a constant viscosity fluid, these tensors are only a function of the particle shape and orientation, characterised by the aspect ratio parameter 
$A_R$ and the orientation vector 
$\boldsymbol {p}$. If a viscosity gradient is present, the tensors will also be a function of the non-dimensional viscosity gradient 
$\boldsymbol {\nabla }\eta = \beta \hat {\boldsymbol {d}}$. Note that the off-diagonal terms of the matrix in (4.1) are transposes of each other as can be proved by the reciprocal theorem (not shown here). Additionally, the mobility tensors 
$(\boldsymbol {A},\boldsymbol {D})$ are symmetric for an arbitrary shaped particle in a viscosity-stratified fluid, as has been shown with due rigour with aid of the reciprocal theorem in earlier literature (Oppenheimer, Navardi & Stone Reference Oppenheimer, Navardi and Stone2016).
In the limit of 
$\beta \ll 1$, we expand the mobility tensors in a Taylor series as follows:
\begin{equation} \{{A}_{ij},{B}_{ij},{D}_{ij}\} =\{{A}^{(0)}_{ij},{B}^{(0)}_{ij},{D}^{(0)}_{ij}\}+\beta\{{A}^{(1)}_{ij},{B}^{(1)}_{ij},{D}^{(1)}_{ij}\}. \end{equation} At leading order (
$O(1)$), the tensors are the same as those for the particle in a constant viscosity fluid. These quantities are well characterized and formulas are given in Appendix B for a general ellipsoid. Specifically, for the case where the particle has an orientation vector 
$\boldsymbol {p}$, they take the form
$$\begin{gather} {A}^{(0)}_{ij} = c_1(\delta_{ij}-{p_i}{p_j})+c_2{p_i}{p_j}, \end{gather}$$
$$\begin{gather}{D}^{(0)}_{ij} = c_3(\delta_{ij}-{p_i}{p_j})+c_4{p_i}{p_j} , \end{gather}$$
$$\begin{gather}{B}^{(0)}_{ij} =0, \end{gather}$$
where 
$c_1,c_2,c_3$ and 
$c_4$ are functions of the aspect ratio parameter and are given in Appendix B.
At 
$O(\beta )$, the motion will be linear in 
$\boldsymbol {\nabla } \eta$. Thus, in non-dimensional form, the mobility tensors take the following structure:
$$\begin{gather} A_{ij}^{(1)} =\alpha_{ijk}\hat{d}_k, \end{gather}$$
$$\begin{gather}D_{ij}^{(1)} =\beta_{ijk}\hat{d}_k , \end{gather}$$
$$\begin{gather}B_{ji}^{(1)} = M_{ikj}\hat{d}_k. \end{gather}$$
Here 
$\hat {d}_k$ is the direction of the viscosity gradient. Therefore, the problem of finding the mobility matrices 
$(A_{ij}^{(1)},B_{ij}^{(1)}, {D}_{ij}^{(1)})$ reduces to the problem of finding 
$\alpha _{ijk}$, 
$\beta _{ijk}$ and 
$M_{ijk}$. For a spheroid, these third-order tensors depend on the orientation product 
$p_ip_j$, since fore-aft symmetry dictates that changing 
$p_i$ to 
$-p_i$ will not alter the results. Noting that 
$(\alpha _{ijk}, \beta _{ijk})$ are third-order true tensors, and such tensors cannot be formed from 
$p_i p_j$, we obtain the result
\begin{equation} \alpha_{ijk} = \beta_{ijk} =0. \end{equation} The above relationship means that at 
$O(\beta )$ the force/velocity coupling and torque/angular velocity coupling are unchanged. However, as we will see next, the force/rotation coupling and torque/velocity coupling will change. Here 
$B_{ji}$ is a pseudo tensor since it connects a pseudo vector (angular velocity) with a true vector (force). Therefore, 
$M_{ikj}$ is a third-order pseudo tensor, which depends on the orientation product 
$p_i p_j$. The general form of 
$M_{ijk}$ is given below as
\begin{equation} M_{ijk}=\lambda_1\epsilon_{ijk}+\lambda_2p_i\epsilon_{jkq}p_q+\lambda_3p_j\epsilon_{ikq}p_q +\lambda_4p_k\epsilon_{ijq}p_q, \end{equation}
where 
$\lambda _1, \lambda _2, \lambda _3, \lambda _4$ are dimensionless coefficients that depend only on the aspect ratio parameter 
$A_R$. One can show that without loss of generality 
$\lambda _2=0$ (see Appendix D) and, therefore, the problem reduces to finding the coefficients 
$\lambda _1,\lambda _3,\lambda _4$. In other words, (4.6) reduces to
\begin{equation} M_{ijk} =\lambda_1\epsilon_{ijk}+\lambda_3p_j\epsilon_{ikq}p_q+\lambda_4p_k\epsilon_{ijq}p_q. \end{equation}We note that at the time of this paper being prepared for publication, another paper solving a similar problem arrived at a similar formulation for the mobility tensor (see Eq. 27 in Gong, Shaik & Elfring Reference Gong, Shaik and Elfring2023).
In summary, the mobility relationships up to 
$O(\beta )$ reduce to
$$\begin{gather} U_i = A_{ij}^{(0)}F_j+\beta M_{jki}\hat{d}_kT_j, \end{gather}$$
$$\begin{gather}\omega_i =\beta M_{ikj}\hat{d}_kF_j+D_{ij}^{(0)}T_j, \end{gather}$$
where 
$A_{ij}^{(0)}$, 
$D_{ij}^{(0)}$ are the known mobility tensors for a spheroid without a viscosity gradient, given by (4.3), while 
$M_{ijk}$ is the cross-coupling term given by (4.7). The unknown coefficients for the tensor 
$M_{ijk}$ are 
$(\lambda _1,\lambda _3,\lambda _4)$, which are functions of the aspect ratio parameter 
$A_R$ for the spheroid. We also note, from (4.7) in conjunction with (4.4c), that 
$\lambda _1,\lambda _3,\lambda _4$ are related to the eigenvalues of 
$\boldsymbol {B}$ and, in some cases of the orientation of the viscosity gradient vector, are exactly equal to the eigenvalues of 
$\boldsymbol {B}$. For a more detailed discussion on this interpretation, please refer to Witten & Diamant (Reference Witten and Diamant2020). The next section discusses how we determine these coefficients.
4.3. Determining coefficients $\lambda _1,\lambda _3,\lambda _4$ for the mobility matrix $M_{ijk}$ (force/rotation and torque/translation coupling)
Figure 5 outlines the simulations we perform to obtain the coefficients 
$(\lambda _1,\lambda _3,\lambda _4)$ for 
$M_{ijk}$ in (4.7). We examine a torque-free particle 
$(T_i = 0)$ and quantify its angular velocity 
$\omega _i$ for the three specific geometries listed below. We note that the angular velocity can cause two effects – it can change the spheroid's orientation or it can keep the orientation the same but spin it along its axis. The rate of change of the orientation is given by
\begin{equation} \frac{{\rm d} p_i}{{\rm d} t} =\epsilon_{ijk}\omega_j p_k, \end{equation}while the rate of spinning is
\begin{equation} \varOmega =\omega_ip_i. \end{equation}These quantities are computed for the cases below.
(i) Case A:
$\boldsymbol {\nabla } \eta \times \boldsymbol {F} =0$. Here, we examine the situation in figure 5(a) where the external force and viscosity gradient are in the same direction, i.e. $\boldsymbol {F} = \hat {\boldsymbol {d}} = \hat {\boldsymbol {z}}$. The particle has its orientation in the $x$–$z$ plane with an angle $\alpha$ with respect to the force direction, i.e. $\boldsymbol {p} =[\sin \alpha,0,\cos \alpha ]$. Using (4.7), (4.8b) and (4.9), one finds the angular velocity to be
(4.11)\begin{equation} {\omega_2=\frac{{\rm d} \alpha}{{\rm d} t} = \frac{1}{2} \beta(\lambda_3+\lambda_4) \sin({2\alpha})}. \end{equation}Thus, performing one simulation at a specific polar angle and viscosity gradient magnitude (say
$\alpha ={\rm \pi} /4, \beta =0.1)$ allows us to obtain $(\lambda _3+\lambda _4)$. To test this theory, we carried out simulations of prolate spheroids undergoing sedimentation at aspect ratios $A_R =3$ and $A_R = 5$. The simulations were performed for different angles $\alpha$ and non-dimensional viscosity gradient $\beta =0.1$, and they found that for a given aspect ratio, the quantity $2({({{\rm d} \alpha }/{{\rm d} t})}/{\beta \sin ({2\alpha })})$ was constant for all values of $\alpha$ (equal to $0.0649$ for $A_R=3$ and $0.0656$ for $A_R =5$). This behaviour is consistent with the expression listed above.(ii) Case B:
$\boldsymbol {\nabla } \eta \boldsymbol{\cdot} \boldsymbol {F} =0, \boldsymbol {p} \times \boldsymbol {F} = 0$. We examine the situation in figure 5(b) where the external force and viscosity gradient are perpendicular to each other, i.e. $\boldsymbol {F}=\hat {\boldsymbol {z}}$, $\hat {\boldsymbol {d}}= \hat {\boldsymbol {x}}$. The particle has its orientation along the force direction, i.e. $\boldsymbol {p}=[0,0,1]$, which corresponds to the polar and azimuthal angles of $\alpha = \phi = 0$ in figure 3. Using (4.7), (4.8b) and (4.9), one finds the angular velocity to be
(4.12)Performing one simulation at a specific value of\begin{equation} \omega_2 = \left.\frac{{\rm d} \alpha}{{\rm d} t}\right|_{\alpha = \phi = 0^{{\circ}}} =-\beta (\lambda_1 + \lambda_4). \end{equation}$\beta$ (e.g. $\beta =0.1$) allows us to obtain $(\lambda _1+\lambda _4)$.(iii) Case C:
$\boldsymbol {\nabla } \eta \perp \boldsymbol {F} \perp \boldsymbol {p}$. We examine the case in figure 5(c) where the orientation, viscosity gradient and force are all perpendicular to each other, i.e. $\boldsymbol {F}=\hat {\boldsymbol {z}}$, $\boldsymbol {\widehat {d }} = \hat {\boldsymbol {x}}$, $\boldsymbol {p}= \hat {\boldsymbol {y}}$. Here, the particle will spin but not change orientation. Using (4.7), (4.8b) and (4.10), we find the spinning rate to be
(4.13)Performing one simulation at a specific value of\begin{equation} \varOmega = \omega_2 =-\beta \lambda_1. \end{equation}$\beta$ (e.g. $\beta =0.1$) allows us to obtain $\lambda _1$.
The three simulations listed above yield a linear system of equations for the coefficients 
$(\lambda _1,\lambda _3,\lambda _4)$ that can be solved. Figure 6 shows the values of the coefficients for different values of the aspect ratio parameter 
$A_R$, for both prolate and oblate spheroids. For comparison, we have also plotted the 
$\lambda$s from Kamal & Lauga (Reference Kamal and Lauga2023), who invoked the slender body theory to solve the problem of sedimenting slender bodies in inertialess flows of Newtonian fluids with spatially varying viscosity. As can be seen from our plots, a good match is observed between our results for prolate spheroids at high 
$A_R$ and that given by Kamal & Lauga (Reference Kamal and Lauga2023).
Figure 6. Computed values of 
$(\lambda _1, \lambda _3, \lambda _4)$ for (a) prolate and (b) oblate spheroids for different values of aspect ratio parameters 
$A_R$. The empty symbols denote the results derived in this paper, while the filled symbols, for prolate spheroids, denote the results from Kamal & Lauga (Reference Kamal and Lauga2023), who used the slender body theory to derive their results.
Once these coefficients are tabulated, one has a general form for the rigid body motion (4.8) for spheroids that can be solved for an arbitrary viscosity gradient, orientation, aspect ratio and external force/torque.
5. Results and illustrative examples
In § 4 we developed a theory to describe the rigid body motion of a spheroid in a spatially varying viscosity field. In the limit of a weak or linearly varying viscosity field, the general form of the translational and rotational velocities is given in (4.8), where 
$A_{ij}^{(0)}$ and 
$D_{ij}^{(0)}$ are the standard mobility tensors for force/translation and torque/rotation in Stokes flow, and 
$M_{ijk}$ is a newly introduced coupling tensor between force/rotation and torque/translation that arises due to viscosity gradients. The tensor 
$M_{ijk}$ is given in (4.7) in terms of three coefficients 
$(\lambda _1, \lambda _3, \lambda _4)$ that are only functions of the spheroid aspect ratio parameter 
$A_R$. These coefficients are estimated numerically using the reciprocal theorem (see § 3 and figure 6).
In this section we investigate the spheroid's dynamics for some special cases and discuss the physics that arise. Details are below. Please note that all the results and examples are valid only up to the first order in 
$\beta$ deviation from the constant viscosity, unless otherwise mentioned.
5.1. Viscosity gradient is along or opposite the force direction
5.1.1. Governing equations
Let us examine the situation in figure 2 where the external force is in the positive 
$z$ direction, and the viscosity gradient is either parallel to the force (positive 
$z$ direction) or anti-parallel to the force (negative 
$z$ direction). In this case, the particle orientation only has one degree of freedom, namely the polar angle 
$\alpha$ measured from the 
$z$ axis. Without loss of generality, we state that 
$\boldsymbol {p}$ lies in the 
$x$–
$z$ plane, and thus, 
$\boldsymbol {p} = [\sin \alpha, 0, \cos \alpha ]$. From our theory ((4.7), (4.8b) and (4.9)), the orientation angle obeys the equation
$$\begin{gather} \frac{{\rm d} \alpha}{{\rm d}t} ={\pm} \frac{1}{2} \beta (\lambda_3 + \lambda_4) \sin (2\alpha), \end{gather}$$
$$\begin{gather}\Rightarrow |\tan (\alpha)|=|\tan (\alpha_0)|\ {\rm e}^{\beta(\lambda_3+\lambda_4)(t-t_0)}, \end{gather}$$
where 
$\pm$ illustrates the cases where the viscosity gradient is parallel (
$+$) or anti-parallel (
$-$) to the force and 
$\alpha _0$ is the initial angle at 
$t=t_0$. The translational motion of the particle obeys
\begin{equation} \frac{{\rm d} x}{{\rm d}t} = \frac{1}{2}(c_2- c_1) \sin(2\alpha), \quad \frac{{\rm d}z}{{\rm d}t} = c_1 \sin^2 \alpha + c_2 \cos^2 \alpha, \end{equation}
where 
$c_1$ and 
$c_2$ are the mobility coefficients for spheroid translation in Stokes flow (given in Appendix B in dimensional form). Major conclusions are given below.
5.1.2. Particle takes a stable orientation depending on its shape and viscosity gradient direction
Figure 7 plots the evolution of 
$\alpha$ with respect to time for prolate and oblate spheroids, for the cases when the force 
$\boldsymbol {F}$ and viscosity gradient 
$\boldsymbol {\nabla } \eta$ are parallel and anti-parallel to each other. For each set of conditions, two curves are given – one arising from the symmetry-based theory (dashed curve (5.1a)), and another from full numerical simulations where the reciprocal theorem is used at every time step (solid curve). The overlap is indistinguishable, thereby validating our theory. The second observation we make is that irrespective of the initial orientation and viscosity gradient direction (parallel or anti-parallel), both prolate and oblate spheroids evolve to a steady configuration of 
$\alpha$. This observation is very different than what is observed in Stokes flow where the orientation stays at its initial angle at all times (Leal Reference Leal2007).
Figure 7. Orientation angle 
$\alpha$ versus time for prolate and oblate spheroids when the external force 
$\boldsymbol {F}$ and viscosity gradient 
$\boldsymbol {\nabla } \eta$ are parallel or anti-parallel to each other. The left figures (a,c) correspond to prolate spheroids with 
$A_R =5$, while the right figures (b,d) correspond to oblate spheroids with 
$A_R =1/5$. The top row (a,b) is the case when the 
$\boldsymbol {F}$ and 
$\boldsymbol {\nabla }\eta$ are in the same direction, while the bottom row (b,d) is the case when they are in opposite directions. The solid curves are from full numerical simulations based on the reciprocal theorem, while the dashed curves are from the symmetry-based theory (solving (5.1a)). The dimensionless viscosity gradient is 
$\beta = 0.1$.
Figure 8 summarizes the steady orientations observed for different particle shapes and viscosity gradient directions. When the external force and the viscosity gradient are parallel to each other, the prolate spheroid adopts a stable configuration where the projector is perpendicular to the external force, while the oblate spheroid orients itself such that the projector is along the same direction as the external force. In both of these cases, the spheroid (whether prolate or oblate) has its shortest axis oriented along the direction of the viscosity gradient. On the other hand, when the spheroid is falling in the direction of decreasing viscosity (i.e. 
$\boldsymbol {F}$ and 
$\boldsymbol {\nabla }\eta$ are anti-parallel), the prolate spheroid attains a stable configuration where the projector is oriented along the force direction, whilst the oblate spheroid orients the projector perpendicular to the force direction. In both these cases, the longest axis of the particle (whether prolate and oblate) will be along the force direction.
Figure 8. Steady configurations attained by (a) prolate and (b) oblate spheroids when the external force 
$\boldsymbol {F}$ and viscosity gradient 
$\boldsymbol {\nabla }\eta$ are co-linear. The top row is for the case when the external force and the viscosity gradient are in the same direction, while the bottom row is when they are in the opposite direction.
To provide a physical understanding of this behaviour, we refer to figure 9. Here, as observed from the reference frame of the particle, the flow around the prolate spheroid bifurcates into two parts about the stagnation point and engenders both a clockwise and counterclockwise hydrodynamic torque. Figure 9(a) illustrates the magnitude of the torques for the case when the viscosity gradient is in the same direction as the force, while figure 9(b) illustrates the case when the viscosity gradient is in the opposite direction. The pictures illustrate that the unequal torques push the particle toward the stable orientations discussed above.
Figure 9. Illustration of unequal torques created on a prolate spheroid when the force and viscosity gradient are co-linear. The left figure (a) is when the viscosity gradient and force are in the same direction, while the right figure (b) is when they are in opposite directions. This schematic is shown in the particle's frame of reference.
Lastly, we note that figure 8 summarizes the unstable, steady orientations that can occur for different combinations of viscosity gradient and particle shape. These orientations only exist if the initial condition is at a specific angle, and can only be observed in exceptionally rare cases.
5.1.3. Particle translation is different than in Stokes flow
Beyond orientational kinematics, we are also interested in the translation of the spheroid. In a constant viscosity fluid with zero inertia, it is well known that the particle stays at its initial orientation (Leal Reference Leal2007). If the initial angle is 
$\alpha = 0$, 
${{\rm \pi} }/{2}$, or 
${\rm \pi}$, the particle will sediment vertically, while if 
$\alpha$ is not these values, the particle will drift in a straight, diagonal path. The direction in which the particle sediments is dictated by the resistances parallel and perpendicular to its orientation vector 
$\boldsymbol {p}$.
When a viscosity gradient is present, the translational velocity 
$\boldsymbol {U}$ obeys the same differential equation as the Stokes flow case 5.2), since we found that the viscosity gradient does not alter the force/translation coupling (see (4.8)). Thus, on the surface, it appears that the particle trajectory may seem unchanged due to the presence of a spatially varying viscosity field. However, upon closer inspection, we see that the differential equation, (4.8), depends on the particle's orientation angle 
$\alpha$, which itself is altered due to the viscosity gradient as discussed in the previous section. Thus, the viscosity gradient plays an indirect role in altering the translational dynamics.
Figure 10 plots the trajectories of oblate and prolate spheroids for different values of the initial orientation angle 
$\alpha _{0}$. For 
$\alpha _0 \neq 0$, 
${{\rm \pi} }/{2}$ and 
${\rm \pi}$, we observe motion in the sedimentation direction (3-direction) as well as a cross-stream drift (1-direction). For the case when no viscosity gradient is present, the particle moves in a straight, diagonal path. When a viscosity gradient is present, the trajectory is no longer a straight line. The cross-stream drift eventually stops when the spheroid reaches a stable orientation, beyond which the spheroid sediments vertically in the 
$3$-direction. Since the spheroid ceases to drift once the stable orientation is reached, a spheroid whose initial orientation is further away from its stable orientation will drift further than a spheroid whose initial orientation is closer to its stable orientation. Therefore, for a prolate spheroid, a particle with initial orientation 
$\alpha _0 ={\rm \pi} /4$ will drift further than one with 
$\alpha _0 = {\rm \pi}/3$, since the stable orientation is 
$\alpha ={\rm \pi} /2$ (see figure 10a). Conversely, for an oblate spheroid, the particle with an initial orientation 
$\alpha _{0} ={\rm \pi} /3$ will drift further than one with 
$\alpha _{0} ={\rm \pi} /4$, since the stable orientation is at 
$\alpha =0$.
Figure 10. Particle trajectories for (a) prolate and (b) oblate spheroids when the external force and viscosity gradient are in the same direction (
$\boldsymbol {F} = \hat {\boldsymbol {z}}, \boldsymbol {\nabla }\eta = \beta \hat {\boldsymbol {z}}$). The dashed curves correspond to when no viscosity gradient is present (
$\beta = 0$), while the solid curve is when a viscosity gradient is present (
$\beta = 0.1$). Different colours and symbols correspond to different initial starting angles 
$\alpha _0$. The prolate spheroid has 
$A_R = 5$ while the oblate spheroid has 
$A_R = 1/5$.
5.2. Viscosity gradient is perpendicular to the external force
5.2.1. Governing equations
We now examine the situation in figure 3 where the external force is in the positive 
$z$ direction (
$\boldsymbol {F} = \hat {\boldsymbol {z}}$), and the viscosity gradient is perpendicular to the force (
$\boldsymbol {\nabla }\eta = \beta \hat {\boldsymbol {x}}$). The spheroid's orientation can point in any direction, and we state it takes the form 
$\boldsymbol {p} = [\sin \alpha \cos \phi, \sin \alpha \sin \phi, \cos \alpha ]$, where 
$\alpha$ and 
$\phi$ are the polar and azimuthal angles, respectively. From our theory ((4.7), (4.8b) and (4.9)), the orientation angles evolve as
$$\begin{gather} \frac{{\rm d} \alpha}{{\rm d}t} =-\beta (\lambda_1 -\lambda_3 \sin^2\alpha + \lambda_4 \cos^2\alpha ) \cos\phi, \end{gather}$$
$$\begin{gather}\frac{{\rm d} \phi}{{\rm d} t} = \beta (\lambda_1 + \lambda_4) \cot \alpha \sin\phi, \end{gather}$$
where 
$\lambda _1$, 
$\lambda _3$ and 
$\lambda _4$ are the force/rotation mobility coefficients determined in § 4. The translation of the particle obeys the following:
\begin{equation} \left.\begin{array}{c@{}} \dfrac{{\rm d} x}{{\rm d}t} = \dfrac{1}{2} (c_2- c_1) \sin(2\alpha) \cos\phi,\\ \dfrac{{\rm d}\kern0.05em y}{{\rm d}t} = \dfrac{1}{2} (c_2- c_1) \sin(2\alpha) \sin\phi,\\ \dfrac{{\rm d}z}{{\rm d}t} = c_1 \sin^2 \alpha + c_2 \cos^2 \alpha. \end{array}\right\} \end{equation}
Here 
$c_1$ and 
$c_2$ are the mobility coefficients for particle translation in Stokes flow (given in Appendix B in dimensional form). Major conclusions are given below.
5.2.2. Particle can take a steady orientation different than the force and viscosity gradient directions
We first discuss the case when the particle starts in the same plane as 
$\boldsymbol {F}$ and 
$\boldsymbol {\nabla } \eta$ – in other words, 
$\phi _0 = 0$. From (5.3b), we see that 
$\textrm {d} \phi /\textrm {d}t = 0$ for this angle, so the particle stays at 
$\phi = 0$ and only the polar angle 
$\alpha$ will change. Figure 11 plots 
$\alpha$ versus time for both prolate and oblate spheroids, for the specific case of 
$A_R = 5$ and 
$A_R = -1/5$, respectively. First of all, we note that the results from the symmetry-based theory (solid curve, (5.3)) are virtually indistinguishable from the full numerical simulation (dashed curve), indicating the validity of our theory. Secondly, for all starting conditions, we observe the particle converges to one steady orientation. However, this steady orientation is not 
$\alpha =0$, 
$\alpha = {\rm \pi}$ or 
$\alpha ={\rm \pi} /2$, which was the case when the force and viscosity gradient vectors were co-linear.
Figure 11. Orientation angle 
$\alpha$ versus time for (a) prolate (
$A_R =5$) and (b) oblate (
$A_R =1/5$) spheroids when the external force and viscosity gradient are perpendicular (
$\boldsymbol {F} = \hat {\boldsymbol {z}}, \boldsymbol {\nabla }\eta = \beta \hat {\boldsymbol {x}}$). The dimensionless viscosity gradient is 
$\beta =0.1$, and the particle initially starts in the plane of 
$\boldsymbol {F}$ and 
$\boldsymbol {\nabla }\eta$ (i.e. 
$\phi _0 = 0$). Solid curves are from full numerical simulations based on the reciprocal theorem, while the dashed curves are from the symmetry-based theory (solving (5.3)).
We elucidate this point more clearly in figure 12. Here, we observe that neither 
$\alpha =0,{\rm \pi} /2$ nor 
${\rm \pi}$ are steady configurations because the counterclockwise torque is different than the clockwise torque at these specific angles. Some general trends are described below for prolate and oblate particles.
Figure 12. Schematic explaining the absence of steady orientations at 
$\alpha =0$ and 
$\alpha ={\rm \pi} /2$ for (a) prolate and (b) oblate spheroids when the external force and viscosity gradient are perpendicular. This schematic is shown in the particle's frame of reference.
(i) For prolate spheroids, we observe from figure 12 that the difference between the counterclockwise and clockwise torques is smaller for
$\alpha =0$ and ${\rm \pi}$ (where the long axis is along the force direction) compared with $\alpha ={\rm \pi} /2$ (where the long axis is along the viscosity gradient direction). Therefore, the steady orientation is closer to $\alpha =0$ and ${\rm \pi}$ than to $\alpha ={\rm \pi} /2$, and continues to approach $\alpha =0$ or ${\rm \pi}$ as the aspect ratio increases. In the limiting case of needle-like particles where $A_R \to \infty$, the steady orientation reaches $\alpha =n{\rm \pi}$. Between the two configurations of $\alpha =0+\varDelta$ and $\alpha ={\rm \pi} -\varDelta$ (where $\varDelta$ is a positive constant depending on aspect ratio), $\alpha ={\rm \pi} -\varDelta$ is the stable configuration, while $\alpha =0+\varDelta$ is unstable (see figure 11a).(ii) For oblate spheroids, the difference in hydrodynamic torques is larger at
$\alpha =0$ compared with $\alpha ={\rm \pi} /2$, because in the former case the longer axis is oriented along the viscosity gradient direction. Therefore, for oblate spheroids, the equilibrium orientation configuration is closer to $\alpha = {\rm \pi}/2$ than to $\alpha = 0$. In the limiting case of a thin disc where $A_R \rightarrow 0$, the stable orientation is at $\alpha ={\rm \pi} /2$. Between the two configurations of $\alpha ={\rm \pi} /2 \pm \xi$ (where $\xi$ is a positive constant depending on the aspect ratio), $\alpha = {\rm \pi}/2 -\xi$ is the stable orientation, while $\alpha = {\rm \pi}/2 +\xi$ is an unstable orientation (see figure 11b).
The results discussed above illustrate the dynamics when the initial particle orientation is co-planar with 
$\boldsymbol {F}$ and 
$\boldsymbol {\nabla } \eta$, i.e. 
$\phi _0 = 0$ or 
$\phi = {\rm \pi}$. Figure 13 plots the orientation angles 
$\phi$ and 
$\alpha$ over time when the starting angle is no longer co-planar with 
$\boldsymbol {F}$ and 
$\boldsymbol {\nabla } \eta$, i.e. 
$\phi _0 \neq 0$ or 
${\rm \pi}$. We see that at long times, the angle 
$\phi \rightarrow 0$ or 
${\rm \pi}$, i.e. the orientation ends up in the same plane as 
$\boldsymbol {F}$ and 
$\boldsymbol {\nabla } \eta$. The angle 
$\alpha$ also converges to the same result as before. Thus, we conclude that the steady orientation angles discussed previously are stable to out-of-plane perturbations.
Figure 13. Orientation angles 
$\alpha (t)$ and 
$\phi (t)$ for (a) prolate and (b) oblate spheroids when the external force and viscosity gradient are perpendicular (
$\boldsymbol {F} = \hat {\boldsymbol {z}}, \boldsymbol {\nabla }\eta = \beta \hat {\boldsymbol {x}}$). The dashed curves show the evolution of 
$\phi$, while the solid curves show the evolution of 
$\alpha$. For all cases, the initial orientation is given by the ordered pair 
$(\phi _0,\alpha _0 )=({\rm \pi} /3,{\rm \pi} /4)$ and the dimensionless viscosity gradient is 
$\beta =0.1$. The results show that 
$\phi \rightarrow 0$ or 
${\rm \pi}$, and hence, the particle becomes co-planar with 
$\boldsymbol {F}$ and 
$\boldsymbol {\nabla }\eta$.
5.2.3. Stable orientation for different aspect ratios
Figure 14 plots the steady orientation angles for prolate and oblate spheroids for different aspect ratio parameters. The steady orientations occur when 
${\textrm {d} \alpha }/{\textrm {d}t} = 0$ and 
${\textrm {d} \phi }/{\textrm {d} t} = 0$ in (5.3), which corresponds to the criterion
$$\begin{gather} \phi = n{\rm \pi}, \quad \lambda_1 - \lambda_3 \sin^2\alpha + \lambda_4 \cos^2 \alpha = 0 \Rightarrow \sin^2{\alpha} =\frac{\lambda_4+\lambda_1}{\lambda_3+\lambda_4}, \end{gather}$$
$$\begin{gather}0<\frac{\lambda_4+\lambda_1}{\lambda_3+\lambda_4} <1. \end{gather}$$
Here (
$\lambda _1, \lambda _3,\lambda _4$) are the mobility coefficients for force/rotation coupling that were calculated in § 4, which are only functions of the aspect ratio parameter 
$A_R$. Figure 14 shows that for a wide range of 
$A_R$, the above criterion is satisfied and a steady angle 
$\alpha _{se}$ exists. The stable orientation 
$\alpha _{se}$ for a prolate spheroid is closer to 
${\rm \pi}$ compared with an oblate spheroid, while the oblate spheroid has a stable angle closer to 
${\rm \pi} /2$.
Figure 14. Stable orientations 
$\alpha _{se}$ for (a) prolate and (b) oblate spheroids of different aspect ratio parameters 
$A_R$ when the external force and viscosity gradient are perpendicular to each other (
$\boldsymbol {F} = \hat {\boldsymbol {z}}$, 
$\boldsymbol {\nabla }\eta = \beta \hat {\boldsymbol {x}}$, 
$\beta = 0.1$).
5.2.4. Translation dynamics
Figure 15 shows the spheroid's translation trajectories for the case when the force and viscosity gradient are perpendicular to each other for both prolate and oblate spheroids. For both these spheroids, when the particle has a stable orientation (
$A_R =5$), the particle at long times will move in a straight, diagonal line, i.e. sediment downwards and also have a component along the viscosity gradient direction. This diagonal motion qualitatively looks similar to the motion when the spheroid is in a constant viscosity fluid (Leal Reference Leal2007). However, in a constant viscosity fluid, the angle of motion is determined by the particle's initial angle, whereas in this case, all particles will eventually move with the same trajectory, regardless of starting angle.
Figure 15. Particle trajectories for spheroids with (a) 
$A_R =5$ and (b) 
$1/A_R =5$ when the external force and viscosity gradient are perpendicular (
$\boldsymbol {F} = \hat {\boldsymbol {z}}, \boldsymbol {\nabla }\eta = \beta \hat {\boldsymbol {x}}$). The dashed curves correspond to when no viscosity gradient is present (
$\beta = 0$), while the solid curve is when a viscosity gradient is present (
$\beta = 0.1$). Different colour curves correspond to different initial starting angles 
$\alpha _0$.
5.3. General case: general direction for viscosity gradient
5.3.1. Governing equations
We now consider the most general case where 
$\boldsymbol {F}$ and 
$\boldsymbol {\nabla }\boldsymbol {\eta }$ are neither parallel or orthogonal to each other, but are inclined at an angle 
$\theta$ to each other. The external force points in the positive 
$z$ direction 
$\boldsymbol {F} = \hat {\boldsymbol {z}}$, while the viscosity gradient is
\begin{equation} \boldsymbol{\nabla}\eta = \beta \hat{\boldsymbol{d}} = \beta \cos \theta \hat{\boldsymbol{z}} + \beta \sin \theta \hat{\boldsymbol{x}}. \end{equation} Similar to before, the orientation vector is 
$\boldsymbol {p} = [\sin \alpha \cos \phi, \sin \alpha \cos \phi, \cos \alpha ]$, where 
$\alpha$ and 
$\phi$ are the polar and azimuthal angles. To determine how these angles evolve over time, we note that the dynamics is a linear superposition of the cases described previously. In other words,
$$\begin{gather} \frac{{\rm d} \alpha}{{\rm d} t} = \left.\frac{{\rm d} \alpha}{{\rm d} t}\right|_{{\parallel}} \cos \theta + \left.\frac{{\rm d} \alpha}{{\rm d} t}\right|_{{\perp}} \sin \theta, \end{gather}$$
$$\begin{gather}\frac{{\rm d} \phi}{{\rm d} t} = \left.\frac{{\rm d} \phi}{{\rm d} t}\right|_{{\parallel}} \cos \theta + \left.\frac{{\rm d} \phi}{{\rm d} t}\right|_{{\perp}} \sin \theta, \end{gather}$$
where 
$({\textrm {d} \alpha }/{\textrm {d} t})_{\parallel }$ and 
$({\textrm {d} \alpha }/{\textrm {d} t})_{\perp }$ are the variations in the polar angle from viscosity gradients parallel and perpendicular to the external force, given by (5.1a) (using the positive sign) and (5.3a), respectively. The corresponding terms 
$({\textrm {d} \phi }/{\textrm {d} t})_{\parallel }$ and 
$({\textrm {d} \phi }/{\textrm {d} t})_{\perp }$ are the same quantities for the azimuthal angle, which is zero for 
$({\textrm {d} \phi }/{\textrm {d} t})_{\parallel }$ and (5.3b) for 
$({\textrm {d} \phi }/{\textrm {d} t})_{\perp }$. The equation for particle translation is the same as (5.4).
5.3.2. Steady orientation angles
If a steady orientation angle exists, it will be in the plane spanned by 
$\boldsymbol {F}$ and 
$\boldsymbol {\nabla } \eta$ as discussed previously, i.e. 
$\phi = 0$. We set 
$\phi = 0$ and determine the conditions under which 
${\textrm {d} \alpha }/{\textrm {d}t} = 0$ in (5.7a). The criterion for a steady orientation angle is
\begin{equation} \tfrac{1}{2} (\lambda_3+\lambda_4) \sin(2\alpha)\cos{\theta}-(\lambda_1-\lambda_3\sin^2\alpha+\lambda_4\cos^2{\alpha})\sin{\theta} = 0. \end{equation} For illustration, figure 16 plots the steady orientation angles 
$\alpha _{se}$ for different values of the angle 
$\theta$ between the external force 
$\boldsymbol {F}$ and 
$\boldsymbol {\nabla }\boldsymbol {\eta }$. We observe that 
$\alpha _{se}$ varies between 
${\rm \pi} /2$ and 
${\rm \pi}$ for prolate spheroids and between 
$0$ and 
${\rm \pi} /2$ for oblate spheroids. As discussed in the previous sections, these limits are the stable orientations for very high aspect ratio spheroids when the viscosity gradients are parallel and perpendicular to the external force. For example, as 
$\theta \to 0$, we see 
$\alpha _{se} \rightarrow {\rm \pi}/2$ for prolate spheroids and 
$0$ for oblate spheroids, which are the stable orientation for these particles when the viscosity gradient is parallel to the external force.
Figure 16. Stable orientation angles 
$\alpha _{se}$ for (a) prolate and (b) oblate spheroids when the viscosity gradient 
$\boldsymbol {\nabla }\eta$ and the external force 
$\boldsymbol {F}$ are inclined at an angle 
$\theta$ to each other.
5.4. Discussion of applicability of model and incorporating disturbance viscosity
In this paper we assumed the viscosity field around the particle is a linear function of space and is independent of the flow and the particle geometry. In reality, however, the viscosity field has a more complicated spatial dependence, as it is linked to a scalar field like temperature or concentration that depends on the aforementioned quantities (see Oppenheimer et al. Reference Oppenheimer, Navardi and Stone2016 for the related problem of a hot sphere in a temperature coupled viscosity field). In this section we make suggestions on how to incorporate these effects into the analysis and what changes can be expected to the main results.
For illustrative purposes, let us consider a particle in a fluid subject to a temperature gradient 
$\boldsymbol {\nabla } T$ far away from the particle. The fluid's viscosity depends linearly on temperature, i.e. 
$\eta - \eta _0 = ({\textit {d} \eta }/{\textit {d}T}) (T - T_0)$, and thus, the viscosity field also varies spatially around the particle. If the thermal Péclet number is small and the temperature profile is steady, the temperature field will satisfy Laplace's equation inside and outside the particle (see Dassios Reference Dassios2012 for details):
\begin{equation} \nabla^2 T^{out} =0, \quad \nabla^2 T^{in} =0. \end{equation} This equation is subject to the following boundary conditions: (a) 
$T^{out} \rightarrow T_0 + \boldsymbol {\nabla } T \boldsymbol{\cdot} \boldsymbol {x}$ far away from the particle (
$|\boldsymbol {x}| \rightarrow \infty )$, and (b) on the particle surface, the temperatures and fluxes are continuous, i.e. 
$T^{in} = T^{out}$ and 
$(\boldsymbol {n} \boldsymbol{\cdot} \boldsymbol {\nabla } T^{out} ) = k_r ( \boldsymbol {n} \boldsymbol{\cdot} \boldsymbol {\nabla } T^{in} )$, where 
$k_r$ is the conductivity ratio between the particle and fluid phase. Once one solves the temperature profile, one can obtain the viscosity field 
$\eta (\boldsymbol {x})$ and then solve for the particle motion in this field. The rigid body motion will still follow the same procedure discussed earlier in the paper, i.e. one performs a perturbation expansion for the viscosity and finds the correction to the rigid body motion via the reciprocal theorem using an extra stress tensor 
$\tau _{ij}^{ex} = (\eta (\boldsymbol {x}) - \eta _0) \gamma _{ij}^{(0)}$. The mobility tensors described in § 4 will take the same form, except the numerical values for the force/rotation mobility coefficients 
$(\lambda _1, \lambda _3, \lambda _4)$ will be different. For the special cases when one neglects the presence of the particle in the transport equation, or if the conductivity ratio is 
$k_r = 1$, the viscosity field will be linear everywhere, and we recover the results described earlier in the paper. Otherwise, and in general, the viscosity field will have a more complicated spatial and even temporal dependence, emanating from the motion of the particle. Such a viscosity field, which is coupled with the particle motion, posits at least a theoretical possibility of a more complicated orientational dynamics than that presented in this paper. However, based on the work carried out previously by Shaik & Elfring (Reference Shaik and Elfring2021) for spheres, we speculate that such a motion-coupled viscosity field will not introduce any novel effects but only bolster the trends in orientational dynamics discussed in our paper. Of course, this hypothesis will need to be tested, via a full fledged analysis. The outlines of such an analysis have been presented in this subsection already. A key mathematical ingredient of such an analysis, if taken up later in a different paper, will include the solution to Laplace's equation around an ellipsoidal particle. In Appendix E we outline how one solves Laplace's equation around an ellipsoidal particle.
6. Conclusion
In this paper we study a spheroid sedimenting in a Newtonian fluid with a viscosity field that varies linearly in space. We employ the principles of linearity, reversibility, symmetry to delineate the mobility relationships for this problem. In the limit of small viscosity gradients, we find that the force/velocity and torque/rotation couplings remain unchanged from the Stokes flow limit. However, the viscosity gradient gives rise to an additional force/rotation and torque/velocity coupling, which is characterized by a third-order tensor 
$M_{ijk}$. The reduced analytical form of this tensor is given by (4.7), up to three undetermined coefficients. The values of these coefficients are determined numerically, under the aegis of a reciprocal theorem-based simulation, for a wide range of particle aspect ratios.
Illustrative examples and specific results of our theory are discussed next. Unlike in Stokes flow where the particle orientation stays at its initial orientation during sedimentation, we find that viscosity gradients alter the orientation over time. When the viscosity gradient is along the external force direction, both prolate and oblate spheroids reach a stable orientation where the longest axis is perpendicular to the viscosity gradient. When the viscosity gradient is opposite to the external force, the spheroids reach a stable orientation where the longest axis is along the viscosity gradient. We also show that for most initial orientation angles, the spheroid acquires a drift in a direction transverse to its (main) sedimentation direction until its orientation stabilizes, at which point it moves downward.
When the viscosity gradient and the external force are perpendicular, the plane defined by the viscosity gradient and force is a plane of stability, and the spheroid, irrespective of its initial orientation, will eventually become co-planar with the force and the viscosity gradient. For the limiting case of a needle-like particle, the prolate spheroid will orient its projector in the direction of the force, while conversely, for the limiting case of a flat disc, the oblate spheroid will orient its projector in the direction of the viscosity gradient. Finally, we note in the general case when the viscosity gradient and external force are neither parallel or perpendicular to each other, the dynamics of the particle is a linear combination of the cases discussed above.
Throughout the analysis, we have neglected the coupling between the viscosity field and the flow or particle motion. Guidelines for incorporating this coupling are presented, using ellipsoidal harmonics to solve the Laplace equation in a low Péclet limit. However, based on previous literature (Shaik & Elfring Reference Shaik and Elfring2021) we believe that such an analysis may not yield any novel results not yet accounted for.
We have taken a perturbative approach to solving this problem, where only the terms first order in viscosity gradient are accounted for. Due to the limitation of the perturbative approach, the analysis here is valid for the following two cases: (a) when the viscosity gradient is weak, so that the terms of 
$O(\beta ^2)$ can be neglected; (b) when the viscosity gradient is linear, but not necessarily weak, so that the terms of 
$O(\beta ^2)$ are identically zero. For this second case, where the viscosity gradient is linear, but not weak, we expect that the steady state behaviour – namely the stable orientation of the spheroids – will remain unchanged. Stronger, but still linear, viscosity gradients will change the rate at which the stable orientation is attained, but not the value of the steady orientation per se. Additionally, a spheroid is a typical axisymmetric particle with no isotropy but fore-aft symmetry. The mechanisms of hydrodynamic torque acting on the spheroid and the resultant stable orientation as discussed in this paper may therefore prove useful in understanding the dynamics of other orientable particles, which are common in nature and in the industry. The current problem may also be a stepping stone towards the analysis of more complex systems, for instance, flows with linear and quadratic components, or with density (in addition to viscosity) stratification, among others. Lastly, we have shown the existence of torque/translation coupling, in addition to the force/rotation coupling, due to the linear viscosity gradient. This coupling opens the exciting possibility of controlling the motion (translation) of the spheroid by subjecting it to external torque, say via a rotating magnetic field (Morozov & Leshansky Reference Morozov and Leshansky2014; Li et al. Reference Li, Li, Morozov, Wu, Xu, Rozen, Leshansky, Li and Wang2017b).
Acknowledgements
The authors would like to thank Prof. Gwynn Elfring and his team (specifically Dr Vaseem Shaik and Mr Jiahao Gong) at the University of British Columbia for constructive feedback and discussion on this paper.
Funding
The authors would like to acknowledge support from the American Chemical Society Petroleum Research Fund (DNI-ACS PRF 61266-DNI9), as well as support from the Michael and Carolyn Ott Endowment at Purdue University.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Disturbance velocity for an ellipsoid in Stokes flow
Consider a reference frame at an ellpsoid's centre of mass with axes aligned along the particle's principle axes. From Kim & Karilla (Reference Kim and Karilla2005, p. 55), the Stokes velocity field around the ellipsoid from external force and torque is
$$\begin{gather} v_{i} = \frac{1}{16{\rm \pi}\eta_0} \sum_{j=1}^3 F_j \left[ \delta_{ij} G_0 - x_j \frac{\partial G_0}{\partial x_i} + \frac{a_j^2}{2} \frac{\partial^2 G_1}{\partial x_i \partial x_j}\right]\!, \end{gather}$$
$$\begin{gather}v_{i} = \frac{3}{64{\rm \pi}\eta_0} \sum_{j=1}^3 \left(\boldsymbol{T}\times \boldsymbol{\nabla}\right)_j \left[ \delta_{ij} G_1 - x_j \frac{\partial G_1}{\partial x_i} + \frac{a_j^2}{4} \frac{\partial^2 G_2}{\partial x_i \partial x_j}\right]\!. \end{gather}$$
In these formulas, no summation is assumed for repeated indices unless explicitly stated. To obtain formulas for 
$v_{ki}^{trans}$ and 
$v_{ki}^{rot}$ in the reciprocal theorem, we substitute into (A1) the force and torque that comes from unit translation and rotation, respectively.
In the above expressions, the expression for 
$G_n$ is
\begin{equation} G_n(x,y,z) = \int_{\lambda}^{\infty} \left( \frac{x^2}{a^2 + t} + \frac{y^2}{b^2 + t} + \frac{z^2}{c^2 + t} - 1 \right)^n \frac{{\rm d}t}{\Delta(t)}, \end{equation}
with 
$\Delta (t) = \sqrt {(a^2 + t)(b^2 + t)(c^2 + t)}$ and 
$\lambda (x,y,z)$ being the positive root of
\begin{equation} \frac{x^2}{a^2 + t} + \frac{y^2}{b^2 + t} + \frac{z^2}{c^2 + t} = 1. \end{equation}Appendix B. Resistance formulae for an ellipsoid in Stokes flow
Let us consider a reference frame with the origin at the centre of mass of an ellipsoid and the Cartesian axes aligned along the principle axes.
We denote the ellipsoid's semi-axes as 
$(a_1,a_2,a_3) = (a,b,c)$. In dimensional form, the resistance tensors 
$R^{FU}$ and 
$R^{T\omega }$ are diagonal, while the cross-coupling term 
$R^{F\omega }=R^{TU}=0$. The diagonal elements are
\begin{equation} R_{11}^{F U}=\frac{12 \eta_0 V}{\chi_0+\alpha_1 a_1^2}, \quad R_{11}^{T \omega}=\frac{4 \eta_0 V(a_2^2+a_3^2)}{\alpha_2 a_2^2+\alpha_3 a_3^2}, \end{equation}
where 
$V=({4 {\rm \pi}}/{3}) a_1 a_2 a_3$ is the particle volume and 
$(\chi _0, \alpha _1, \alpha _2 \alpha _3)$ are elliptic integrals defined as
$$\begin{gather} \chi_0=\frac{3}{4 {\rm \pi}} V \int_0^{\infty} \frac{{\rm d} t}{\Delta(t)}, \end{gather}$$
$$\begin{gather}\alpha_i=\frac{3}{4 {\rm \pi}} V \int_0^{\infty} \frac{{\rm d} t}{(a_{\mathrm{i}}^2+\mathrm{t}) \Delta(t)}, \end{gather}$$
$$\begin{gather}\Delta(t)=\sqrt{(a_1^2+t)(a_2^2+t)(a_3^2+t)}. \end{gather}$$The other elements of the diagonal tensors are obtained by index cycling.
The mobility matrix is the inverse of the resistance matrix, and hence, given by the inverse of the diagonal elements above. For the special case when the particle is a spheroid with 
$a_1 \neq a_2=a_3$, the coefficients 
$c_1-c_4$ for the mobility matrix in (4.3a) and (4.3b) are
\begin{equation} c_1=\frac{1}{R_{22}^{F U}}, \quad c_2=\frac{1}{R_{11}^{F U}}, \quad c_3=\frac{1}{R_{22}^{T \omega}}, \quad c_4=\frac{1}{R_{11}^{T \omega}}. \end{equation}These coefficients have analytical formulae (see p. 64 and 68 in Kim and Karilla). Using the notation in this paper, we obtain the following for prolate and oblate spheroids.
(i) For prolate spheroids,
(B6a)$$\begin{gather} c_1 = \frac{1}{6{\rm \pi} \eta_0 a} \frac{1}{Y_A},\quad Y_A = \frac{16}{3} e^3 [2e + (3e^2 -1)L]^{-1}, \end{gather}$$(B6b)$$\begin{gather}c_2 = \frac{1}{6{\rm \pi} \eta_0 a} \frac{1}{X_A},\quad X_A = \frac{8}{3} e^3 [-2e + (1+e^2)L]^{-1}, \end{gather}$$(B6c)$$\begin{gather}c_3 = \frac{1}{8{\rm \pi} \eta_0 a^3} \frac{1}{Y_C},\quad Y_C = \frac{4}{3} e^3 (2-e^2) [-2e + (1+e^2)L]^{-1} , \end{gather}$$(B6d)where$$\begin{gather}c_4 = \frac{1}{8{\rm \pi} \eta_0 a^3} \frac{1}{X_C},\quad X_C = \frac{4}{3} e^3 (1-e^2) [2e - (1-e^2)L]^{-1}, \end{gather}$$$e = \sqrt {1 - {b^2}/{a^2}}$ is the spheroid's eccentricity and $L = \ln ({(1+e)}/{(1-e)})$. To get the non-dimensional form used in the paper, we multiply $c_1$ and $c_2$ by $6{\rm \pi} \eta _0 R = 6{\rm \pi} \eta _0 (ab^2)^{1/3}$, and multiply $c_3$ and $c_4$ by $6{\rm \pi} \eta _0 R^3 = 6{\rm \pi} \eta _0 ab^2$.
(ii) For oblate spheroids,
(B7a)$$\begin{gather} c_1 = \frac{1}{6{\rm \pi} \eta_0 b} \frac{1}{Y_A},\quad Y_A = \frac{8}{3} e^3 [(2e^2 + 1) C - e\sqrt{1-e^2}]^{-1}, \end{gather}$$(B7b)$$\begin{gather}c_2 = \frac{1}{6{\rm \pi} \eta_0 b} \frac{1}{X_A}, \quad X_A = \frac{4}{3} e^3 [(2e^2 - 1) C + e\sqrt{1-e^2}]^{-1}, \end{gather}$$(B7c)$$\begin{gather}c_3 = \frac{1}{8{\rm \pi} \eta_0 b^3} \frac{1}{Y_C}, \quad Y_C = \frac{2}{3} e^3 (2-e^2) [e\sqrt{1-e^2} - (1 - 2e^2) C ]^{-1} , \end{gather}$$(B7d)where$$\begin{gather}c_4 = \frac{1}{8{\rm \pi} \eta_0 b^3} \frac{1}{X_C}, \quad X_C = \frac{2}{3} e^3 [C - e\sqrt{1-e^2}]^{-1}, \end{gather}$$$e = \sqrt {1 - {a^2}/{b^2}}$ is the spheroid's eccentricity and $C = \cot ^{-1} ({\sqrt {1-e^2}}/{e} )$. To get the non-dimensional form used in the paper, we multiply $c_1$ and $c_2$ by $6{\rm \pi} \eta _0 R = 6{\rm \pi} \eta _0 (ab^2)^{1/3}$, and multiply $c_3$ and $c_4$ by $6{\rm \pi} \eta _0 R^3 = 6{\rm \pi} \eta _0 ab^2$.
Appendix C. Reciprocal theorem and $O(\beta )$ solution
To delineate the 
$O(\beta )$ correction to the particle kinematics, i.e. obtain the solution for (
$U_i^{(1)}, \omega _i^{(1)}$), there are two approaches possible. The brute force approach is to solve the velocity and stress field around the particle, and then integrate the stress on the particle's surface to find the viscosity-stratified force and torque. However, this approach is tedious and analytically intractable for complicated geometries. Instead, we circumvent the calculation of the velocity and the stress field around the particle and directly obtain the viscosity-stratified force and torque using the reciprocal theorem (Leal Reference Leal1980).
First, we note that the fluid stress field at 
$O(\beta )$ has two parts:
\begin{equation} \sigma_{ij}^{(1)} = \dot{\gamma}_{ij}^{(1)} -p^{(1)}\delta_{ij}+\tau_{ij}^{ex}. \end{equation}
One is the Newtonian part given by 
$\dot {\gamma }_{ij}^{(1)} -p^{(1)}\delta _{ij}$. The other part is viscosity stratified, denoted as 
$\tau _{ij}^{{ex}}$ and given by
\begin{equation} \tau_{ij}^{ex} = (\hat{d}_kx_k)\dot{\gamma}_{ij}^{(0)}. \end{equation}
We note the important observation that the viscosity-stratified stress at 
$O(\beta )$ depends on the strain rate at leading order.
In the spirit of the reciprocal theorem, we define an auxiliary problem wherein the same particle, at the same location and same orientation, is sedimenting in a Newtonian fluid with a constant (spatially invariant) viscosity. The quantities pertaining to the auxiliary problem are denoted by the 
${aux}$ superscript. Therefore, the external force and the torque acting on the particle in the auxiliary problem is given by 
$F_i^{aux},T_i^{aux}$ and its rigid body motion is given by 
$U_i^{aux},\omega ^{aux}_i$. The flow field around the particle is 
$v_i^{aux}$, while the stress field is 
$\sigma _{ij}^{aux}$, expressed as
\begin{equation} \sigma_{ij}^{aux} = \dot{\gamma}_{ij}^{aux} - p^{aux}\delta_{ij}. \end{equation} Since the stress field of the auxiliary problem and that of the 
$O(\beta )$ problem are divergence free:
\begin{equation} \frac{\partial \sigma_{ij}^{aux}}{\partial x_j} =\frac{\partial \sigma_{ij}^{(1)}}{\partial x_j} =0, \end{equation}or
\begin{equation} v_i^{(1)}\frac{\partial \sigma_{ij}^{aux}}{\partial x_j} =v_i^{aux}\frac{\partial \sigma_{ij}^{(1)}}{\partial x_j} =0. \end{equation}Using the product rule, the above equation reduces to
\begin{equation} \frac{\partial v_i^{(1)}\sigma_{ij}^{aux}}{\partial x_j} - \frac{\partial v_i^{aux}\sigma_{ij}^{(1)}}{\partial x_j} =\sigma_{ij}^{aux}\frac{\partial v_i^{(1)}}{\partial x_j}-\sigma_{ij}^{(1)}\frac{\partial v_i^{aux}}{\partial x_j}. \end{equation}
We now substitute the expressions 
$\sigma _{ij}^{aux} = \dot {\gamma }_{ij}^{aux} - p^{aux}\delta _{ij}$ and 
$\sigma _{ij}^{(1)} = \dot {\gamma }_{ij}^{(1)}-p^{(1)}\delta _{ij}+\tau _{ij}^{ex,(0)}$ to the right-hand side. Using the identities 
${\partial v_i}/{\partial x_i} = {\partial v_i^{aux}}/{\partial x_i} = 0$ and 
$\dot {\gamma }_{ij}^{(1)}({\partial v_i^{aux}}/{\partial x_j}) = \dot {\gamma }_{ij}^{aux}({\partial v_i^{(1)}}/{\partial x_j})$, we obtain
\begin{equation} \frac{\partial v_i^{(1)}\sigma_{ij}^{aux}}{\partial x_j} - \frac{\partial v_i^{aux}\sigma_{ij}^{(1)}}{\partial x_j} =-\tau_{ij}^{ex,(0)}\frac{\partial v_i^{aux}}{\partial x_j}. \end{equation}Next, we integrate the above equation over the volume outside the particle and use the divergence theorem. This procedure yields
\begin{equation} \int_{\mathcal{S}} n_j v_i^{(1)}\sigma_{ij}^{aux}\,{\rm d}S = \int_{\mathcal{S}}n_jv_i^{aux}\sigma_{ij}^{(1)}\,{\rm d}S + \int_{V}\tau_{ij}^{ex,(0)}\frac{\partial v_i^{aux}}{\partial x_j}\,{\rm d}V, \end{equation}
where 
$\mathcal {S}$ is the surface of the particle and 
$n_j$ is the normal to the particle surface pointing inside the fluid. On particle surface, 
$v_i^{(1)}$ and 
$v_i^{aux}$ are rigid body motion, i.e. 
$v_i^{(1)} = U_i^{(1)} + \epsilon _{ijk} \omega _j^{(1)}x_k$ and 
$v_i^{aux} = U_i^{aux} + \epsilon _{ijk} \omega _j^{aux}x_k$. Substituting these expressions into the surface integrals yield
\begin{equation} -F_i^{aux} U_i^{(1)}-T_i^{aux}\omega_i^{(1)} = \int_{V}\tau_{ij}^{ex,(0)}\frac{\partial v_i^{aux}}{\partial x_j}\,{\rm d}V. \end{equation}
Note that when deriving the above expression, we made use of the fact that the force and torque acting on the particle at 
$O(\beta )$ is zero (
$F_i^{(1)} = T_i^{(1)} = 0$). Lastly, let us write the auxillary force and torque as a linear combination of the rigid body velocities using the resistance tensors for the particle:
$$\begin{gather} F_i^{aux} = R_{ij}^{FU} U_j^{aux} + R_{ij}^{F\omega} \omega_j^{aux}, \end{gather}$$
$$\begin{gather}T_i^{aux} = R_{ij}^{TU} U_j^{aux} + R_{ij}^{T\omega} \omega_j^{aux}. \end{gather}$$
Here in the above equation, the resistance tensors satisfy the following symmetry relationships: 
$\boldsymbol {R}^{FU} = (\boldsymbol {R}^{FU})^\textrm {T}$, 
$\boldsymbol {R}^{T\omega } = (\boldsymbol {R}^{T\omega })^\textrm {T}$ and 
$\boldsymbol {R}^{F\omega } = (\boldsymbol {R}^{TU})^\textrm {T}$. We also write the auxillary velocity field in the volume integral for (C9) as a linear combination of the rigid body motions:
\begin{equation} v_i^{aux} = v_{ik}^{trans} U_k^{aux} + v_{ik}^{rot} \omega_k^{aux}. \end{equation}
Here 
$v_{ik}^{rot}$ and 
$v_{ik}^{rot}$ are the velocity fields in the 
$i$ direction induced by unit translation or rotation in the 
$k$ direction. Substituting (C10) and (C11) into (C9) and eliminating 
$U_i^{aux}$ and 
$\omega _i^{aux}$ yields the final result, (3.6), stated in the paper.
Appendix D. Simplification of mobility tensor $M_{ijk}$
Here, we show that (4.6) is equivalent to (4.7). To that end, we rewrite (4.6) as
\begin{equation} M_{ijk} =\lambda_1\underbrace{\epsilon_{ijk}}_{{Term\ 1}}+\lambda_2\underbrace{p_i\epsilon_{jkq}p_q}_{{Term\ 2}} +\lambda_3\underbrace{p_j\epsilon_{ikq}p_q}_{{Term\ 3}} +\lambda_4\underbrace{p_k\epsilon_{ijq}p_q}_{{Term\ 4}}. \end{equation}
Without any loss of generality, we assume a particular orientation of the projection vector 
$p_i$ namely 
$p_i =\delta _{i1}$ (and so on). Therefore, (D2) may be written as
\begin{equation} M_{ijk} =\lambda_1\underbrace{\epsilon_{ijk}}_{{Term\ 1}}+\lambda_2\underbrace{\delta_{i1}\epsilon_{jk1}}_{{Term\ 2}} +\lambda_3\underbrace{\delta_{j1}\epsilon_{ik1}}_{{Term\ 3}} +\lambda_4\underbrace{\delta_{k1}\epsilon_{ij1}}_{{Term\ 4}}. \end{equation}Say,
$$\begin{gather} M_{ijk}^{(1)} = \text{Term 1} =\epsilon_{ijk}, \end{gather}$$
$$\begin{gather}M_{ijk}^{(2)} = \text{Term 2} =\delta_{i1}\epsilon_{jk1}, \end{gather}$$
$$\begin{gather}M_{ijk}^{(3)} = \text{Term 3} =\delta_{j1}\epsilon_{ik1}, \end{gather}$$
$$\begin{gather}M_{ijk}^{(4)} = \text{Term 4} =\delta_{k1}\epsilon_{ij1}. \end{gather}$$ To expand the different tensors, term by term, we find that the only non-zero terms in 
$M_{ijk}^{(1)}$ are
$$\begin{gather} M_{123}^{(1)} = M_{312}^{(1)} = M_{231}^{(1)} =1, \end{gather}$$
$$\begin{gather}M_{132}^{(1)} = M_{321}^{(1)} = M_{213}^{(1)} =-1. \end{gather}$$ Similarly, the non-zero terms in 
$M_{ijk}^{(2)}$ are given as
$$\begin{gather} M_{123}^{(2)} =1, \end{gather}$$
$$\begin{gather}M_{132}^{(2)} =-1; \end{gather}$$
and the non-zero terms in 
$M_{ijk}^{(3)}$ are given as
$$\begin{gather} M_{213}^{(3)} =-1, \end{gather}$$
$$\begin{gather}M_{312}^{(3)} =1; \end{gather}$$
whilst the non-zero terms in 
$M_{ijk}^{(4)}$ are given as
$$\begin{gather} M_{231}^{(4)} = 1, \end{gather}$$
$$\begin{gather}M_{321}^{(4)} =-1. \end{gather}$$From the visual inspection of (D4), (D5), (D6), (D7), we obtain the following relationship:
\begin{equation} M_{ijk}^{(1)} = M_{ijk}^{(2)}+M_{ijk}^{(3)}+M_{ijk}^{(4)}, \end{equation}
which means out of the four terms 
$M_{ijk}^{(1)}$, 
$M_{ijk}^{(2)}$, 
$M_{ijk}^{(3)}$ and 
$M_{ijk}^{(4)}$, only three are linearly independent, and therefore, without loss of generality, we can remove 
$M_{ijk}^{(2)}$ (
$=\textrm {Term 2}$) from (D2) (or (4.6)), which leads to, with a slight change in notation, (4.7).
Appendix E. Solving Laplace equation around an ellipsoid with a far-field temperature gradient
Here we outline how to solve Laplace's equation around an ellipsoidal particle. We use ellipsoidal harmonics, a technique widely used in electrostatics, and the results in Sten (Reference Sten2006) directly apply here. Let us consider a frame of reference where the Cartesian coordinate system 
$(x,y,z)$ aligns with the semi-major axes 
$(a,b,c)$ of the ellipsoid, with 
$a \geq b \geq c$. If the far-field temperature is
\begin{equation} T^{\infty}(\boldsymbol{x}) = T_0 + \frac{\partial T}{\partial x} x + \frac{\partial T}{\partial y} y + \frac{\partial T}{\partial z} z, \end{equation}the solution outside the ellipsoid takes the following form:
\begin{equation} T(\boldsymbol{x}) = T^{\infty}(\boldsymbol{x}) + \sum_{p=1}^3 B_{1p} \mathcal{F}_1^p(\boldsymbol{x}). \end{equation}
In the above equation, 
$\mathcal {F}_1^p(\boldsymbol {x}) = F_1^p(\xi ) E_1^p(\mu ) E_1^p(\nu )$ are decaying ellipsoidal harmonics using the ellipsoidal coordinate system 
$(\xi, \mu,\nu )$, where 
$\xi = a$ denotes the surface of the ellipsoid. The functions 
$E_1^p$ and 
$F_1^p$ are Lame functions of the first and second kind, defined in Sten (Reference Sten2006). In (E2) the coefficients 
$B_{1p}$ are
$$\begin{gather} B_{11} = \frac{abc}{3} \frac{1}{kh} \frac{ (1 - k_r)}{1 + L_1^1(a) (k_r -1)} \frac{\partial T}{\partial x}, \end{gather}$$
$$\begin{gather}B_{12} = \frac{abc}{3} \frac{1}{h\sqrt{k^2 - h^2}} \frac{ (1 - k_r)}{1 + L_1^2(a) (k_r -1)} \frac{\partial T}{\partial y}, \end{gather}$$
$$\begin{gather}B_{13} = \frac{abc}{3} \frac{1}{k \sqrt{k^2 - h^2}} \frac{ (1 - k_r)}{1 + L_1^3(a) (k_r -1)} \frac{\partial T}{\partial z}, \end{gather}$$
where 
$k = \sqrt {a^2 - c^2}$ and 
$h = \sqrt {a^2 - b^2}$. The geometric factors 
$L_1^{1,2,3}(a)$ take the following form (Sten Reference Sten2006):
\begin{equation} L_1^{1,2,3}(a) = abc \int_a^{\infty} \frac{ {\rm d} \xi'}{ [ E_1^{1,2,3}(\xi')]^2 \sqrt{( \xi'^{2} - h^2 ) ( \xi'^{2} - k^2)} }. \end{equation}
Open access
