Sedimentation of spheroids in Newtonian fluids with spatially varying viscosity

Abstract This paper examines the rigid body motion of a spheroid sedimenting in a Newtonian fluid with a spatially varying viscosity field. The fluid is at zero Reynolds number, and the viscosity varies linearly in space in an arbitrary direction with respect to the external force. First, we obtain the correction to the spheroid's rigid body motion in the limit of small viscosity gradients, using a perturbation expansion combined with the reciprocal theorem. Next, we determine the general form of the particle's mobility tensor relating its rigid body motion to an external force and torque. The viscosity gradient does not alter the force/translation and torque/rotation relationships, but introduces new force/rotation and torque/translation couplings that are determined for a wide range of particle aspect ratios. Finally, we discuss results for the spheroid's rotation and centre-of-mass trajectory during sedimentation. A steady orientation arises at long time whose value depends on the viscosity gradient direction and particle shape. These results are significantly different than when no viscosity gradient is present, where the particle stays at its initial orientation for all times. We summarize the observations for prolate and oblate spheroids for different viscosity gradient directions and provide plots for the orientation and centre-of-mass trajectory versus time. We also provide guidelines to extend the analysis when the viscosity gradient exhibits a more complicated spatial behaviour.


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 (Howard Berg 2004), and therefore exhibit variable viscosity, either with (Esparza López et al. 2021) or without sharp viscosity gradients (Du et al. 2012).Similarly, gradients in temperature, salinity, or concentration may induce spatial variation in viscosity, most commonly observed in marine ecosystems (Arrigo et al. 1999).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. 2004;Rafaï et al. 2010).
In this manuscript, we will examine an idealized problem of a single spheroid sedimenting in a spatially varying viscosity field.We will discuss the dynamics that are observed, and how they differ from other situations studied in the literature.By now, it is wellknown 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 (Fig. 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 (Auguste et al. 2013;Cox 1965;Khayat & Cox 1989), or if the suspending fluid has normal stresses due to the presence of polymers (Galdi 2000;Galdi et al. 2011;Kim 1986).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 et al. 2015).Conversely, fluid viscoelasticity orients the spheroid such that its longest axis is along the force directioni.e., an "edge wise" configuration (Dabade et al. 2015;Kim 1986).These effects markedly change the particle trajectory as well as the sedimentation speed (Fig. 1), since the particle's drag coefficient is a function of orientation and is minimized when the longest axis is along the force direction.
Another way in which symmetry breaking could occur is if there is a stratified fluidi.e., variations in density, viscosity, or other fluid properties that alter the force and torque on the particle (More & Ardekani 2022).This area of research is relatively modern, and most of the efforts have examined the effect of density stratification on particle dynamics (Ardekani et al. 2017;Doostmohammadi et al. 2014).When density increases along the gravity direction, it is found that the drag on a sphere is enhanced as confirmed by theory (Mehaddi et al. 2018), experiments (Lofquist & Purtell 1984;YICK et al. 2009) and simulations (Hanazaki et al. 2009;More et al. 2021).The buoyancy force also leads to continuous deceleration and absence of a terminal velocity (Doostmohammadi et al. 2014).For anisotropic particles like spheroids, there has been some research to understand their settling behavior in density stratified fluids.Using a reciprocal theorem based approach, Varanasi and Subramanian (Varanasi & Subramanian 2022) 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 et al (Dandekar et al. 2020).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. 2021;Varanasi & Subramanian 2022).For example, density is often linked to a scalar field like temperature, which depends on a convection-diffusion equation.Depending on the Peclet number, the density around the particle may or may not be coupled with the fluid flow.In the low Peclet number limit, this additional torque is opposite the hydrostatic torque (More et al. 2021;Varanasi & Subramanian 2022).
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 (Dandekar & Ardekani 2020;Jacquemin et al. 2006).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. 2004;Liebchen et al. 2018;Rafaï et al. 2010;Sokolov & Aranson 2009), as well as roboticists who design microrobots in such fluids (Asghar et al. 2020;Kim et al. 2016;Li et al. 2017;Nelson et al. 2010;Palagi et al. 2013;Zhuang et al. 2017).Some questions that arise are: how does a spatially varying fluid viscosity affect the common swimming speed, propulsion, and efficiency (Gagnon & Arratia 2016)?Do microswimmers orient themselves in preferable positions in response to the viscosity gradients (Takabe et al. 2017)?The common approach is to leverage a prototypical swimmer model (squirmers (Datt & Elfring 2019;Shaik & Elfring 2021), swimming sheet (Dandekar & Ardekani 2020;Eastham & Shoele 2020), Purcell's swimmer (Qin & Pak 2023), cilia (Asghar et al. 2020;Palagi et al. 2013)) and then couple it to the Stokes flow field with a variable viscosity.Currently, work has been performed on the the motion of a single sphere in a viscosity varying fluid (Datt & Elfring 2019), 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 will investigate in this paper.
With this motivation in mind, this manuscript will 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).Sec. 2 outlines the particle geometry and equations of motion.Sec. 3 numerically solves for the particle's rigid body motion in the limit of weak viscosity gradient using the reciprocal theorem.Sec. 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 arbitrary set of forcing, viscosity gradient direction, and particle geometry.Sec. 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 vs. oblate spheroid), and particle aspect ratio, the spheroid can take on different steady orientation angles, and sometimes experience no steady orientation.The section concludes on how to extend the analysis to more complicated situations, followed by Sec. 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 behavior 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 (Butler & Shaqfeh 2002;Herzhaft & Guazzelli 1999;Koch & Shaqfeh 1989, 1991;Kuusela et al. 2003;Nicolai et al. 1998;Shin et al. 2006Shin et al. , 2009;;Vishnampet & Saintillan 2012).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.

