Active particles crossing sharp viscosity gradients

Active particles (living or synthetic) often move through inhomogeneous environments, such as gradients in light, heat or nutrient concentration, that can lead to directed motion (or taxis). Recent research has explored inhomogeneity in the rheological properties of a suspending fluid, in particular viscosity, as a mechanical (rather than biological) mechanism for taxis. Theoretical and experimental studies have shown that gradients in viscosity can lead to reorientation due to asymmetric viscous forces. In particular, recent experiments with Chlamydomonas Reinhardtii algae swimming across sharp viscosity gradients have observed that the microorganisms are redirected and scattered due to the viscosity change. Here we develop a simple theoretical model to explain these experiments. We model the swimmers as spherical squirmers and focus on small, but sharp, viscosity changes. We derive a law, analogous to Snell’s law of refraction, that governs the orientation of active particles in the presence of a viscosity interface. Theoretical predictions show good agreement with experiments and provide a mechanistic understanding of the observed reorientation process.

diffusiophoretic Janus particles and ciliated organisms like Paramecium or Opalina. In this case the interaction of the spatially varying viscosity with the active slip boundary conditions on a squirmer generically resulted in negative viscotaxis 40,41 . A different swimmer, Taylor's swimming sheet speeds up while moving along or against the gradients 42 . These theoretical models all rely on weak, diffuse gradients in the fluid viscosity. However recent experiments have shown very interesting particle dynamics in sharp viscosity gradients both for synthetic 43 and biological active particles 25 . In particular, we are interested here in experiments which probed the motion of Chlamydomonas Reinhardtii swimming across a sharp jump (or interface) in viscosity between miscible fluids 25 . Among other results, it was found that the algae would be quickly reoriented by the interface in viscosity and if the organism approached the interface at a sufficiently shallow angle could be reflected by the interface if going from low to high viscosity (see Fig. 1). Here, we develop a simple fluid dynamical model to unravel the physics underlying these experiments. We model the swimmers as spherical squirmers and focus on small, but sharp, viscosity changes. We show that the reorientation process is always in the direction of lower viscosity and derive a law, analogous to Snell's law of refraction, that governs the orientation of active particles in the presence of a viscosity interface. In analogy to ray optics, the refraction of the trajectory is always towards the medium of lower resistance. As we will show below, our theory (for pullers) matches well with experimental observations of Chlamydomonas Reinhardtii algae swimming across sharp viscosity gradients 25 . Our results are also quite similar to recent theoretical work modeling gliders moving across a substrate features a jump in frictional properties. In particular, the functional form of the reorientation law we find is identical to that found for gliders 44 . However that work, and other studies where the propulsive force is similarly fixed 39 , shows reorientation towards higher viscosities as one might expect due to the modulation of drag alone.
We organize the paper as follows. In the following section, we provide the essential details of our model and the resulting particle dynamics. We then interpret the implications of our model and compare them with experimental observations in "Results" section. We then provide some concluding remarks and finally the technical details concerning the mathematical methods used are left to "Methods" section.