Problem Geometry
The schematic of our system is shown in Figs. 2 and 3. We consider a torque-free spheroid under an external force F in a Newtonian fluid with a constant viscosity gradient ∇η.The force is in the positive 3-direction.The viscosity gradient ∇η can be co-linear with the force (Fig. 2, where ∇η is in the ±3-direction) or perpendicular to the force (Fig. 3, where ∇η is in the +1-direction).The spheroid has three semi-major axes of lengths (a, b, c), with a = b = c.The initial center of mass of the spheroid is (x 01 , x 02 , x 03 ) = (0, 0, 0).
We will define the spheroid's orientation vector p as the direction along its unequal axis (i.e., the a-axis).Two different cases arise.A prolate spheroid has p along its longest axis, while an oblate spheroid has 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.
The two systems of parameterization are connected by the following relationship: Evidently, a prolate particle has its aspect ratio parameter A R > 1, while an oblate particle has its aspect ratio parameter A R < 1. Figs. 2 and 3 describe the polar and azimuthal angles α ∈ [0, π] and φ ∈ [0, 2π) for the particle orientation.The next section discusses the equations of motion and the rheology of the fluid.

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 = (ρ f U L max )/η 0 ≈ 0. Here, ρ 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 η 0 is 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: where σ 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 following form: where p is the pressure, γij = ∂vj ∂xi + ∂vi ∂xj is twice the strain rate tensor and η is the viscosity of the medium.In this problem, the viscosity is independent of the strain rate, unlike shear thinning (Anand 2014;Anand & Christov 2019b) and viscoelastic (Anand 2016)fluids, but exhibits a spatial dependence.The viscosity field is: In the above equation, η 0 is the viscosity at the origin and ∇η = η0 R β d is a constant viscosity gradient with dimensionless magnitude β and unit direction d.
The goal of the problem is to solve Eqs.(2.3), (2.4), and (2.5) for the stress and velocity around the particle.The equations have to be solved with the following boundary conditions: where (U i , ω i ) are the rigid body velocities of the particle, S p is the particle surface, x cm k is the center of mass, and ǫ ijk is the Levi-Civita symbol.An additional constraint is that the particle's external force and torque are specified.These are: 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 small Peclet 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 manuscript (Sec.5.4).
Irrespective of the rheology of the fluid, due to the introduction of the particle, the flow around the particle changes to satisfy the no slip boundary condition on the surface of the particle.The flow , in turn, applies hydrodynamic force (and torque) on the particle, thereby affecting the translation and rotation of the particle, and if the particle is soft, also its deformation.This interaction between the fluid flow and the particle means that the current problem may also be interpreted as a fluid-structure interaction (FSI) problem .FSI problems have already been studied extensively in the case of deformable channels (Anand et al. 2019;Venkatesh et al. 2022)and tubes (Anand & Christov 2019a, 2020, 2021) conveying Newtonian and non Newtonian fluids in steady as well as transient conditions.

Non dimensionalization, dimensionless numbers and perturbation expansion
Unless otherwise noted, all quantities from here on out will be written in nondimensional form.Lengths will be scaled by the average particle size R, forces by its magnitude F , and viscosities by its value at the origin η 0 .Velocities will be scaled by the Newtonian sedimentation velocity U = F 6πη0R , times by t c = R/U , strain rates and rotational velocities by γc = 1/t c , stresses and pressures by η 0 γc , and torques (if present) by F R.
The dynamics of the spheroid will depend on the following dimensionless quantitiesthe particle aspect ratio parameter A R , the particle orientation p (characterized by angles α and φ), and the non-dimensional viscosity gradient ∇η (characterized by magnitude β and direction d): In dimensionless form, the viscosity of the fluid in Figs. 2 and 3 is the following: where the first case corresponds to the case where the viscosity gradient is parallel (+3 direction) or anti-parallel (−3 direction) to the external force, and the second case where the viscosity gradient is perpendicular to the external force.For a general viscosity gradient ∇η, 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 small viscosity gradient: The above condition indicates that one can neglect fluid inertia and perform a regular perturbation expansion in β.We will solve for the rigid body motion up to O(β), both numerically and semi-analytically using symmetry arguments listed in the next sections.

Reciprocal theorem
We will determine the rigid body motion of the spheroid by performing a perturbation expansion in the non-dimensional viscosity gradient β ≪ 1.We perturb the dependent variables as follows: i , U i , ω i , U i , ω (1) i + . . .
and solve for the momentum and mass balances Eqs.(2.3) -(2.7) at each order in β.At leading order, the spheroid sediments in a zero Reynolds number fluid with a constant, non-dimensional viscosity η = 1 and a non-dimensional external force F = 1: The solution to the above problem is given in many classical texts (for example see (Kim & Karilla 2005)).The velocity field is presented in Appendix A, while the rigid body motion satisfies the classical resistance relationship: Here the body force is due to the spatially varying viscosity field: In the above Eq.(3.5), τ ex ij is the extra stress tensor, and γ(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(β) problem.This theorem has a storied history in the Stokes flow community, as 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 Eq. (3.3), except the forces and torques are replaced by an effective polymeric force and torque: The polymeric force and torque are given as follows: where the integrals are evaluated over the volume V out outside the particle.The quantities v trans ki and v trans ki 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.

Numerical implementation:
The volume integrals in the Eq.(3.7) are difficult to evaluate analytically.A custommade 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 et al. 2020), 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 that the origin is at the particle's center of mass and the Cartesian coordinate axes align with the particle's principle axes.We then evaluate the volume integrals in Eq. (3.7) for the polymeric force and torque, using an elliptical coordinate system and performing Gaussian quadrature via Legendre polynomials.We then solve the matrix equations Eq. (3.3) and Eq.(3.6) for the rigid body motions at O(1) and O(β), and transform back to the laboratory frame.The particle's center of mass and orientation are evolved by solving the rigid body dynamics: We use a forward Euler scheme with ∆t = 0.01.More details are found in our prior publications (Anand & Narsimhan 2023;Wang et al. 2020).

Verification of code
For the case of a sphere sedimenting in a linear, imposed viscosity gradient, we refer to the work by (Datt & Elfring 2019).Specifically, Eqs.(7, 8) in (Datt & Elfring 2019) are the resistance relationships for the external force F and torque T on a sphere of radius a in a fluid with a constant viscosity gradient ∇η, with translational velocity U and rotational velocity ω.For convenience, these equations are reproduced in dimensional form here: (i) Spatial variation in y direction: For a torque-free (T = 0) sphere sedimenting in the x-direction where the dimensional viscosity gradient is along the y-direction ∇η = β η0 a ŷ, the above equations give us: The analytical result for ω z in Eq. (3.10) is compared against the results of the numerical simulation and the comparison is shown in Fig. 4(a) showing an accurate match.
(ii) Spatial variation in x direction: Similarly, for a torque-free sphere sedimenting in the x-direction where the dimensional viscosity gradient is along the x-direction ∇η = β η0 a x, the above equations give: We compare the results of the simulation against Eq.(3.11 where a good match is seen.

Introduction and motivation
The simulations described in the previous section solve the rigid body motion of the particle, but are computationally intensive.At each timestep, one has to evaluate six volume integrals in Eq. (3.7) to obtain the polymeric 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 below.

General form of mobility tensor
The governing momentum and continuity Eqns.(2.3) -(2.7) are linear in the external force and torque (F , T ).Thus, the translational and rotational velocities (U , ω) are also linear in these quantities and obey the following relationship: Here, (A, B, D) are mobility tensors that are non-dimensionalized by 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 p.If a viscosity gradient is present, the tensors will also be a function of the non-dimensional viscosity gradient ∇η = β d.Note: the the off-diagonal terms of the matrix in Eq. (4.1) are transposes of each other as can be proved by the reciprocal theorem (not shown here).
In the limit of β ≪ 1, we expand the mobility tensors in a Taylor series as follows: 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 p, they take the form: where c 1 , c 2 , c 3 and c 4 are functions of the apsect ratio parameter and are given in Appendix B.
At O(β), the motion will be linear in ∇η.Thus, in non-dimensional form, the mobility tensors take the following structure: where dk is the direction of the viscosity gradient.Therefore the problem of finding the mobility matrices (A ij ) reduces to the problem of finding α ijk ,β ijk and M ijk .For a spheroid, these third order tensors depend on the orientation product p i p j , since fore-aft symmetry dictates that changing p i to −p i will not alter the results.Noting that (α ijk , β ijk ) are third order true tensors, and such tensors cannot be formed from p i p j , we obtain the result: The above relationship means that at O(β), the force-velocity coupling and torqueangular velocity coupling are unchanged.However, as we will see next, the force-rotation coupling and torque-velocity coupling will change.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: where λ 1 , λ 2 , λ 3 , λ 4 are dimensionless coefficients that depend only on the aspect ratio parameter A R .One can show that without loss of generality λ 2 = 0 (see Appendix D) and therefore the problem reduces to finding the coefficients λ 1 , λ 3 , λ 4 .In other words, Eq. (4.6) reduces to In summary, the mobility relationships up to O(β) reduce to: ) where ij are the known mobility tensors for a spheroid without a viscosity gradient, given by Eq. (4.3), while M ijk is the cross-coupling term given by Eq. (4.7).The unknown coefficients for the tensor M ijk are (λ 1 , λ 3 , λ 4 ), which are functions of the aspect ratio parameter A R for the spheroid.The next section discusses how we determine these coefficients.4.3.Determining coefficients λ 1 , λ 3 , λ 4 for the mobility matrix M ijk (force-rotation and torque-translation coupling) Fig 5 outlines the simulations we perform to obtain the coefficients (λ 1 , λ 3 , λ 4 ) for M ijk in Eq. (4.7).We examine a torque-free particle (T i = 0) and quantify its angular velocity ω 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: while the rate of spinning is: These quantities are computed for the cases below:  its orientation in the x − z plane with an angle α with respect to the force directioni.e., p = [sin α, 0, cos α].Using Eqs.(4.7), (4.8b), and (4.9), one finds the angular velocity to be: Thus, performing one simulation at a specific polar angle and viscosity gradient magnitude (say α = π/4, β = 0.1) allows us to obtain (λ 3 + λ 4 ).In Fig 6, we plot simulations of 2 dα dt sin 2α for many values of angles α and non-dimensional viscosity gradients β.This quantity is constant for all values of α and β, but is a function of the aspect ratio parameter A R , which is consistent with the expression listed above.(4.7), (4.8b), and (4.9), one finds the angular velocity to be: Performing one simulation at a specific value of β (e.g., β = 0.1) allows us to obtain (λ 1 − λ 3 + λ 4 ).
(iii) Case C: ∇η ⊥ F ⊥ p: We examine case in Fig 5c where the orientation, viscosity gradient, and force are all perpendicular to each other -i.e., F = ẑ, d = x, p = ŷ.Here, the particle will spin but not change orientation.Using Eqs.(4.7), (4.8b), and (4.10), we find the spinning rate to be: Performing one simulation at a specific value of β (e.g., β = 0.1) allows us to obtain λ 1 .
The three simulations listed above yield a linear system of equations for the coefficients (λ 1 , λ 3 , λ 4 ) that can be solved.Fig 7 shows the values of the coefficients for different values of the aspect ratio parameter A R , for both prolate and oblate spheroids.Once these coefficients are tabulated, one has a general form for the rigid body motion (Eq.(4.8)) for spheroids that can be solved for arbitrary viscosity gradient, orientation, aspect ratio, and external force/torque.