A model for active particles crossing sharp viscosity gradients
We consider an active particle immersed in an otherwise quiescent fluid, moving near and across a region where the fluid has a relative sharp change of viscosity. This change in viscosity can be due to a corresponding variation in fluid temperature, salinity, or a nutrient dissolved in the fluid. Regardless of the origin, one expects sharp viscosity gradients to vanish due to diffusion over long times, but during the short time scales over which the particle crosses the interface, relatively sharp gradients can be stable 25 . For instance, Chlamydomonas Reinhardtii algae, with a characteristic size of ≈ 10 µ m, traveling at a body length per second take O(10 s) to approach and cross the interface while the salinity gradients take O(10 3 s) to vanish 25 . We assume that the fluid viscosity is prescribed and steady, and not significantly disturbed by the presence and activity of the moving particle as in previous work [39][40][41] ; however, this is an uncontrolled approximation because we assume that the sharp gradient persists, nevertheless we will show that this reasonably captures the experimental observations 25 . We note that Figure 1. Trajectories of Chlamydomonas Reinhardtii going from low to high (top) and high to low (bottom) viscosities from recent experimental work 25 . The trajectories on the left are relatively steep (closely aligned to the interface-normal) while those on the right are relatively shallow and display scattering at the interface when going from low to high viscosities. Image is from the paper by Coppola  www.nature.com/scientificreports/ the viscosity gradients considered here are distinct from the long-lasting viscosity differences that can exist at the interface between immiscible fluids with non-negligible surface tension that can dramatically affect (and even prevent) the particle crossing [45][46][47][48] . For simplicity we assume changes in only one direction and choose a coordinate with the z−axis oriented in the direction of change, such that η = η(z) . The viscosity changes from one uniform viscosity η(z → −∞) = η 0 to another η(z → ∞) = η 1 , and define a relative change in viscosity ε = (η 1 − η 0 )/η 0 . We will first assume the viscosity jumps (discontinuously) from η 0 to η 1 at z = 0 , as shown in Fig. 2. In this case the viscosity field may be written where H(z) is the Heaviside function. This representation is of course an idealization (as finite diffusivity would instantly smooth any discontinuity); however, we will later show that relaxing this assumption to smooth changes in viscosity leaves our results unchanged. To make mathematical progress we focus on small changes in viscosity such that |ε| ≪ 1 , but we believe the main physical picture holds for any ε.
The fluid flow generated by the active particle satisfies the incompressible Stokes equations where u is the velocity field and stress in a Newtonian fluid, σ = −pI + ηγ , where p is the pressure, γ = ∇u + (∇u) T is the fluid strain-rate tensor. The active particle swims with a translational velocity U and an angular velocity due to its activity. Thus, the velocity of the fluid on the surface of the particle S p , can be decomposed as where u s is the boundary velocity of the activity alone, r = x − x c , and x c = x c , y c , z c denotes the particle (center) position and ẋ c = U . Far from the particle the fluid remains quiescent, hence Here we prescribe the activity of the particle u s , but the translational and rotational velocity are fixed by the dynamic constraints on the particle. The particle is inertialess and neutrally buoyant and with no external forcing acting on it therefore the hydrodynamic force and torque on the particle must vanish where n p is a unit normal to the particle surface.
as |r| → ∞. www.nature.com/scientificreports/ Specifically, we model the active particle as a spherical squirmer of radius a. In the squirmer model, the details of the surface activity of the swimmer are coarse-grained into a prescribed tangential slip velocity on the surface of a spherical particle [49][50][51] . This model is particularly well-suited for ciliated microorganisms like Paramecium and Opalina that propel by synchronously beating numerous very small cilia on their surface, or diffusiophoretic Janus particles, which propel due to the motion of a thin layer of fluid on their surface as a result of chemical gradients 15 . The slip velocity is generally decomposed into Legendre polynomials (called squirming modes) in the form where p is the particle orientation, P n is the Legendre polynomial of degree n and B n represents the coefficients of the squirming modes. In homogeneous Newtonian fluids, the B 1 mode alone determines the swimming velocity (we assume here B 1 > 0 ), whereas B 2 mode is the slowest decaying contribution to the far-field flow, furthermore the second mode determines whether propulsion is primarily from the front or the back of the swimmer. Organisms, such as E. coli, which produce propulsion from their rear end are called pushers, have B 2 < 0 whereas those that pull the fluid in front of them using their flagella are called pullers, such as Chlamydomonas Reinhardtii, and have B 2 > 0 . Swimmers, with propulsion that is not distinctly fore or aft, such as Volvox carteri with flagella uniformly distributed on its surface, are called neutral and are well described by setting B 2 = 0 . Here we look at only the effects of the first two modes, and while one is generally only well justified in neglecting higher-order modes in the far-field, we find dynamics here to be well captured by just the B 1 mode.
The fluid flow field and particle velocity can be determined simultaneously by solving the Stokes equations for a force and torque-free particle. Here, we take a perturbative approach; when ε → 0 the viscosity is uniform and the solution to a single squirmer is well known, we then obtain the leading order correction in terms of the viscosity jump ε , by means of the reciprocal theorem (see "Reciprocal theorem" subsection in "Methods" for more technical details).
At the leading order in ε , the particle moves through a homogeneous Newtonian fluid of viscosity η 0 and its velocity is well known 49,50 The viscosity variations relative to η 0 are captured at the next order and the particle swims due to these variations at velocity where we have assumed a regular expansion in ε , U(ε) = U 0 + εU 1 + O(ε 2 ) , and (ε) = 0 + ε 1 + O(ε 2 ) . Here n = e z is the interface normal pointing from fluid of viscosity η 0 to that of viscosity η 1 while the functions A(z c ) , B(z c ) , C(z c ) , D(z c ) , f (z c ) , and g(z c ) depend on the particle's separation from the interface and are given in "Methods" section. The piecewise behavior of these functions and of the velocities U 1 , 1 is due to the fact that the particle is in contact with the viscosity interface when |z c | ≤ a and otherwise not. It is important to note that for |z c | > a the particle is still affected by the presence of the viscosity change due to hydrodynamic interactions mediated by the fluid at a distance.
We show later that the reorientation caused by 1 is always towards the low viscosity fluid. The physical reason for this is similar to that observed in weak gradients 40 . Both the rigid-body drag felt by the particle and the propulsive thrust generated by surface activity are altered by the presence of viscosity changes in the fluid (both due to contact and hydrodynamic interactions). Changes in the thrust tend to cause rotation towards regions of lower viscosity, while the opposite is true for changes to the drag, and in all instances the former dominates the dynamics (for spherical squirmers). The speed changes represented by U 1 here are somewhat more complex than those observed in weak gradients, with even the neutral swimmer speeding up or slowing down depending on its position and orientation relative to the interface. However, these speed changes do not affect leading order trajectory of the particle. In order to quantify particle reorientation, we simply integrate particle velocities. Noting that ẋ c = U and ṗ = × p , we substitute the leading order results for the particle velocities and project onto the interface normal direction n = e z to obtain where the angle between the particle direction and the interface is defined by p · n = cos θ . Combining these equations gives www.nature.com/scientificreports/ where β = B 2 /B 1 . This differential equation entirely captures the leading order effect on particle orientation θ of a viscosity jump at a distance z c from the particle center. We note that any effect of viscosity on the translational velocity of the particle would enter at O ε 2 in (13) and so is negligible compared to the leading order terms for ε ≪ 1 . As we will show, Eq. (13) is straightforward to integrate analytically for neutral squirmers, β = 0 , and easily integrated numerically for pushers ( β < 0 ) and pullers ( β > 0 ), and the rest of this paper are results and discussion that arise out of it.