Results and illustrative examples
In Sec. 4, we developed a theory to describe the rigid body motion of a spheroid in a spatially varying viscosity field.The general form of the translational and rotational velocities is given in Eq. (4.8), where A (0) ij and D (0) ij 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 Eq. (4.7) in terms of three coefficients (λ 1 , λ 3 , λ 4 ) that are only functions of the spheroid aspect ratio parameter A R .These coefficients are estimated numerically using the reciprocal theorem (see Sec. 3 and Fig. 7).
In this section, we investigate the spheroid's dynamics for some special cases and discuss the physics that arise.Details are below.5.1.Viscosity gradient is along or opposite the force direction 5.1.1.Governing equations Let us examine the situation in Fig. 2 where the external force is in the positive zdirection, 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 α measured from the z-axis.Without loss of generality, we will state that p lies in the x − z plane, and thus p = [sin α, 0, cos α].From our theory (Eqs.(4.7), (4.8b), and (4.9)), the orientation angle obeys the following equation: where ± illustrates the cases where the viscosity gradient is parallel (+) or anti-parallel (−) to the force.The translational motion of the particle obeys: 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 Fig. 8 plots the evolution of α with respect to time for prolate and oblate spheroids, for the cases when the force F and viscosity gradient ∇η are parallel and anti-parallel to each other.For each set of conditions, two curves are given -one arising from the reduced order theory (dashed curve, Eq. (5.1)), 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 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 α.This observation is very different than what is observed in Stokes flow where the orientation stays at its initial angle at all times (Leal 2007).
Fig. 9 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 external force, while the oblate spheroid orients itself such 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., F and ∇η 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.
To provide a physical understanding of this behavior, we refer to Fig. 10.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 counter-clockwise hydrodynamic torque.Fig. 10(a) illustrates the magnitude of the torques for the case when the viscosity gradient is in the same direction as the force, while Fig. 10(b) illustrates the case when the viscosity gradient is in the opposite direction.
The pictures illustrate that the the unequal torques push the particle toward the stable orientations discussed above.Lastly, we note that Fig. 9 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.

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 2007).If the initial angle is α = 0, π 2 , or π, the particle will sediment vertically, while if α 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 p.
When a viscosity gradient is present, the translational velocity U obeys the same differential equation as the Stokes flow case (Eq.(5.2)), since we found that the viscosity gradient does not alter the force/translation coupling (see Eq. (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 (Eq.(4.8)) depends on the particle's orientation angle α, 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.The dashed curves correspond to when no viscosity gradient is present (β = 0), while the solid curve is when a viscosity gradient is present (β = 0.1).Different color curves correspond to different initial starting angles α 0 .The prolate spheroid has A R = 5 while the oblate spheroid has A R = 1/5.Fig. 11 plots the trajectories of oblate and prolate spheroids for different values of the initial orientation angle α 0 .For α 0 = 0, π 2 , and π, 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 α 0 = π/4 will drift further than one with α 0 = π/3, since the stable orientation is α = π/2 (see Fig. 11(a)).Conversely, for an oblate spheroid, the particle with an initial orientation α 0 = π/3 will drift further than one with α 0 = π/4, since the stable orientation is at α = 0.

Governing equations
We will now examine the situation in Fig. 3 where the external force is in the positive z-direction (F = ẑ), and the viscosity gradient is perpendicular to the force (∇η = β x).The spheroid's orientation can point in any direction, and we state it takes the form p = [sin α cos φ, sin α sin φ, cos α], where α and φ are the polar and azimuthal angles, respectively.From our theory (Eqs.(4.7), (4.8b), and (4.9)), the orientation angles evolve as follows:   where λ 1 , λ 3 , and λ 4 are the force-rotation mobility coefficients determined in Sec. 4. The translation of the particle obeys the following: (5.4) where 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.

Particle can take a steady orientation different than the force and viscosity gradient directions
We will first discuss the case when the particle starts in the same plane as F and ∇η -in other words φ 0 = 0. From Eq. (5.3b), we see that dφ/dt = 0 for this angle, so the particle stays at φ = 0 and only the polar angle α will change.The dashed curves show the evolution of φ, while the solid curves show the evolution of α.For all cases, the initial orientation is given by the ordered pair (φ 0 , α 0 ) = (π/3, π/4) and the dimensionless viscosity gradient is β = 0.1.The results show that φ → 0 or π, and hence the particle becomes co-planar with F and ∇η.
. A R = −1/5, respectively.First of all, we note that the results from the reduced order theory (solid curve, Eq. ( 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 α = 0, α = π, or α = π/2, which was the case when the force and viscosity gradient vectors were co-linear.
We elucidate this point more clearly in Fig. 13.Here, we observe that neither α = 0, π/2, or π are steady configurations because the counter-clockwise torque is different than the clockwise torque at these specific angles.Some general trends are described below for prolate and oblate particles: • Prolate spheroids: For prolate spheroids, we observe from Fig. 13 that the difference between the counter-clockwise and clockwise torques is smaller for α = 0 and π (where the long axis is along the force direction) compared to α = π/2 (where the long axis is along the viscosity gradient direction).Therefore, the steady orientation is closer to α = 0 and π than to α = π/2, and continues to approach α = 0 or π as the aspect ratio increases.In the limiting case of needle like particles where A R → ∞ the steady orientation reaches α = nπ.Between the two configurations of α = 0 + ∆ and α = π − ∆ (where ∆ is a positive constant depending on aspect ratio), α = π − ∆ is the stable configuration, while α = 0 + ∆ is unstable (see Fig. 12(a)).
• Oblate spheroids: For oblate spheroids, the difference in hydrodynamic torques is larger at α = 0 compared to α = π/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 α = π/2 than to α = 0.In the limiting case of a thin disc where A R → 0, the stable orientation is at α = π/2.Between the two configurations of α = π/2 ± ξ (where ξ is a positive constant depending on the aspect ratio), α = π/2 − ξ is the stable orientation, while α = π/2 + ξ is an unstable orientation (see Fig. 12(b)) The results discussed above illustrate the dynamics when the initial particle orientation is co-planar with F and ∇η -i.e., φ 0 = 0 or φ = π.Fig. 14 plots the orientation angles φ and α over time when the starting angle is no longer co-planar with F and ∇η -i.e., For the aspect ratio parameter shown in this figure (A R = 13 for prolate and A R = 1/13 for oblate), there is no stable orientation and the spheroids continue to tumble.
φ 0 = 0 or π.We see that at long times, the angle φ → 0 or π -i.e., the orientation ends up in the same plane as F and ∇η.The angle α 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.

Not all spheroids have a steady orientation
Fig. 15 plots the steady orientation angles for prolate and oblate spheroids for different aspect ratio parameters.The steady orientations occur when dα dt = 0 and dφ dt = 0 in Eq. (5.3), which corresponds to the criterion: In the above equation, (λ 1 , λ 3 , λ 4 ) are the mobility coefficients for force-rotation coupling that were calculated in Sec. 4, which are only functions of the aspect ratio parameter A R .Fig. 15 show that for a wide range of A R , the above criterion is satisfied and a steady angle α se exists.The stable orientation α se for a prolate spheroid is closer to 0 and π compared an oblate spheroid, while the oblate spheroid has a stable angle closer to π/2.We also observe that for certain values of the aspect ratio parameter A R , The dashed curves correspond to when no viscosity gradient is present (β = 0), while the solid curve is when a viscosity gradient is present (β = 0.1).Different color curves correspond to different initial starting angles α 0 .Plot (a) illustrates a case when the spheroid attains a steady orientation, while (b) illustrates a case when the spheroid tumbles.
no steady orientation is reached.These situations occur for prolate spheroids between A R = 9 to 21, and oblate spheroids between 1/A R = 11 to 1/A R = 91 (see Fig. 15).At these aspect ratio parameters, the spheroid keeps tumbling and does not reach a steady state.This trend is illustrated vividly in Fig. 16 for both prolate (A R = 13) and oblate (A R = 1/13) spheroids.

Translation dynamics
Fig. 17 shows the spheroid's translation trajectories for the case when the force and viscosity gradient are perpendicular to each other.Two different dynamics occur depending on whether the spheroid obtains a stable orientation or not.In Fig. 17(a) 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 2007).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 (see Fig. 17(a)).
Conversely, in Fig. 17 (b) when the particle is at an aspect ratio that does not have a steady orientation, the particle will tumble throughout its sedimentation.In this case, the particle's motion will sediment in the gravity direction, but its trajectory will oscillate in the viscosity gradient direction (1-direction), with the oscillation period scaling with the tumbling time.

Governing equations
We now consider the most general case where F and ∇η are neither parallel or orthogonal to each other, but are inclined at an angle θ to each other.The external force points in the positive z-direction F = ẑ, while the viscosity gradient is as follows: Similar to before, the orientation vector is p = [sin α cos φ, sin α cos φ, cos α], where α and φ are the polar and azimuthal angles.To determine how these angles evolve over time, we note that the dynamics are a linear superposition of the cases described previously.
In other words, where dα dt and dα dt ⊥ are the variations in the polar angle from viscosity gradients parallel and perpendicular to the external force, given by Eq. (5.1) (using the positive sign) and Eq.(5.3a), respectively.The corresponding terms dφ dt and dφ dt ⊥ are the same quantities for the azimuthal angle, which is zero for dφ dt and Eq.(5.3b) for dφ dt ⊥ .The equation for particle translation is the same as Eq.(5.4).

Steady orientation angles
If a steady orientation angle exists, it will be in the plane spanned by F and ∇η as discussed previously -i.e., φ = 0. We set φ = 0 and determine the conditions under which dα dt = 0 in Eq. (5.7a).The criterion for a steady orientation angle is: (5.8) For illustration, Fig. 18 plots the steady orientation angles α se for different values of the angle θ between the external force F and ∇η.The results are plotted for prolate and oblate spheroids with aspect ratio parameter A R = 5 and A R = 1/5, respectively.We observe that α se varies between π/2 and π for prolate spheroids and between 0 and π/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 external force.For example, as θ → 0, we see α se → π/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.
Lastly, Fig. 19 provides a phase diagram that describes when a steady orientation exists for different particle shapes and viscosity gradient directions.When the viscosity gradient is parallel (θ = 0) or anti-parallel (θ = π) to the force, there always exists a stable, steady orientation, whereas when the viscosity gradient is perpendicular to the Here θ is the angle between the viscosity gradient ∇η and external force F , while A R is the aspect ratio parameter.
force (θ = π/2), there is a range of aspect ratio parameters A R where steady behavior does not exist.At other angles, we observe intermediate behavior between the two limits as illustrated in the figure.

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.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 ∇T far away from the particle.The fluid's viscosity depends linearly on temperature -i.e., η − η 0 = dη dT (T − T 0 ), and thus the viscosity field also varies spatially around the particle.If the thermal Peclet number is small and the temperature profile is steady, the temperature field will satisfy Laplace's equation inside and outside the particle: (see (Dassios 2012) for details): This equation is subject to the following boundary conditions: (a) T out → T 0 +∇T •x far away from the particle (|x| → ∞), and (b) on the particle surface, the temperatures and fluxes are continuous -i.e., T in = T out and n • ∇T out = k r n • ∇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 η(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 τ ex ij = (η(x) − η 0 )γ (0) ij .The mobility tensors described in Sec. 4 will take the same form, except the numerical values for the force/rotation mobility coefficients (λ 1 , λ 3 , λ 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 will recover the results described earlier in the manuscript.Otherwise, the viscosity field will have a more complicated spatial dependence, but the qualitative trends will likely remain the same for the steady orientations and the shape of the particle trajectories.We currently do not have any quantitative results for such an analysis (perhaps to be taken up later in a different paper).But as shown in (Shaik & Elfring 2021) for swimming spheres, we expect that the disturbance temperature field will only bolster the effects of spatial variation in viscosity, and may not lead to any novel effects.In Appendix E, we outline how one solves Laplace's equation around an ellipsoidal particle.

Conclusion
In this paper, we study a spheroid sedimenting in 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 Eq. (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 the external force, the spheroids reach a stable orientation where the longest axis is along to the viscosity gradient.We also show that for most initial orientation angles, the spheroid aquires 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.Depending on the particle aspect ratio, the spheroid may continue to rotate in this plane or reach a steady orientation.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 disk, 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 low P e limit.However, based on previous literature (Shaik & Elfring 2021) we believe that such an analysis may not yield any novel results not yet accounted for.
Finally, we remark that even though we employed a perturbative approach to the solution in the limit of weak viscosity gradients, we expect that the steady state behavior -namely the stable orientation of the spheroids -will remain unchanged even when the viscosity gradients become stronger.Stronger viscosity gradients will change the rate at which the stable orientation is attained, but not the value of the steady orientation per se.Lastly, a spheroid is a typical axisymmetric particle with no isotropy but fore-aft symmetry.We believe the qualitative results here will hold for other orientable particles with fore-aft symmetry, and thus can be a model representation of several systems in nature and in industry.Additionally, the current problem may 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.

Appendix A -Disturbance velocity for an ellipsoid in Stokes flow
Consider a reference frame at an ellpsoid's center of mass with axes aligned along the particle's principle axes.From (Kim & Karilla 2005, pg. 55), the Stokes velocity field around the ellipsoid from external force and torque is the following: In these formulas, no summation is assumed for repeated indices unless explicitly stated.
To obtain formulas for v trans ki and v rot ki in the reciprocal theorem, we substitute into Eq.(7.1) the force and torque that comes from unit translation and rotation, respectively.
In the above expressions, the expression for G n is: with ∆(t) = (a 2 + t)(b 2 + t)(c 2 + t) and λ(x, y, z) being the positive root of

Appendix B -Resistance formulae for an ellipsoid in Stokes flow
Let us consider a reference frame with the origin at the center 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 F U and R T ω are diagonal, while the cross-coupling term R F ω = R T U = 0.The diagonal elements are: where V = 4π 3 a 1 a 2 a 3 is the particle volume, and (χ 0 , α 1 , α 2 α 3 ) are elliptic integrals defined below: 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 = a 2 = a 3 , the coefficients c 1 − c 4 for the mobility matrix in Eq. (4.3a) and Eq.(4.3b) are: These coefficients have analytical formulae (see pgs 64 and 68 in Kim and Karilla).Using the notation in this paper, we obtain for prolate and oblate spheroids: where e = 1 − b 2 a 2 is the spheroid's eccentricity and L = ln 1+e 1−e .To get the nondimensional form used in the manuscript, we multiply c 1 and c 2 by 6πη 0 R = 6πη 0 (ab 2 ) 1/3 , and multiply c 3 and c 4 by 6πη 0 R 3 = 6πη 0 ab 2 .

Appendix C -Reciprocal theorem and O(β) solution
To delineate the O(β) correction to the particle kinematics -i.e., obtain the solution for (U (1) i , ω (1) i ) -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 polymeric 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 polymeric force and torque using the reciprocal theorem (Leal 1980).
First, we note that the fluid stress field at O(β) has two parts: One is the Newtonian part given by γ(1) ij − p (1) δ ij .The other part is polymeric, denoted as τ ex ij and given by: We note the important observation that the polymeric stress at O(β) 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 aux  (7.14)where 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 Note when deriving the above expression, we made use of the fact that the force and torque acting on the particle at O(β) is zero (F (1) i = T (1) i = 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: where in the above equation, the resistance tensors satisfy the following symmetry relationships: R F U = (R F U ) T , R T ω = (R T ω ) T , and R F ω = (R T U ) T .We will also write the auxillary velocity field in the volume integral for Eq. ( 7.15) as a linear combination of the rigid body motions: (7.17) where v rot ik and v rot ik are the velocity fields in the i direction induced by unit translation or rotation in the k direction.Substituting Eqs.(7.16) and (7.17) into (7.15) and eliminating U aux i and ω aux i yields the final result (Eq.(3.6)) stated in the manuscript.7.4.Appendix D -Simplification of mobility tensor M ijk Here, we show that Eq. (4.6) is equivalent to Eq. (4.7).To that end, we re-write Eq. (4.6) as: +λ 2 p i ǫ jkq p q Term2 +λ 3 p j ǫ ikq p q Term3 +λ 4 p k ǫ ijq p q Term4 (7.18)Without any loss of generality, we assume a particular orientation of the projection vector p i namely p i = δ i1 (and so on).Therefore, Eq. (7.19) may be written as: ellipsoid.The functions E p 1 and F p 1 are Lame functions of the first and second kind, defined in publication (Sten 2006).In Eq. (7.27), the coefficients B 1p are the following: ∂T ∂x (7.28a) ∂T ∂y (7.28b) ∂T ∂z (7.28c) where k = √ a 2 -c 2 and h = √ a 2 -b 2 .The geometric factors L 1,2,3 1 (a) take the following form (Sten 2006):

Figure 1 :
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 behavior when viscosity stratification is present -i.e., case (d)

Figure 2 :
Figure 2: Schematic of a prolate and oblate spheroid falling under an external force acting in the 3-direction.The viscosity gradient is along the 3-direction (parallel or antiparallel).The particle's orientation vector p makes a polar angle α ∈ [0, π] with respect to the sedimentation direction.

Figure 3 :
Figure 3: Schematic of a prolate and oblate spheroid falling under an external force F acting in the 3-direction, while the viscosity varies spatially in the 1-direction.The particle's orientation p makes a polar angle α ∈ [0, π] with respect to the 3-direction, and makes an azimuthal angle φ ∈ [0, 2π) in the 1-2 plane.
are the resistance tensors for a spheroid, which are given in Appendix B. The external force and torque are given in Eq. (3.2) .At the next order of approximation O(β), the momentum and mass balance equations become Stokes flow with an extra fluid body force b i :

Figure 4 :
Figure4: 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 T = 0.The radius and fluid viscosity are a = 1 and η 0 = 1, respectively.The results of the theory are from(Datt & Elfring 2019), expanded in Sec.3.2.1.

Figure 6 :
Figure 6: Verification of theory by plotting of Eq. (4.11) for different values of α, for a prolate spheroid with external force F and viscosity gradient ∇η in the positive zdirection.This situation corresponds to Case A shown in Fig. 5a.

Figure 8 :
Figure 8: Orientation angle α vs. time for prolate and oblate spheroids when the external force F and viscosity gradient ∇η are parallel or anti-parallel to each other.The left figures (a,c) correspond to prolate spheroids with A R = 5, while those the right figures (b,d) correspond to oblate spheroids with A R = 1/5.The top row (a,b) is the case when the F and ∇η 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 reduced order theory (solving Eq. (5.1)).The dimensionless viscosity gradient is β = 0.1.

Figure 9 :
Figure 9: Steady configurations attained by (a) prolate and (b) oblate spheroids when the external force F and viscosity gradient ∇η 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.

Figure 10 :
Figure 10: 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.

Figure 11 :
Figure 11: Particle trajectories for (a) prolate and (b) oblate spheroids when the external force and viscosity gradient are in the same direction (F = ẑ, ∇η = β ẑ).The dashed curves correspond to when no viscosity gradient is present (β = 0), while the solid curve is when a viscosity gradient is present (β = 0.1).Different color curves correspond to different initial starting angles α 0 .The prolate spheroid has A R = 5 while the oblate spheroid has A R = 1/5.

Figure 12 :
Figure12: Orientation angle α vs. time for prolate (A R = 5) and oblate (A R = 1/5) spheroids when the external force and viscosity gradient are perpendicular (F = ẑ, ∇η = β x).The dimensionless viscosity gradient is β = 0.1, and the particle initially starts in the plane of F and ∇η (i.e., φ 0 = 0).Solid curves are from full numerical simulations based on the reciprocal theorem, while the dashed curves are from the reduced order theory (solving Eq. (5.3)).

Figure 13 :
Figure 13: Schematic explaining the absence of steady orientations at α = 0 and α = π/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.

Figure 14 :
Figure14: Orientation angles α(t) and φ(t) for prolate and oblate spheroids when the external force and viscosity gradient are perpendicular (F = ẑ, ∇η = β x).The dashed curves show the evolution of φ, while the solid curves show the evolution of α.For all cases, the initial orientation is given by the ordered pair (φ 0 , α 0 ) = (π/3, π/4) and the dimensionless viscosity gradient is β = 0.1.The results show that φ → 0 or π, and hence the particle becomes co-planar with F and ∇η.

Figure 15 :Figure 16 :
Figure15: Stable orientations α se for prolate and oblate spheroids of different aspect ratio parameters A R when the external force and viscosity gradient are perpendicular to each other (F = ẑ ∇η = β x, β = 0.1).Regions that do not have data points are regions where the particle tumbles and does not exhibit a stable orientation.

Figure 17 :
Figure 17: Particle trajectories for prolate spheroids with (a) A R = 5 and (b) A R = 11 when the external force and viscosity gradient are perpendicular (F = ẑ, ∇η = β x).The dashed curves correspond to when no viscosity gradient is present (β = 0), while the solid curve is when a viscosity gradient is present (β = 0.1).Different color curves correspond to different initial starting angles α 0 .Plot (a) illustrates a case when the spheroid attains a steady orientation, while (b) illustrates a case when the spheroid tumbles.

Figure 18 :
Figure 18: Stable orientation angles α se for prolate and oblate spheroids when the viscosity gradient ∇η and the external force F are inclined at an angle θ to each other.

Figure 19 :
Figure19: Phase diagram demarcating the region in (θ, A R ) space where a stable orientation is reached (blue circles) and where the spheroid tumbles without reaching any stable orientation (yellow triangles).Here θ is the angle between the viscosity gradient ∇η and external force F , while A R is the aspect ratio parameter.