Results
We begin first with analytical results for neutral squirmers, β = 0 , before proceeding to present numerical results for pushers ( β < 0 ) and pullers ( β > 0 ) and a comparison to recent experiments for pullers.
Neutral squirmers. The reorientation of a neutral swimmer can be immediately understood by examining the instantaneous rotational dynamics. To leading order equation (12) with B 2 = 0 simplifies to Noting that B 1 and f (z c ) are both positive, when ε > 0 we see that θ = 0 is an unstable fixed point and all orientations flow towards θ = ± π . Conversely, for ε < 0 all orientations flow to θ = 0 . This means that no matter the orientation or position, the particle is always reorienting to align along n = e z and point in the direction of the lower viscosity, consistent with results for squirmers in weak viscosity gradients 40,41 . Because f ∝ z −4 c , the reorientation rate decreases very quickly with distance from the interface, and thus the reorientation process is ultimately dominated by contact with the interface. One consequence of these dynamics is that a particle going from low to high viscosity can be scattered off the interface depending on its incident orientation.
To quantify the reorientation we note that (13) is separable when β = 0 , integrating we obtain where θ i is the orientation at an initial position z i and likewise θ f is the final orientation at z f . We define the 'total' reorientation caused by the interface as the particle crosses from far on one side to far on the other to be the limit when z i → −∞ and z f → ∞ (when the particle goes from η 0 to η 1 ). In this case the integral simply equals 1/3 and the total reorientation is given by the formula This formula bears a striking similarity to Snell's law of refraction, except here the 'relative refractive index' is given by the exponentiated relative viscosity difference. The reorientation is independent of the speed of the particle due to the linearity of the Stokes equations, in this case both the thrust generated by the particle and the drag felt by the particle would be proportional to the B 1 mode. This form of reorientation law, sin θ f = e α sin θ i , was found for gliders moving across a substrate featuring a jump in frictional properties 44 . In that case α = −2aζ rt /ζ rr where ζ rr , ζ rt are torque-rotation and torque-translation resistance coefficients of the particle respectively. The similarities arise because in both cases the particles are subject to linear drag laws.
Unlike the refraction of light, or gliders on a substrate, squirmers interact (hydrodynamically) with the interface from any point in space, but the functional form of the interaction, given by f (z c ) changes upon contact. Because of this we integrate equation (15) in multiple stages, separately accounting for the particle's approach to the interface (z c = −∞ → −a , θ = θ i → θ −a ) , crossing the interface (z c = −a → +a, θ = θ −a → θ a ) and the departure from the interface (z c = +a → +∞ , θ = θ a → θ f ) as illustrated in Fig. 3a. In this way we can quantify the starting (θ start ) and ending (θ end ) orientation in each of these stages where for θ start = {θ i , θ −a , θ a } , and θ end = θ −a , θ a , θ f , we find α = ε 32 , 7ε 16 , ε 32 , in that order in each of the three stages. We can also relate the particle's initial and final orientations by combining Eq. (17) in all three stages. We notice that the amount of reorientation during the approach and departure from the interface is the same. However, the large value of α means that the reorientation process is dominated during contact with the interface (θ −a → θ a ) . And, as discussed earlier, the reorientation process is always in the direction of lower viscosity. In analogy to ray optics, the refraction of the trajectory is always towards the medium of lower resistance. A similar preference was also shown by active particles in linear or diffuse viscosity gradients 24,40,41 , and as we will show below (for pullers), matches well with experimental observations of Chlamydomonas Reinhardtii algae swimming across sharp viscosity gradients 25 . In contrast, studies done where the propulsive force is fixed, both for swimmers in diffuse viscosity gradients and gliders across a frictional substrate, show reorientation towards higher viscosities as one might expect due purely to the modulation of drag.  (16) we assume that the trajectory from one side to the other is physically realizable. This is not always the case. If the active particle is swimming towards a higher viscosity η 1 > η 0 , with a sufficiently shallow angle it may be reoriented back (reflected by the interface). But note, due to hydrodynamic interactions, the particle may be reoriented back before even coming into contact with the interface, or even after completely crossing the interface. To examine these phenomena we find the limit of validity of (16), which we define θ i = θ crit which occurs when θ f = π/2 , that is the swimmer is tangent to the interface at z f → ∞ , this case we find Therefore, for η 1 > η 0 Eq. (16) is valid only when θ i ≤ θ crit , as a particle with an initial angle, θ i > θ crit (a sufficiently shallow angle of approach to the interface), will be reflected back (it will not reach z f → ∞ ). We can likewise define critical initial angles such that the particle does not cross the interface θ crit,a = arcsin exp [15(η 0 − η 1 )/(32η 0 )] , or even touch the interface θ crit,−a = arcsin exp [(η 0 − η 1 )/(32η 0 )] , using (17). We see that, much like the total reorientation, θ crit in (18) is dominated by particles which are scattered at the interface, with only a narrow set of initial angles that lead to reflection before or after contact with the interface. Regardless of where the particle is reflected, the entire scattering process is symmetric (about the point when θ = π/2 as shown in Fig. 3b) and hence obeys the reflection law for all particles when θ i > θ crit , as previously shown for gliders 44 .
We have thus far assumed a physically unrealistic discontinuous viscosity change and a simple neutral squirmer in order to derive these simple formulas. In the following sections, we relax these assumptions and find that neither significantly impact the reorientation process.

Smooth gradients.
To investigate the reorientation and scattering of the particle due to a viscosity change that varies smoothly (due to the effects of diffusion), instead of a Heaviside function in (1) we say H(z) = (1 + tanh(kz))/2 where k > 0 and 1/k is the effective length scale over which the viscosity varies between η 0 and η 1 . Hence, this viscosity variation approaches the discontinuous profile used previously as k → ∞ . Calculation proceeds similarly to that in the sharp gradients. The angular velocity required for this calculation is found by substituting the tanh viscosity profile in the reciprocal theorem (27) and its expression looks the same as that found in sharp viscosity gradients (10) except for the functions f (z c ) and g(z c ) which are given in "Methods" section.
Despite the differences in viscosity profile and angular velocity, we find no difference in the overall reorientation (θ i → θ f ) between smooth and sharp viscosity gradients. The law governing the reorientation in smooth gradients is identical to that found in sharp gradients (16). This implies that the critical orientation required for scattering in smooth and sharp viscosity gradients is also the same. It appears that, as with the refraction of light, the interface between fluids of differing viscosity can be smoothed out and the total reorientation remains (18) Figure 3. Schematic showing the reorientation of an active particle as it ( a ) crosses a viscosity interface or ( b ) gets reflected by the interface. This reorientation depends largely on the viscosity difference η 1 − η 0 and the particle is reflected only when going from low to high viscosity if its initial orientation is sufficiently shallow, θ i > θ crit . www.nature.com/scientificreports/ unchanged. This is surprising, because unlike light, the active particle interacts non-locally with the entire medium at once at all times due to hydrodynamics. (13) is not separable and so we compute the reorientation numerically for β = 0 . We find that the reorientation and scattering of pushers or pullers are similar to those of neutral swimmers with only a weak dependence on the squirming ratio β . See Fig. 4a for the reorientation of the swimmers crossing the interface and Fig. 4b for the critical orientation required for scattering, obtained from the numerical solution of (13). Hence, the pushers and pullers, like the neutral swimmers, orient towards regions of lower viscosity. Going from high to low viscosity ( ε < 0 ) pullers rotate slightly less while the pushers rotate slightly more than the neutral swimmer. Conversely going from low viscosity to high viscosity ( ε > 0 ) pullers rotate more while the pushers rotate less than a neutral swimmer thereby very weakly changing θ crit as shown in Fig. 4b. Part of the reason for the weak dependence of reorientation on the B 2 mode occurs because the rotation caused by this mode before and after crossing the interface is in the opposite direction (as g(z c ) is an odd function) and hence counteract each other (they do not cancel due to the cosine term in (13)). Conversely, as discussed previously the rotation caused by the B 1 mode is always in the same direction before and after crossing the interface (as f (z c ) is an even function).

Pushers and pullers. For pushers and pullers the differential equation governing the reorientation
The weak dependence of reorientation on β can be leveraged to find the reorientation experienced by the pushers and pullers analytically. This is achieved by expanding the leading order (in ε ) orientation θ in terms of β and solving Eq. (13) at each order in β as shown in "Methods" section. In principle, such a perturbation holds for only |ε| ≪ |β| ≪ 1 but the weak functional dependence yields accurate results up to |β| ≈ 10 for any |ε| ≪ 1.

Comparison to experiment.
We now compare our theory with recent experiments conducted with both wild-type (wt) and short-flagellated (sfl) Chlamydomonas Reinhardtii (CR), swimming across sharp viscosity gradients 25 . Just as we have predicted above, the CR were found generically to reorient towards lower viscosities and there was a critical angle, going from low to high viscosity, past which the swimmers would be reflected by the viscosity interface. The initial and final orientations of the algae were recorded 1s before reaching and 1s after crossing the interface in the experiments. Using experimental velocities this equates to z i ≈ −3a and z f ≈ 3a where a = 5µ m is the approximate swimmer radius, and we use these values in our theory for comparison. In the experiments, the viscosity of one fluid (water) was held constant ( η = 10 −3 Pa s) and a variety of different viscosities were used for the other fluid from 2 to 62× greater by dissolving varying concentrations of methylcellulose in water, resulting in relative viscosity differences |ε| = 0.5-61. We compare our asymptotic theory, which assumes ε ≪ 1 , only to the smallest values ε = −0.5, 1 representing particle motion from high to low and low to high viscosities respectively. We note that the experimental data indicate small but systematic reorientation even in homogeneous Newtonian fluid 'control' experiments. In order to remove this effect we subtracted the reorientation reported in homogeneous fluids from that in finite viscosity gradients and compared the difference with the theory. The experiments reported different swim speeds for wt and sfl algae, but found similar reorientation in viscosity gradients. This is consistent with our theory; our model predicts that the reorientation of a neutral swimmer is independent of speed (see Eqs. (15)-(18)), while the effect of the squirming ratio β is very weak for pushers and pullers, as shown in Fig. 4. In light of this, we represent both algae by a single squirmer whose B 1 mode is known from the known swim velocity of wt algae in homogeneous fluid 2B 1 /3 ≈ 100 µm/s and B 2 value from the stresslet exerted by same algae 10 pN×10 µ m ≈ 4πηa 2 B 2 found in other experiments 52 , ultimately yielding the squirming ratio β = 2 . Lastly, the algae in experiments displayed diffusive dynamics at large time scales, with translational and rotary diffusivities D T , D R respectively. But here we neglect both diffusivities as the former does not affect the reorientation while the latter is small compared to angular velocity for the viscosity jumps of interest ε = O(1) , where D R ≈ 1 s −1 < | | which is at least 4 s −1 .  www.nature.com/scientificreports/ Our theoretical model matches experimental observations reasonably well, in the aforementioned parameter regime. In Fig. 5a, we show the reorientation ( θ f vs θ i ) for swimming from high to low viscosity ( ε = −0.5 ) while Fig. 5b shows swimming from low to high viscosity ( ε = 1 ). Figure 5c shows particles reflected by the viscosity interface. In all cases our model somewhat over-predicts the amount of reorientation but captures nicely the general qualitative features observed in experiments. The largest discrepancy between theory and experiments occurs at shallow angles of approach to the interface (θ i > π/5) , when particles undergo large changes in orientation. In particular, experimental data does not well satisfy a symmetric reflection law for particles that are reflected at the interface as shown in Fig. 5c, but the data in this sensitive regime is somewhat limited and further data may clarify this discrepancy. Over-predicting the reorientation naturally leads to a critical angle (θ crit ≈ π/5) found in theory that is lower than that reported in the experiments (θ crit = π/3) , where a shallower approach to the interface is needed to scatter (see Fig. 5d).
Quantitative differences between experiments and theory are not surprising as our perturbative approach assumes ε ≪ 1 while in experiments at best we have ε = O(1) . Another possible cause of quantitative discrepancy may be due to confinement of the algae, in a microfluidic channel of height 20 µ m, in the experiments unlike the free-space assumption made in the theory. We approximated Chlamydomonas by a spherical squirmer whose gait remains fixed in viscosity gradients, while in reality alga has a spherical body with two flagella in front, whose beating pattern (or gait) varies with viscosity 25 ; we also assumed that the swimmer does not stir the viscosity field due to its motion and any mixing of the fluid in experiments is likely to weaken the effect of viscosity differences on reorientation and we expect these issues in particular to be exacerbated with shallower approachs to interface (θ i > θ crit ) , when particles spend a relatively large amount of time close to the interface. Finally, we neglected any density gradients that inevitably exist in viscosity gradients, and which preferentially alter motion of particle going from high to low viscosity (unlike in the opposite direction) due to relatively large volume of fluid dragged by the particle 43 .   25 , we developed a simple analytical model for active particles swimming across sharp changes in the viscosity of the suspending fluid. We found that pushers, pullers and neutral swimmers all interact similarly with the interface. Swimmers are generically reoriented towards the region of lower viscosity (as found in previous studies with weak gradients 40,41 ). As a result, if active particles approach a viscosity interface at a sufficiently shallow angle they can be reflected if swimming from low to high viscosity; otherwise, they simply cross the interface undergoing a degree of reorientation set by the relative viscosity difference. This is similar to the refraction or reflection of the light due to a change in refractive index and the law we derive governing the reorientation of neutral swimmers similar to Snell's law of ray optics (as previously shown for gliders on a frictional substrate 44 ). Our theory compares very well with experimental observations 25 and provides a simple model for the dynamics of active particles in fluids with inhomogeneous viscosity. These results suggest that tailoring the mechanical properties of fluids can be an effective method to control particle speed [53][54][55][56] and orientation 57-60 , ultimately organize active matter systems.

Methods
Reciprocal theorem. The dynamics of a force-free and torque-free active particle in a fluid medium of arbitrary rheology is given by 61 where U = [U ] T is a six-dimensional vector containing the swimmer's translational and angular velocities, and likewise F = [F L] T contains force and torque. R FU is the resistance tensor for the particle in a fluid of uniform viscosity η 0 , F s is the thrust force and torque due to particle activity in a homogeneous Newtonian fluid of viscosity η 0 , while the additional force F NN accounts for the changes in the rheological properties (viscosity) of the fluid. The formulas are obtained by using the reciprocal theorem, by projecting onto a known auxiliary flow as an adjoint solution (denoted by a hat), where V denotes the entire fluid volume outside the particle and τ NN = σ + pI − η 0γ = (η − η 0 )γ represents the extra (deviatoric) stress due to changes in viscosity from η 0 . T U , Ê U are linear operators relating the stress and strain-rate in the fluid to particle velocity, σ =T U ·Û and γ = 2Ê U ·Û in a fluid of homogeneous viscosity η 0 .
To facilitate the evaluation of swimming velocity in Eq. (20), we assume small relative viscosity differences ε ≪ 1 and regular perturbation expansion for any functional dependence on ε , for example h(ε) = h 0 + εh 1 + . . . . In this way, because the extra stress due to viscosity changes is O(ε), the swimmer is moving through a homogeneous Newtonian fluid of viscosity η 0 to leading order and its velocity is well known The effects of viscosity variations relative to η 0 are captured at the next order, where the swimming velocity is Here, τ NN,1 = η 0 H(z)γ 0 and γ 0 = ∇u 0 + (∇u 0 ) T is the rate of strain tensor associated with the leading order flow u 0 . These small viscosity variations εη 0 H(z) alter the velocity of swimmer in homogeneous fluid U 0 , 0 by a small correction U 1 , 1 . An evaluation of integrals in Eqs. (26), (27) with discontinuous viscosity jump (H is a Heaviside function) yields this correction as When we assume a smooth viscosity profile H(z) = (1 + tanh(kz))/2 , we obtain Pushers and pullers, the effect of β. In order to find the leading order effect of the second squirming mode for pushers and pullers, we assume β ≪ 1 and perform a regular perturbation expansion of the orientation in β Here, we assumed |β| ≫ |ε| and retained the terms at O(β) unlike those at O(ε) . We substitute this expansion in (13) and solve the resulting equation at each order in β . At zeroth order, pushers or pullers become neutral swimmers and the orientation is given by (16). Any deviations relative to the reorientation of the neutral swimmer are captured at the next order where The initial condition is θ 1 = 0 as z c → −∞ . We solve (37) in multiple stages, separately accounting for the reorientation during the interface approach, crossing and departure. We find as the particle touches the interface. Then (28) U 1 = B 1 A(z c )n + B(z c )p + B 2 C(z c )n + D(z c )(pp · n) + E(z c )(n · p) 2 n , (29) 1 = B 1 f (z c ) + B 2 g(z c )(n · p) (n × p), (30) A(z c ) = a 3 (−a 2 + 3z 2 c ) 24z 5 c , B(z c ) = a 3 (−a 2 + z 2 c ) 24z 5 www.nature.com/scientificreports/ as the particle crosses the interface and eventually to as the particle departs away from the interface. Accounting for the leading order reorientation, the final orientation of pushers or pullers as they cross and go far ahead of the interface is where θ 0f = θ 0 | z c →∞ follows from (16) as sin θ 0f = e ε/2 sin θ i .

Data availability
The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.