Slippery rheotaxis: new regimes for guiding wall-bound microswimmers

Abstract

The near-surface locomotion of microswimmers under the action of background flows has been studied extensively, whereas the intervening effects of complex surface properties remain hitherto unknown. Intending to delineate the shear-driven dynamics near a planar slippery wall, we adopt the squirmer model of microswimmers and employ a three-dimensional analytical-numerical framework in bispherical coordinates. It is interpreted that both the self-propulsion and the external shear flow are redistributed due to hydrodynamic slippage, followed by modulations in the thrust torque on the microswimmer. Phase portraits of the quasi-steady dynamics indicate that the stable upstream swimming states, known as ‘rheotaxis’, are significantly modulated by the slip length compared with the no-slip case. For puller swimmers, an intricate interplay among the modulated interfacial friction near a slippery surface, velocity gradients of the shear flow and the strength of the squirmer parameter promotes a critical shear rate beyond which a wide range of new rheotactic states exist. Consequently, an escaping microswimmer may exhibit rheotaxis or an existing rheotactic state annihilates due to crashing. Although stable states are absent for pushers without steric interactions, transitions from escaping and undamped oscillations to ‘rheotaxis’ occur in the presence of wall repulsion, but only until the other characteristics are overwhelmed by escape due to enhanced shear. Disclosing the ability of hydrodynamic slippage in broadening the scope of migration against a background flow for a wide range of parameters, the present work paves the way for investigations on the entrapment of microswimmers near complex pathways or sorting using selective rheotaxis.

1. Introduction

The booming importance of microswimmers, mainly in the biomedical sector and environmental remediation, has fueled the research on both living and synthetic microswimmers surrounded by various challenging real-life environments. In addition to developing an enhanced understanding of common biophysical processes, e.g. swimming of spermatozoa through female oviduct (Guidobaldi et al. 2015; Ishimoto & Gaffney 2015) and propagation of bacterial infection (Harkes, Dankert & Feijen 1992), the potential applications of microswimmers range from lab-on-a-chip devices (Denissenko et al. 2012; Bechinger et al. 2016), assisted fertilization (Magdanz et al. 2017; Bunea & Taboryski 2020), targeted drug delivery using bio-hybrid microbots (Wang & Gao 2012; Park et al. 2017; Li et al. 2022) to water treatment using Chlamydomonus sp. (Escudero et al. 2020) and other artificial microswimmers (Poddar, Bandopadhyay & Chakraborty 2019; Chen et al. 2021; Yuan et al. 2022).

Microswimmers often have to swim through complex physiological environments by overcoming the effects of a counterflow, thus preventing their downstream transportation. Conversely, a background flow can assist the entrapment of bacteria near surfaces and facilitate surface attachment, which is an essential step for biofilm formation (Rusconi, Guasto & Stocker 2014). The background conditions have also been found to influence the migration of droplets (Poddar et al. 2018, 2019; Mantripragada & Poddar 2022). The capability of certain microswimmers to change their orientation in response to velocity gradients is known as rheotaxis (Bretherton & Rothschild 1961; Miki & Clapham 2013; Ishimoto & Gaffney 2015; Sharan et al. 2022). A common example of this motion characteristic exists during the fertilization process where sperms have to be navigated through a long distance to reach the female eggs (Roberts 1970; Miki & Clapham 2013; Kantsler et al. 2014). Positive rheotaxis of motile spermatozoa in a Poiseuille flow was reported by Bretherton & Rothschild (1961). Later, Rothschild (1963) explained the boundary accumulation of sperms due to the difference of drag forces on the head and tail, similar to a weather vane. Under a similar background flow, Kantsler et al. (2014) observed that the sperm cells migrate differently for varying flow conditions, e.g. upstream spiraling along the bounding substrates and downstream advection for low and high shear velocities, respectively. These motion attributes were described as a combined effect of shear flow, chirality of flagellar beat and steric interaction with the substrate.

Rheotaxis was also reported in bacteria Bacillus subtilis (Marcos et al. 2012), Escherichia coli (Kaya & Koser 2012), and artificial microswimmers (Palacci et al. 2015). Subsequently, the theoretical analysis of Uspal et al. (2015) revealed that rheotaxis could also occur for spherical active particles, which do not possess shape asymmetry similar to an elongated micro-organism, such as spermatozoa. They showed that the rheotaxis of spherical microswimmers stems from a mechanism of shear-induced rotation near a hard surface, leading to a constrained motion in the shear plane towards the upstream direction at a steady height and orientation. In contrast, theoretical investigations on a virtual monoflagellate Leishmania mexicana promastigote (Walker et al. 2018) predicted no general stable guided taxis under background shear in the absence of a steric contact force at the wall. However, Ishimoto (2017) showed that the inclusion of wall repulsion results in rheotactic states for disc (two dimensional) squirmers, while the same effect further stabilizes the motion of spherical (three dimensional) squirmers of different types.

Although rheotaxis is possible in bulk fluid (Marcos et al. 2012; Kumar & Ardekani 2019), various attributes of rheotaxis were found to be greatly influenced in the neighbourhood of a liquid–solid interface (Hill et al. 2007; Kaya & Koser 2012; Mathijssen et al. 2019). The microfluidic experiments of Hill et al. (2007) illustrated that the bacteria at the channel centreline have an enhanced tendency of a cross-stream drift towards the channel walls while orienting towards the upstream direction. An enhanced accumulation of micro-organisms at the walls in a uniform shear and a parabolic flow was also reported through computer simulations Chilukuri, Collins & Underhill (2014). Similarly, the work of Mathijssen et al. (2019) provided insights into different rheotaxis regimes near a surface, e.g. upstream shifting of swimming orientation, oscillatory rheotaxis and coexistence of rheotactic migration up or against the vorticity direction of the background flow, and also identified the transitional shear rates. Nevertheless, the near-surface accumulating tendency of microswimmers has also been observed even in a quiescent environment (Li & Tang 2009; Tailleur & Cates 2009; Li et al. 2011; Elgeti & Gompper 2013; Kantsler et al. 2013). In addition, different other intriguing near-surface phenomena, such as bacterial entrapment for biofilm formation (Costerton et al. 1987), directional circular motion of cells (Lauga et al. 2006; Di Leonardo et al. 2011), pairwise dancing of Volvox (Drescher et al. 2011) and tumbling trajectories of bacteria E. coli (Kantsler et al. 2013), etc. were reported in the absence of a background flow. On the other hand, a significant volume of theoretical studies (Berke et al. 2008; Shum, Gaffney & Smith 2010; Crowdy 2011, 2013; Ishimoto & Gaffney 2013; Li & Ardekani 2014; Poddar, Bandopadhyay & Chakraborty 2020) has shed light on the crucial role of hydrodynamic interaction with the confining substrates in affecting the biophysical dynamics of microswimmers.

The interfacial properties at the solid–fluid or fluid–fluid interfaces were reported to have a substantial contribution in regulating the motion characteristics of a microswimmer near a confining boundary (Lemelle et al. 2013; Lopez & Lauga 2014; Hu et al. 2015; Pimponi et al. 2016). The deviation from the no-slip condition at a boundary leads to modulations in the interfacial friction at the micro- and nano-scale (Chakraborty 2008; Das et al. 2015; Maduar et al. 2015). The said deviation is often quantified by a slip length, which is the fictitious distance below the physical wall where the no-slip condition would apply. While, for the hydrophilic surfaces, the slip length is negligible, the same parameter can be significantly higher for the hydrophobic surfaces (Huang et al. 2008; Bocquet & Charlaix 2010). Different physical sources behind high slip length has been proposed in the literature. For instance, the slip length can go upto micrometres when the surface is in contact with a bacterial polymeric solution or due to a coating of monolayers of hydrophobic molecules (Tretheway & Meinhart 2002, 2004; Lauga, Brenner & Stone 2007). Also, air bubble entrapment between the asperities of micro- and nano-structured surfaces can be treated as an effective continuous partial slip boundary condition (Choi & Kim 2006; Joseph et al. 2006; Lee et al. 2008; Nizkaya et al. 2015).

In relation to microswimmers, Hu et al. (2015) predicted a transition of model E. Coli bacteria from circular to snaking trajectories due to alterations in the slip length. Similarly, the experimental observation of reversed circular motion of E. Coli due to added polymeric inclusions were estimated to be an effect of intensified slip (Lemelle et al. 2013). Lopez & Lauga (2014) employed a far-field analysis based on a force dipole swimmer incorporating non-zero slip lengths and reported that slip induces an additional rotation towards the wall, leading to attraction of pusher-type microswimmers. Very recently, Poddar et al. (2020) provided theoretical insights into the effects of high slip on the near-wall trajectories of different types of spherical micro-organisms. They also unveiled that slip length may be chosen as an effective control mechanism for switching from a scattering trajectory to wall entrapment. In another work, Ketzetzi et al. (2020) experimentally found that artificial microswimmers show augmented swimming speeds near a hydrophobic surface. These studies considered a quiescent flow condition only, and the effects of different background flows were never looked into. At the other extreme, the literature on rheotactic migration near confinement is limited to the no-slip condition at the wall. In view of the effects brought in by the hydrodynamic slip, it is anticipated that the coupling of the slip condition with the velocity gradients of an external flow (Loussaief, Pasol & Feuillebois 2015) may non-trivially alter the conditions of stable swimming states, thereby triggering unexplored attributes of rheotaxis.

In the present work we attempt to address the above shortcomings in the literature by formulating a mathematical model of a shear-driven spherical microswimmer in the vicinity of a slippery plane wall. The Navier slip condition has been employed to track the surface wettability condition at the wall. By incorporating the effects of an arbitrary slip length in a squirmer model, the present study stands apart from the widely employed models with asymptotically small slip lengths (Swan & Khair 2008; Willmott 2008). In addition, using the bispherical coordinate system for obtaining an exact solution of the creeping flow problem, the present model aptly captures the hydrodynamics at any distance from the wall, outside the scope of an image-singularity-based analysis (Lopez & Lauga 2014). Furthermore, the effects of the Navier slip condition both on the shear flow and the self-propulsion make the outcome of the study unpredictable beyond a simple linear superposition of the slip effects and the corresponding effects with a no-slip boundary condition (Uspal et al. 2015). We have further performed a detailed analysis of the quasi-steady dynamics of both puller- and pusher-type microswimmers and investigated the effects of different important dimensionless parameters, e.g. dimensionless slip length, shear rate and the squirmer parameter. We also discuss the effect of steric contact interactions on rheotaxis by employing a repulsive force at the wall. It has been found that beyond a critical shear rate, enhancement of slip length can either create new rheotactic states or destroy them depending on a critical interplay between the slip-induced effects and shear flow.

2. Mathematical formulation

In the present problem, we consider a microswimmer immersed in a background pure shear flow $\tilde {\boldsymbol {u}}^{(\boldsymbol {ex})}_{\infty }=\dot {\gamma } (\tilde {l}_s+\tilde {z}) \boldsymbol {e}_x$ near a planar wall that obeys the Navier slip condition (Navier 1823), as depicted in figure 1. The microswimmer is considered to have a spherical body of radius $a$, and its centre ($O$) is located at a vertical distance $\tilde {h}$ from the neighbouring slippery wall. The slip length ($\tilde {l}_s$) can be interpreted as the distance below the plain wall where the extrapolated velocity vanishes. The dimension of the microswimmer is assumed in the range ${O}(10^1\unicode{x2013}10^2) \ \mathrm {\mu } \mathrm {m}$. The orientation of the microswimmer is represented by the director $\hat {\boldsymbol {p}}$, defined as $\hat {\boldsymbol {p}} = p_x \boldsymbol {e}_x + p_y \boldsymbol {e}_y + p_z \boldsymbol {e}_z$, where $p_x = \cos (\theta _p)\cos (\phi _p)$, $p_y = \cos (\theta _p)\sin (\phi _p)$ and $p_z = -\sin (\theta _p)$. Here, we make an assumption that the slip length is uniform along the wall and that the microswimming properties are unaffected by surface texture. Motility of puller- and pusher-type microswimmers are illustrated in the inset, where the blue arrows indicate the mobility direction of the surrounding fluid and green arrows orient towards the local forcing direction of the model swimmer to the fluid if observed from the fixed frame. At this point, we would like to highlight the significance of choosing a three-dimensional (3-D) squirmer model instead of a two-dimensional (2-D) model (Ishimoto & Crowdy 2017). The physical phenomenon of upstream swimming or rheotaxis, which is the central theme of the present work, cannot be captured by a 2-D model, while a 3-D model was shown to be successful in this context (Uspal et al. 2015). It was also shown by Ishimoto & Crowdy (2017) that a 2-D squirmer can only show stable swimming if the Hamiltonian symmetry is broken by perturbations such as weak fluid viscoelasticity (Yazdi, Ardekani & Borhan 2015) or repulsive potential at the wall. Thus, it is essential in the present scenario to adopt a 3-D model of the squirmer to investigate the slip-modulated rheotaxis even in the absence of such external perturbations.

Figure 1. Illustration of a spherical microswimmer in a background pure shear flow $\tilde {\boldsymbol {u}}^{(\boldsymbol {ex})}_{\infty } = \dot {\gamma } (\tilde {l}_s+\tilde {z}) \boldsymbol {e}_x$ adjacent to a slippery planar wall satisfying the Navier slip condition and having a slip length $\tilde {l}_s$. The microswimmer is of radius $a$, its orientation vector is designated with the vector $\hat {\boldsymbol {p}}$ and the centre (point $O$) is located at a height of $\tilde {h}$ from the wall. Here, $O x' y' z'$ denote the body-fitted reference frame for the microswimmer. The vector $\hat {\boldsymbol {p}}$ makes an angle $\theta _p$ (also known as the pitching angle) with the wall, while $\phi _p$ is its angle with the external flow direction, measured in the plane of the wall. The flow around the model microswimmers in the laboratory frame, for both puller and pusher types, has been schematically shown as an inset.

2.1. Governing equations and boundary conditions

Due to the low Reynolds numbers encountered in flow around microswimmers, the flow field can be described by the Stokes equation (Lauga & Powers 2009; Michelin & Lauga 2014). In addition to this, application of the incompressibility condition leads to the following governing equation for fluid flow:

(2.1a,b)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \tilde{\boldsymbol{u}}=0 \quad \text{and} \quad -\boldsymbol{\nabla} \tilde{p}+ \mu \nabla ^2 \tilde{\boldsymbol{u}}=0. \end{equation}

Here $\tilde {\boldsymbol {u}}$ denotes the velocity vector and $\tilde {p}$ is pressure.

The hydrodynamic slip velocity $\tilde {\boldsymbol {u}}_{||}$ at the slippery plane wall is related to the shear rate at the wall by the Navier slip condition (Navier 1823) as

(2.2a,b)\begin{equation} \text{at } \tilde{z}= 0, \quad \tilde{\boldsymbol{u}}_{||}= \tilde{l}_{\!_S}\boldsymbol{n}_w\boldsymbol{\cdot}(\boldsymbol{\nabla} \tilde{\boldsymbol{u}}+(\boldsymbol{\nabla} \tilde{\boldsymbol{u}})^{\rm T})(\mathbb{I}-\boldsymbol{n}_w\boldsymbol{n}_w), \end{equation}

where $\boldsymbol {n}_w$ stands for the unit normal at the plane boundary directed towards the fluid domain and $\mathbb {I}$ represents the identity tensor. The presently developed mathematical model remains applicable for an arbitrary magnitude of the slip length at the plane wall $(\tilde {l}_s)$, in stark contrast to the earlier works with asymptotically small slip lengths (Swan & Khair 2008; Willmott 2008). Moreover, the present model is applicable for a wide range of heights of the microswimmer above the wall, starting from an unbounded domain to the lubrication regime. Hence, the results of the present work cannot be obtained just by considering a linear superposition of the slip-induced effects with those of a no-slip problem. Rather, the slip effects are intrinsically coupled with the hydrodynamic problem, the effects of which can only be visualized through a detailed analysis, as performed subsequently.

The linearity property of the Stokes flow (2.1a,b) and the boundary condition at the microswimmer surface allow us to decompose the full flow problem in two sub-problems ‘sq’ and ‘ex’, to be discussed subsequently. Thus, different flow variables ($\psi \in [ \tilde {\boldsymbol {u}}, \tilde {\boldsymbol {V}}, \tilde {\boldsymbol {\varOmega }}]$) can be expressed as

(2.3)\begin{equation} \psi=\psi^{(sq)}+\psi^{(ex)}, \end{equation}

where $\tilde {\boldsymbol {V}}$ and $\tilde {\boldsymbol {\varOmega }}$ denote the translational and rotational velocity components of the microswimmer, respectively.

2.1.1. Sub-problem ‘sq’

Diverse swimming appendages, like cilia or flagella, create surface distortions and work behind the motility of microswimmers. We model this swimming action by the ‘squirmer’ model proposed by Lighthill (1952) and Blake (1971), and extensively used in literature (Uspal et al. 2015; Ishimoto 2017; Poddar et al. 2020) related to self-propelling microswimmers. Accordingly, we impose a tangential surface velocity, given by

(2.4)\begin{equation} \tilde{\boldsymbol{u}}^{(sq)}_s=\left(\frac{\hat{\boldsymbol{p}}\boldsymbol{\cdot} \boldsymbol{r} }{|\boldsymbol{r}|}\frac{\boldsymbol{r}}{|\boldsymbol{r}|}-\hat{\boldsymbol{p}}\right) \sum_{n=1}^{\infty} \frac{2}{n(n+1)} B_n P'_n\left(\frac{\hat{\boldsymbol{p}} \boldsymbol{\cdot} \boldsymbol{r} }{|\boldsymbol{r}|}\right). \end{equation}

Here, $\boldsymbol {r}$ denotes the position vector of points on the microswimmer surface with reference to its centre, $B_n$ denotes the amplitude of the $n{{\rm th}}$ squirming mode and $P'_n({\hat {\boldsymbol {p}}\boldsymbol {\cdot } \boldsymbol {r} }/{|\boldsymbol {r}|})$ is the derivative of the Legendre polynomial with respect to the argument ${\hat {\boldsymbol {p}}\boldsymbol {\cdot } \boldsymbol {r} }/{|\boldsymbol {r}|}$. For an unbounded creeping flow, the first squirming mode contributes solely to the propulsion speed, whereas the second mode quantifies the strength of the stresslet exerted by the squirmer (Ishikawa, Simmonds & Pedley 2006; Li & Ardekani 2014; Chisholm et al. 2016; Pedley 2016; Pietrzyk et al. 2019). Thus, similar to a host of earlier works (Ishikawa et al. 2006; Li & Ardekani 2014; Uspal et al. 2015; Shaik & Ardekani 2017; Shen, Würger & Lintuvuori 2018; Yazdi & Borhan 2017), we retain the first two squirming modes ($B_1$ and $B_2$) only to capture the essential physics of squirming motion. The ratio of the first two squirming modes $\beta$ arises as an important parameter in the problem, and the parameter helps categorizing the different members of the microswimmer family as puller ($\beta >0$), pusher ($\beta <0$) and neutral ($\beta =0$) types. We further adopt a non-dimensionalization scheme based on reference values chosen for different variables as: length $\sim$ radius of microswimmer $(a)$, velocity $v_{ref} \sim B_1$, time $t_{ref} \sim a/B_1$ and pressure $p_{ref} \sim \mu B_1/a$, and subsequently, remove the symbol ‘$\tilde {}$’ to represent the corresponding dimensionless quantities. Thus, the boundary condition at the microswimmer surface can be expressed as

(2.5)\begin{equation} \text{at} \ r=1, \quad \boldsymbol{u}^{(sq)}=\boldsymbol{V}^{(sq)}+\boldsymbol{\varOmega}^{(sq)} \times \boldsymbol{r} + \boldsymbol{u}^{(sq)}_s. \end{equation}

2.1.2. Sub-problem ‘ex’

In this sub-problem we segregate the effects of a background pure shear flow of the form $\boldsymbol {u}^{(ex)}_{\infty } = \mathcal {S} (z+l_s) \boldsymbol {e}_x$ on the locomotion of an inert sphere, disregarding the squirming action. Here, $\mathcal {S}$ denotes the dimensionless shear rate of the background flow, defined as $\mathcal {S} = \dot {\gamma } a / v_{ref}$. The sphere undergoes rigid body motion with velocities $\boldsymbol {V}^{(ex)}$ and $\boldsymbol {\varOmega }^{(ex)}$ in the background flow field. In the presence of the sphere, the disturbed velocity can be written as a superposition of the ambient flow $(\boldsymbol {u}^{(ex)}_{\infty } )$ and the perturbation velocity $(\boldsymbol {u}^{(ex)})$ as $\boldsymbol {u}^{(ex)} + \boldsymbol {u}^{(ex)}_{\infty }$. Since the no-slip condition holds true at the particle surface for the disturbed flow, the perturbation velocity satisfies the following boundary condition:

(2.6)\begin{equation} \text{at} \ r=1, \quad \boldsymbol{u}^{(ex)}=\boldsymbol{V}^{(ex)}+\boldsymbol{\varOmega}^{(ex)} \times \boldsymbol{r} - \boldsymbol{u}^{(ex)}_{\infty}. \end{equation}

The unknown velocity components ($\boldsymbol {V} \ \text {and}\ \boldsymbol {\varOmega }$) are evaluated by considering a neutrally buoyant microswimmer within the flow field, leading to the following force and torque-free conditions:

(2.7a,b)\begin{equation} \boldsymbol{F}= \iint\limits_{S_p} {\boldsymbol{\sigma}}\boldsymbol{\cdot} \boldsymbol{n_p}\, {\rm d}S=0 \quad \text{and} \quad \boldsymbol{L}= \iint\limits_{S_p} \boldsymbol{r} \times ({\boldsymbol{\sigma}}\boldsymbol{\cdot} \boldsymbol{n_p})\, {\rm d}S=0. \end{equation}

Now, the thrust force or torque experienced by the microswimmer originates from the squirming action $(\boldsymbol {F}^{(sq)}_{(Thrust)},\boldsymbol {L}^{(sq)}_{(Thrust)})$ as well as the externally applied flow $(\boldsymbol {F}^{(ex)}_{(Thrust)},\boldsymbol {L}^{(ex)}_{(Thrust)})$. The resultants of these thrusts are further balanced by the hydrodynamic resistance on the rigid sphere $(\boldsymbol {F}_{(Drag)},\boldsymbol {L}_{(Drag)})$. Again exploiting the linearity of the problem, two components of the thrust force can be linearly added, reducing (2.7a,b) to

(2.8a,b)\begin{equation} \boldsymbol{F}^{(sq)}_{(Thrust)}+\boldsymbol{F}_{(Thrust)}^{(ex)}+\boldsymbol{F}_{(Drag)}=0 \quad \text{and} \quad \boldsymbol{L}^{(sq)}_{(Thrust)}+\boldsymbol{L}_{(Thrust)}^{(ex)}+\boldsymbol{L}_{(Drag)}=0. \end{equation}

2.2. Solution strategy

In order to solve the above system of governing equations and boundary conditions, together with the force-free constraint, we use eigenfunction expansion of the Stokes flow problem in the bispherical coordinates $(\xi, \eta, \phi )$ (Happel & Brenner 1983). In the bispherical system the plane boundary is located at $\xi =0$ and the spherical swimmer surface corresponds to $\xi =\xi _0$ (Behera, Poddar & Chakraborty 2023; Poddar 2023). In this solution method the expressions for the velocity components contain a set of unknown coefficients $A_n^m,B_n^m, C_n^m, E^m_n, F_n^m, G_n^m$ and $H^m_n$ (details in Appendix A). To solve for these constants, we employ the different boundary conditions ((2.2a,b), (2.5) and (2.6)), the no penetration condition at the solid surfaces and the incompressibility condition, and apply the orthogonality property of the Legendre polynomials. Due to the decaying nature of these constants, we truncate the infinite series solution of the flow field at large values of $N$ that give an accuracy of $\boldsymbol {O}(10^{-6})$ between successive values of each of the constants considered. The linear algebraic equations to be solved simultaneously for the unknown coefficients have been arranged as a banded matrix of size $7N\times 7N$, to be solved numerically. Higher values of $l_s$ causes convergence issues and an increasingly higher number of terms have to be retained before they are solved numerically. It is to be noted that the simplicity of a no-slip boundary condition at the plane wall allows one to explicitly relate all the other constants in terms of a single constant (O'Neill 1964) and only a matrix of size $N \times N$ has to be inverted to obtain all the $7N$ desired unknown constants. In contrast, the slip boundary condition complicates the numerical task by demanding a matrix inversion of size $4N\times 4N$. Here, the conversion of an original $7N \times 7N$ system to a smaller matrix size is performed by considering the following compatibility condition (Loussaief et al. 2015):

(2.9)\begin{equation} \frac{\partial u_z}{\partial z} = l_s \frac{\partial^2 u_z}{\partial z^2}. \end{equation}

The constants $X_n^m,Y_n^m$ and $Z_n^m$ (A15)–(A16), associated with the boundary conditions on the swimmer surface, take different forms for different sub-problems. They were derived by Shaik & Ardekani (2017) in relation to the ‘sq’ sub-problem. We derive the corresponding constants for the ‘ex’ sub-problem as

(2.10a–c)\begin{align} X_n^1=0, \quad Y_n^1={-}2 \sqrt {2} \, \mathcal{S}\sinh \left( \xi_0 \right) \left( 2n+1 \right) \exp({-( n+1/2 ) \xi_0}) \quad \text{and} \quad Z_n^1=0. \end{align}

The axisymmetry of the truncated squirmer surface velocity model (2.4) confines the orientation vector $(\hat {\boldsymbol {p}})$ of the body-fitted $x'z'$ plane as shown in figure 1. Consequently, its rotation is fixed along the $y'$ axis. Thus, the squirmer velocities can be written as

(2.11a,b)\begin{equation} \boldsymbol{V}^{(sq)}= V^{(sq)}_{z'} \boldsymbol{e}_{z'} + V^{(sq)}_{x'} \boldsymbol{e}_{x'} \quad \text{and} \quad \boldsymbol{\varOmega}^{(sq)}= \varOmega^{(sq)}_{y'} \boldsymbol{e}_{y'}. \end{equation}

In contrast, the symmetry of the pure shear flow imparts a velocity to the inert sphere parallel to the flow direction only, i.e. $\boldsymbol {V}^{(ex)}= V^{(ex)}_{x} \boldsymbol {e}_{x}$, while the rotational motion of the sphere is triggered along the vorticity direction, i.e. $\boldsymbol {\varOmega }^{(ex)}= \varOmega ^{(ex)}_{y} \boldsymbol {e}_{y}$. Thus, the resultant microswimmer velocities in the fixed frame take the form

(2.12a)$$\begin{gather} \boldsymbol{V}= (V^{(sq)}_{x'} \cos(\phi_p) + V^{(ex)}_x )\boldsymbol{e}_{x} + V^{(sq)}_{x'} \sin(\phi_p)\boldsymbol{e}_{y} + (V^{(sq)}_{z'} + V^{(ex)}_{z})\boldsymbol{e}_{z}, \end{gather}$$

(2.12b)$$\begin{gather}\boldsymbol{\varOmega}={-}\varOmega^{(sq)}_{y'}\sin(\phi_p) \, \boldsymbol{e}_{x} + (\varOmega^{(sq)}_{y'}\cos(\phi_p)+\varOmega_{y}^{(ex)}) \boldsymbol{e}_{y}. \end{gather}$$

The details of the ‘squirmer’ thrust components (sub-problem ‘sq’), as well as the hydrodynamic resistance factors (common to both sub-problems ‘sq’ and ‘ex’), can be found in Poddar et al. (2020). The thrust components associated with the external flow (sub-problem ‘ex’) have been evaluated by solving the Stokes problem of pure shear flow around a fixed sphere (i.e. only a part of sub-problem ‘ex’), i.e. considering the boundary condition

(2.13)\begin{equation} \text{at} \ r=1, \quad \boldsymbol{u} ={-} \boldsymbol{u}^{(ex)}_{\infty}. \end{equation}

The thrust force and torque due to the external flow are thus obtained as

(2.14a)\begin{equation} F^{(ex)}_{(Thrust,x)}={-}{\sqrt{2}{\rm \pi}}\,\mathcal{S}(l_s+\cosh(\xi_0))\sinh(\xi_0)\sum_{n=0}^{\infty} [G_n^1-H_n^1+n(n+1)(A_n^1-B_n^1)], \end{equation}

(2.14b)\begin{align} L^{(ex)}_{(Thrust,y)} &={\sqrt{2}{\rm \pi}}\,\mathcal{S}\sinh^2(\xi_0)\sum_{n=0}^{\infty} [\coth(\xi_0)\{n(n+1)(A_n^1-B_n^1)+(G_n^1-H_n^1)\} \nonumber\\ &\quad -2n(n+1)C_n^1 -(2n+1)(G_n^1-H_n^1)]. \end{align}

It is to be noted that the thrust due to the external flow can be alternatively calculated by using the Lorentz reciprocal theorem (LRT) as outlined in Poddar et al. (2020) for motion near a slippery plane, thus bypassing the solution to the full Stokes problem. We use LRT only to verify the results obtained by the full solution technique, as presented above.

2.3. Swimming trajectories

The quasi-steady dynamics of the microswimmer (Spagnolie & Lauga 2012; Uspal et al. 2015; Mozaffari et al. 2016; Walker et al. 2018) can be fully described by simultaneously determining the location of the microswimmer in space $\boldsymbol {r} (t)$ along with its preferential orientation with respect to the plane wall, represented by $\hat {\boldsymbol {p}} (t)$. Thus, the trajectories can be obtained by simultaneously solving the following set of coupled ordinary differential equations:

(2.15a,b)\begin{equation} \frac{{\rm d}\boldsymbol{r}(t)}{{\rm d}t}=\boldsymbol{V} (\boldsymbol{r}(t),\hat{\boldsymbol{p}}(t))\quad \text{and} \quad \frac{{\rm d}\boldsymbol{\hat{p}(t)}}{{\rm d}t}=\boldsymbol{\varOmega}(\boldsymbol{r}(t),\hat{\boldsymbol{p}}(t)) \times \hat{\boldsymbol{p}}(t) \end{equation}

for a given set of initial conditions $(\boldsymbol {r}_0,\hat {\boldsymbol {p}}_0)$. The different translational and rotational velocity components at each time instant can be obtained by a combination of their self-propulsion and external flow counterparts, as discussed in (2.12). We neglect the effects of stochastic forces on the microswimmer motion and compute the trajectories considering the deterministic forces only (Shum et al. 2010; Spagnolie & Lauga 2012; Mozaffari et al. 2016; Poddar, Bandopadhyay & Chakraborty 2021).

3. Results and discussion

In this section we illustrate the combined interaction of wall slip and a background shear flow in dictating the locomotion characteristics of both puller- and pusher-type microswimmers. The dimensionless analysis presented above can fully describe the rheotactic swimming near a slippery plane using the parameters dimensionless shear rate $\mathcal {S}$, slip length $l_s$ and squirmer parameter $\beta$; in addition to the positional $(\boldsymbol {r})$ and orientational variables of the microswimmer $(\hat {\boldsymbol {p}})$. In the following subsections we discuss the 3-D trajectories under the influence of different dimensionless parameters involved. Subsequently, we summarize these effects in the form of regime maps and illuminate on the governing physics behind contrasting motion characteristics.

In order to estimate the practical values of the parameter $\mathcal {S} = \dot {\gamma } a / v_{ref}$, we consider realistic ranges of the dimensional quantities in different microfluidic experiments related to microswimmers (Kantsler et al. 2014; Ohmura et al. 2021), as stated below: shear rate $\dot {\gamma } = 0.1\ {\rm s}^{-1}$ to $20\ {\rm s}^{-1}$, velocity of a typical microswimmer $v_{ref} = 10 \ \text {to} \ 100 \ \mathrm {\mu } {\rm m}\ {\rm s}^{-1}$ and the length scale of the microswimmer $a = 10 \ \text {to} \ 100 \ \mathrm {\mu } {\rm m}$. Although $\mathcal {S}$ ranges in $\boldsymbol {O}(10^{-2}\unicode{x2013}10^{2})$, a high value of the same parameter amounts to sweeping away of the microswimmer along the external flow. Consequently, the competitive effects of the shear flow and self-propulsion remain obscure. Thus, motivated by the earlier theoretical investigations (Uspal et al. 2015; Walker et al. 2018), we choose $\mathcal {S}$ between 0 to 1. Similarly, considering previous experimental observations, the dimensionless slip length $(l_{s})$ is varied between 0 and 10 (Zhu & Granick 2001; Tretheway & Meinhart 2002; Huang et al. 2008).

3.1. Swimming states for puller microswimmers

Here, we investigate the modulations in the rheotactic states of puller microswimmers brought in by the near-wall hydrodynamic slip. Consideration of only hydrodynamic forces can resolve the microswimmer dynamics only upto a finite gap with the solid flat surface (Shum et al. 2010; Spagnolie & Lauga 2012; Uspal et al. 2015) due to the requirement of infinite computational resources. Thus, the present trajectory simulations based on forces of pure hydrodynamics origin have been performed by considering a minimum distance between the microswimmer surface and the wall as $\delta = 0.01$. Consequently, any swimming state indicating a downward descend below this gap is considered as a ‘crashing’ or ‘collision’ state (Uspal et al. 2015). This mathematical treatment is justified since below this small gap, nanoscale interaction forces other than a hydrodynamic origin (Klein, Clapp & Dickinson 2003) become prominent and are expected to influence the motion characteristics. However, to facilitate a direct comparison of the present results with the previously reported no-slip cases (Uspal et al. 2015), we have not considered any non-hydrodynamic repulsive force at the plain surface. On the other hand, a microswimmer going beyond a height of $h= 15$ marks its ‘escape’ from the wall, similar to earlier works (Ishimoto & Gaffney 2013; Poddar et al. 2020).

The 3-D trajectories of the microswimmer moving in a background pure shear flow with or without near-wall slippage for different initial orientations ($\theta _0$) have been compared in figure 2. Figure 2(a) demonstrates that near a no-slip wall the microswimmer swims against the flow with damped amplitude oscillations in the vertical direction and finally shows a sliding motion after reaching an out-of-plane angle $\phi _p=83^{\circ }$, keeping a constant height $h=1.23$ and orientation $\theta _p=27.5^{\circ }$. These motion characteristics have been termed as ‘upstream rheotaxis’ in the literature (Uspal et al. 2015; Ishimoto 2017; Walker et al. 2018).

Figure 2. Trajectories of a puller microswimmer $(\beta = 7)$ in background pure shear flow ($\mathcal {S}=0.1$), starting from the same set of initial conditions $x_0=0, y_0=0, h_0=2, \phi _0 = 270^{\circ }$ but different pitch angles, (a,b) $\theta _0=10^{\circ }$ and (c,d) $\theta _0=150^{\circ }$. The effects of the slippery $(l_s=1)$ and no-slip boundary conditions on the trajectories are studied. Red solid trajectories are 3-D representations of motion. Blue and black dashed lines are the projections of the 3-D trajectories on the $xz$ and $yz$ planes, respectively. Arrowheads represent the directions of motion of the microswimmer.

In stark contrast, in the presence of wall slip $(l_s=1)$, the microswimmer comes in close proximity to the wall and descends below $\delta = 0.01$, which indicates a ‘collision’ state in the absence of an additional contact force (shown in figure 2b). It is worth mentioning that the cutoff distance for trajectory simulations leaves the possibility for the mathematical model to predict more stable swimming states against the flow (rheotaxis) instead of the collision states presented here if the numerical simulations were performed for wall gaps below $\delta =0.01$. Thus, the collision states predicted through a pure hydrodynamic analysis may not be observed if a non-hydrodynamic repulsive interaction at the wall is considered.

However, an altered orientation angle ($\theta _0=150^{\circ }$) results in an escaping trajectory for the no-slip case, as presented in figure 2(c). Near the boundary, the microswimmer experiences a strong counterclockwise (CCW) reorientation torque due to hydrodynamic interaction with the wall. This enhanced torque has a tendency to rotate the director away from the wall and the vertical velocity switches from $V_z<0$ to $V_z>0$ at a point, resulting in a rapid reorientation of the microswimmer.

The collision states created by the wall slip in figures 2(b) and 2(d) suggest that the said torque is weakened by the slip effects, thus failing to supply the required reorientation for rheotactic sliding or escape. Further comparing the two collision states in the presence of the wall slip, it is found that collision occurs much earlier with $\theta _0=10^{\circ }$ (at $t_{c}=2.8$) than $\theta _0=150^{\circ }$ (at $t_{c}=12$), despite the director $\hat {p}$ initially tilting more towards the wall in the latter scenario. Here, $t_{c}$ refers to the collision time, which is the duration until the microswimmer descends to $\delta =0.01$ within the flow field. This is due to the contrasting consequences of the slip-induced torque in the two situations. While for $\theta _0={10^{\circ }}$, the said torque assists rotation towards the wall, it favours a director movement away from the wall for $\theta _0={150^{\circ }}$. However in the latter case, this slip-triggered torque is not sufficient for escaping or rheotactic sliding since the cutoff height $\delta =0.01$ has already been encountered.

A comprehensive understanding of the non-trivial motion characteristics due to the diverse plausible combinations of the parameters $l_s$, $\beta$ and initial conditions demands a large number of individual long-time trajectory simulations in three dimensions, calling for a massive computational time. However, the analysis can be greatly simplified by using the theory of dynamic systems and considering the symmetries in the system. In this regard, we align the plane of microswimmer motion $(x'z')$ along the plane of the background flow $(xz)$, and thus preventing rotation out of the shearing plane. Hence, the pitch angle $\theta _p$ fully parametrizes the angular orientation of the director $\hat {\boldsymbol {p}}$, without loss of generality. As a result, the dynamics of the microswimmer can be described by a plane autonomous system, i.e.

(3.1a,b)\begin{equation} \frac{{\rm d}h(t)}{{\rm d}t}=V^{(sq)}_{z'} + V^{(ex)}_{z} \quad \text{and} \quad \frac{{\rm d}\theta_p(t)}{{\rm d}t}=\varOmega^{(sq)}_{y'}+\varOmega_{y}^{(ex)}. \end{equation}

The relevance of the above dimensionally reduced system in predicting the behaviour of the full system was thoroughly examined by Walker et al. (2018), and long-time simulations of the full system (2.15) were found to be in accordance with the restricted system.

3.1.1. Annihilation of rheotactic states

Here we provide a concise representation of puller dynamics by analysing the phase portraits obtained from (3.1), the results of which are shown in figures 3(a)–3( f). The justification behind the consideration of the phase-plane dynamics instead of analysing all possible out-of-plane dynamics has been discussed in Appendix B. The results corresponding to the no-slip wall in figures 3(a) and 3(d) are in perfect agreement with the work of Uspal et al. (2015). As shown in figure 3(a), in a quiescent environment $(\mathcal {S}=0)$, two stable dynamical attractors (black square markers) appear between $\theta _p=0^{\circ }$ and ${180^{\circ }}$ at mirror symmetric locations about $\theta _p={90^{\circ }}$. These points indicate the final swimming states of the microswimmer as sliding at a steady height and fixed orientation but in opposite directions along the $x$ axis. However, the invariance of the dynamic system along $x$ in the absence of a background flow leads to the same sense of these swimming states, similar to the reported results of Uspal et al. (2015). This stable swimming state is a sole consequence of the propulsive torque generated beyond a critical value of $\beta$ (Li & Ardekani 2014; Poddar et al. 2020). In addition, two unstable fixed points (red triangles) and a saddle point (denoted with an orange cross) are observed in the same phase portrait. The slip condition at the wall $(l_s>0)$ severely influences the dynamics, as can be found by comparing figures 3(b) and 3(c) with figure 3(a). For a low value of slip length $l_s = 0.36$, the unstable fixed points disappear from the phase portrait and the steady-state height $(h^*)$ corresponding to the attractors comes downward. While this trend of downward shifting fixed points retained for higher slip lengths (e.g. $l_s=1$ in figure 3c), the microswimmer gradually descends below a height $\delta = 0.01$, and finally, collides against it, thus wiping out the attractors from the phase portrait. As a concurrent effect, the unstable fixed points are also suppressed.

Figure 3. Phase space diagrams ($\theta _p \text {-} h$ plot) for the reduced dynamical system in the $x$–$z$ plane for a puller microswimmer $(\beta = 7)$. Absence of external flow presented in (a–c) for $l_s=0$ , $0.36$ and $1$ slip lengths. The dimensional shear rate has been chosen as $\mathcal {S}=0.1$ in (d–f) along with the same environment for the microswimmer. Green circular markers indicate the rheotactic attractor whereas black squares imply the non-rheotactic attractor that is unstable for $\phi \neq 0$. (g) Range of the rheotaxis zone is shown for different slip lengths before annihilation by the plotting vertical distance ($h^*$) of the rheotactic attractor for different slip lengths. (h) Behaviour of the saddle point (near $\theta ^*>270^{\circ }$) is plotted for different slip lengths $l_s$. The green and red superimposed lines in (d, f), respectively, denote the sample phase space trajectories starting from $h_0=2$ and $\theta _{p,0} =153.8^{\circ }$.

The physics behind the above observations can be described by analysing the time variations of different velocity components of the microswimmer, as shown in figure 4. Figure 4(c) shows that the rotational velocity due to squirming action ($\varOmega _y^{(sp)}$) has a clockwise (CW) magnitude at the initial times, leading to a rotation of the director $\hat {\boldsymbol {p}}$ towards the wall. A further illustration of figure 4(d) reveals that with the increase in time, the CW magnitude strengthens for a slip length of $l_s=1$. On the other hand, consideration of the other source of microswimmer rotation, i.e. the background shear flow, reveals that $\varOmega _y^{(ex)}$ remains unaffected by the slip length in the far field $(\delta \to \infty )$ and attains a constant magnitude $\varOmega _y^{(ex)}=\mathcal {S}/2$. However, $\varOmega _y^{(ex)}$ becomes a function of $l_s$ in the wall-adjacent region. It was shown that the magnitude of $\varOmega _y^{(ex)}$ decreases with increasing slip length, with a maximum change (for $\delta = 0.01$) of $3.7\, \%$ and $11.7\, \%$ for $l_s=1$ and 10, respectively (Loussaief et al. 2015).

Figure 4. Time variation of different microswimmer velocity components corresponding to the sample trajectories highlighted in figure 3. Panels (a,c) correspond to no-slip results, while (b,d) correspond to the case with $l_s=1$. In each panel, the contributions from shear flow and squirming action are also presented. Panel (e) shows the behaviour of $V_z$ for $l_s=0$ (dashed green line) and for $l_s=1$ (red line). The blue dashed line indicates the collision time of the microswimmer.

The velocity component in the vertical direction $V_z$ remains unaffected by the background shear flow. In the absence of wall slip, the vertical velocity is less reduced as compared with the $l_s=1$ condition (see figure 4d), which retains the microswimmer at a greater height at the collision time $`t_c$’. A simultaneous strong CCW rotation of the director $(\varOmega _y>0)$ lifts off the microswimmer from the collision zone and subsequently imparts a vertically upward velocity $(V_z>0)$. After the microswimmer attains a certain height, the transition from CCW to CW rotation takes place, which facilitates switching of the vertical motion $(V_z>0\ \text {to}\ V_z<0)$. This cycle continues with a damped amplitude of oscillations and finally leads to the rheotactic attractor. In comparison, stronger magnitudes of the CW rotation $(\varOmega _y<0)$ and downward movement $(V_z<0)$ results in the presence of wall slip. Therefore, at time $t_c$, the height from the wall $(h)$ is also reduced under the action of slip, and, as a consequence, the microswimmer does not face enough force required to run away from the trapped condition of crashing against the wall.

Although the velocity component of the microswimmer parallel to the wall $V_x$ does not directly influence the phase portraits, the time variation of the microswimmer position along its trajectory is highly dependent on it. In the absence of shear flow, the wall slip causes drastic changes in the near-wall self-propulsion velocity of the squirmer $(V_x^{(sq)})$, with $V_x^{(sq)}$ becoming higher or lower than the far-field velocity $V_x^{(sq)}|_{z\to \infty } = \cos (\theta _p)$, depending on the influence of the wall slip on the propulsive thrust and the resistance factors (Poddar et al. 2020). Now, a background shear flow always contributes a positive velocity $V_x^{(ex)}$ for $\mathcal {S}>0$ irrespective of the orientation of the director or the slip length. However, the magnitude of $V_x^{(ex)}$ is highly dependent on slip length in the near-wall zone (Loussaief et al. 2015). It is found that the magnitude of $V_x^{(ex)}$ enhances with a corresponding increase in slip length. For example, the said enhancement is $257.62\,\%$ for $l_s=1$ and $1989.68\,\%$ for $l_s=10$ at the minimum wall gap considered, i.e. $\delta \to 0.01$. Thus, different combinations of $l_s$ and $\mathcal {S}$ would cause drastic modulations in the overall velocity $V_x$ of the microswimmer up or against the flow, the direction of which is decided by the relative importance of the $V_x^{(sq)}$ and $V_x^{(ex)}$ for the situation under consideration. For the case demonstrated in figure 4(a), the sign of $V_x$ remains unchanged $({<}0)$ under the additional shear component $V_x^{(ex)}$ but faces significant changes in magnitude and modulation of the amplitudes of damped oscillations in $V_x$.

The inclusion of background shear flow retains the attractors for puller dynamics at almost the same locations, but the dynamic properties of these attractors become contrasting in nature (compare figures 3a and 3d). The attractor located in the region $90^{\circ }<\theta _p<180^{\circ }$ (green dot in figure 3d) was shown to be a global attractor displaying stability even against out-of-plane $(xz)$ perturbations (Uspal et al. 2015). This attractor represents a sliding state of the microswimmer with its orientation vector directed against the background flow and towards the wall. This motion behaviour of the microswimmer is known as the upstream rheotaxis. In stark contrast, the attractor at $\theta _p<90^{\circ }$ is unstable for out-of-plane perturbations ($\phi \neq 0$) and remains stable for in-plane perturbations only (Uspal et al. 2015). The rheotactic attractor comes with an accompanying saddle point just above it (black cross-marker in the figure). In addition, the existing saddle point, which appears at ${270^{\circ }}$ for $\mathcal {S}=0$, shifts towards the right of the phase portrait due to the action of the shear flow. For example, with $\mathcal {S}=0.1$, the orientation angle corresponding to this saddle point shifts from $270^{\circ }$ to $283^{\circ }$.

The slip-induced modulations to the near-wall rheotaxis discussed above are also highlighted in figures 3(e) and 3( f) for two different values of slip length $l_s=0.36$ and 1, respectively. It is found that the wall slip causes a downward shift of the rheotactic attractor and finally leads to its annihilation, keeping its angular location unaffected. The effect of increasing slip length on the height of the rheotactic attractor $h^*$ has been portrayed in figure 3(g). It depicts that $h^*$ decreases monotonically until a threshold value of the slip length $l_s \approx 0.5$ is reached and causes the disappearance of the rheotactic state beyond this threshold. It is interesting to note that the shear-induced saddle point above the rheotactic attractor hardly faces any change in its location in the phase portrait for non-zero slip lengths. However, the shear-induced rightward shifting of the other saddle is suppressed with increased slip length, as shown in figure 3(h). It can also be observed from the same figure that the saddle point orientation saturates to a value of ${\approx }273^{\circ }$ for increasing slip lengths. Furthermore, similar to the attractors, this saddle point also moves downward.

3.1.2. Emergence of new rheotactic states

Here, we elaborate on the slip-induced transitions in puller dynamics at a higher shear strength $\mathcal {S} = 0.15$. The phase portraits in figures 5(a) to 5(c) highlight three different cases of slip length $l_s=0$, 0.137 and 1, respectively. In addition, figure 5(d) displays three contrasting in-plane trajectories exhibited by microswimmers starting from the same initial condition $h_0=2$ and $\theta _0=150.69^{\circ }$ but for varying slip lengths at the wall. The corresponding phase space trajectories $(\theta _p \ {\rm vs}\ h)$ have also been superimposed on the phase portraits. Unlike the situation presented for $\mathcal {S}=0.1$ in § 3.1.1, upstream rheotaxis is not observed for a no-slip wall when $\mathcal {S}=0.15$. This phenomenon can be substantiated by considering the relative importance of the two driving mechanisms involved, namely, the squirming action and the background flow. The term proportional to $\beta$ in the expression of squirmer surface velocity (2.4) represents the intensity of vorticity generated by the self-propulsion mechanism. It was identified that the hydrodynamic rotation of a squirmer in wall proximity is considerable only when $\beta$ is sufficiently high (Ishimoto & Gaffney 2013). Thus, increasing the shear strength $(\mathcal {S})$ calls for a higher value of the squirmer parameter (here $\beta >7$) so that the necessary competitive behaviour of the two flows can be achieved. As a result, the rotational dynamics imparted by the shear flow dominate, leading to the disappearance of the rheotactic attractor, and the microswimmer escapes from the near-wall region.

Figure 5. Phase space diagrams for (a) $l_s=0$, (b) $l_s=0.14$ and (c) $l_s=1$. (d) In-plane trajectories ($x$ vs $z$) of the microswimmer for these slip lengths when released from location $h_0=2$ and $\theta _0=150.7^{\circ }$. The corresponding phase space trajectories have been highlighted in (a–c). Other parameters $\beta =7$ and $\mathcal {S}=0.15$ are common for all the panels. Initial launching conditions are highlighted with circular yellow dots. Note that we have shown the phase portraits only in the region $90^{\circ } \le \theta _p \le 270^{\circ }$ and $h \le 4$ in order to focus on the rheotactic attractor.

Although the interaction between the shear flow and squirming action in this case is not sufficient to result in a rheotactic attractor, the said interaction is non-trivially intervened by the wall slip, and a new rheotactic attractor comes into existence. This outcome demonstrated for a low value of slip length $l_s=0.14$ in figures 5(d) and 5(b), constitutes the most counterintuitive result of the present work. In this case, the microswimmer reaches a steady sliding state opposite to the flow direction with a fixed height $h^*=1.12$ and orientation $\theta _p^*=153.46^{\circ }$. The saddle point accompanying the rheotactic attractor has been highlighted with a black cross in 5(b). While annihilation of the rheotactic attractor can be attributed to the mathematical limitation of predicting motion behaviour below a cutoff height $(\delta =0.01)$, the ‘emergence’ is a novel physical consequence of the combined interplay between slip and shear.

In order to bring out the contrasting physics of the in-plane dynamics with or without the wall slip, we appeal to the sample trajectories highlighted in figure 5(d). For the no-slip case, the combination of the squirming action and shear flow results in a strong CCW rotation throughout the trajectory of the microswimmer. This causes the director $\hat {\boldsymbol {p}}$ to point away from the wall even if the swimmer is very close to the wall. As an ensuing effect, the downward velocity $(V_z<0)$ also reduces in magnitude, and switching in the direction of the vertical motion sets in at $\theta _p \approx {153^{\circ }},\ h=1.21$. The lack of a mechanism capable of triggering an opposing propulsive torque causes the microswimmer to finally escape from the wall. Now, in the presence of wall slip, both the propulsive torques due to the squirming action and shear flow are greatly modified. As a net effect, the microswimmer gains a tendency to rotate its director in the opposite direction (i.e. CW), with an accompanying downward movement of the microswimmer. Subsequently, a switching of vertical motion takes place, followed by an oscillatory motion with damped amplitude, leading to a stable rheotactic state, similar to the discussion in § 3.1.1.

With further enhancement in the slip length $(l_s)$, however, the height $(h^*)$ of the rheotactic attractor reduces, and finally, a collision state results when a critical slip length is exceeded (see figure 5c). This result unravels the fact that the slip-induced rheotactic attractors are observed only for some intermediate range of slip lengths pertaining to a specific value of the dimensionless shear rate $\mathcal {S}$.

The slip-induced modulations to the near-wall rheotaxis discussed above are also highlighted in figures 3(e) and 3( f) for two different values of slip length $l_s=0.36$ and 1, respectively. It is found that the wall slip causes a downward shift of the rheotactic attractor and finally leads to its annihilation, keeping its angular location unaffected. The effect of increasing slip length on the height of the rheotactic attractor $h^*$ has been portrayed in figure 3(g). It depicts that $h^*$ decreases monotonically until a threshold value of the slip length $l_s \approx 0.5$ is reached and causes the disappearance of the rheotactic state beyond this threshold. It is interesting to note that the shear-induced saddle point above the rheotactic attractor hardly faces any change in its location in the phase portrait for non-zero slip lengths. However, the shear-induced rightward shifting of the other saddle is suppressed with increased slip length, as shown in figure 3(h). It can also be observed from the same figure that the saddle point orientation saturates to a value of ${\approx }273^{\circ }$ for increasing slip lengths. Furthermore, similar to the attractors, this saddle point also moves downward.

3.1.3. Summary of slippery rheotaxis for pullers

In order to obtain a comprehensive understanding of the combined effects of the key dimensionless parameters, i.e. slip length $l_s$ and shear rate $\mathcal {S}$, we analyse the in-plane phase portraits described by (3.1) for all possible combinations of $l_s$ and $\mathcal {S}$ in the regime $0\le l_s \le 10$ and $0 \le \mathcal {S} \le 1$, and the results are summarized in figure 6. In this figure we show regime maps demarcating different motion characteristics for a wide range of values of positive $\beta$.

Figure 6. Regime maps demarcating different motion characteristics of puller microswimmers with (a) $\beta =3$, (b) $\beta = 4.5$, (c) $\beta =6$, (d) $\beta =7$ and (e) $\beta =30$. The filled green circles represent parametric combinations $(l_s, \mathcal {S})$ giving rise to upstream rheotaxis, while the grey triangular markers correspond to non-rheotactic states, i.e. either collision or escape. The blue dashed lines mark the critical dimensionless shear rates $\mathcal {S}_{cr}$ (${=}\,0.04$, $0.1$ and $0.14$ for (b), (c) and (d), respectively) required for achieving a rheotactic state caused exclusively by the presence of slip. The red dots indicate the mirror symmetric sliding states in the absence of shear flow. Also, A and B zones in each figure denote two types of non-rheotactic zones.

In the absence of the background flow, mirror symmetric sliding states (shown with red dots) are observed, as described earlier in figures 3(a) and 3(b). The transition from the red to grey zone corresponding to the $\mathcal {S}=0$ case signifies the suppression of attractors by wall slip, the details of which have been shown in figure 3(c). One of these mirror symmetric sliding states, located in the pitch angle range $90^{\circ }\leq \theta _p \leq 180^{\circ }$, transforms into the upstream rheotactic state upon strengthening of the external shear $(\mathcal {S} > 0)$. Thus, a band of rheotactic states is developed in the $(l_s, \mathcal {S})$ map, dividing two types of non-rheotactic zones (marked as A and B in the regime maps) in figure 6. Zone B corresponds to suppression of the rheotactic attractors and occurrence of collision states due to high slip length. Zone A develops above a critical shear rate $(\mathcal {S}_{cr})$, where external shear washes out the rheotactic states to result in escape states in the absence of adequate wall slip. The influence of the squirmer parameter $\beta$ can be observed by comparing the regime maps for low to high values of the same parameter. While in the limit of a low value of $\beta$ (e.g. $\beta = 3$ in figure 6a) zone B is absent, the other extreme condition of very high $\beta$ (e.g. $\beta = 30$ in figures 6e) does not exhibit the characteristics of zone A. The sequential appearance/disappearance of the non-rheotactic zones A and B as well as the different spans of the rheotactic zone in these panels provide a general understanding of the impact of $\beta$ on the rheotactic states in the presence of wall slip.

It is noteworthy that below $\mathcal {S}_{cr}$ the rheotactic attractors are present at lower values of slip length for intermediate magnitudes of $\beta$ (see figure 6b–d), even for a no-slip condition. Here, the insufficient strength of the background shear flow is unable to wash out the steady-state motion of the microswimmer. However, upstream rheotaxis is reported for higher values of the slip length above the critical shear rate (e.g. for $\beta =7$, at $\mathcal {S}=0.15$ the squirmer exhibits upstream steady rheotactic motion at $l_s= 0.07$). Thus, increasing the shear rate shifts the rheotactic zone towards the right, which implies that strengthening the external shear requires an intense hydrodynamic interaction between the squirming action and the boundary slip for rheotaxis to occur.

The critical shear rates $\mathcal {S}_{cr}$ increase with the squirmer parameter $\beta$. For instance, a comparative analysis of figures 6(b)–6(d) brings out that $\mathcal {S}_{cr}$ is raised from 0.04 for $\beta =4.5$ to 0.14 for $\beta =7$. Moreover, a lower range of the squirmer parameter widens the span of the rheotactic zone along the $l_s$ axis.

3.2. Swimming states for pusher microswimmers

The phase portrait of a pusher, in the absence of a non-hydrodynamic repulsive force at the wall, can be obtained by simply taking a mirror image of the phase portrait of a corresponding puller and reversing the streamlined directions (not shown for brevity). This is due to the phenomenon of puller–pusher duality associated with time reversal (Ishimoto 2017; Poddar et al. 2020). Similar to Ishimoto (2017), we find a fixed point $(\theta ^*_p,h^*)$ for a puller with a squirmer parameter $\beta$ also indicating the existence of a fixed point at $(-\theta ^*_p,h^*)$ for a corresponding pusher with a squirmer parameter $-\beta$. However, the stability property of the fixed point is inverted for the pusher. This indicates the absence of any stable swimming state in the case of pusher microswimmers, reaffirming the observations of Walker et al. (2018) related to virtual Leishmania promastigotes. In contrast to the scenario depicted above, stable swimming trajectories were reported for pusher swimmers under the action of short-range, non-hydrodynamic repulsive interactions with the nearby wall (Spagnolie & Lauga 2012; Lintuvuori et al. 2016; Walker et al. 2019; Poddar et al. 2020). The distinctive near-wall swimming characteristics induced by the additional repulsive interactions irrespective of the background flow were termed as thigmotaxis (Ishimoto 2017).

In an effort to investigate the combined effects of background flow and the slip condition on the contrasting pusher trajectories, we incorporate the following form of repulsive force (Spagnolie & Lauga 2012): $\boldsymbol {F}_{rep}= ({\alpha _1 \exp {(-\alpha _2 \delta )}}/{1-\exp {(-\alpha _2 \delta )}}) \boldsymbol {e}_z$, where $\alpha _1=100,\alpha _2=100$ have been chosen to ensure a minimum wall-swimmer distance of $\delta = 0.01$. The repulsive nature of wall breaks the puller–pusher duality (Ishimoto 2017), and the pusher trajectories can no longer be predicted from the puller trajectories presented in § 3.1. This additional repulsive force has the effect of stabilizing the motion for both pullers and pushers, and thus, stable fixed points can be observed for both swimmer types. The pushers also show the unique trend of attraction towards the wall even if the director is initially tilted away from the wall $(\theta _0>180^{\circ })$ for both no-slip (Lintuvuori et al. 2016) and slip conditions (Poddar et al. 2020).

Pusher microswimmers have been found to exhibit three types of trajectories – escape, upstream oscillations with decaying amplitudes (attractors) and periodic oscillations (limit cycles) in the upstream direction for different combinations of shear rate $\mathcal {S}$ and slip length $l_s$. The regime map in figure 7(a) describes the different swimming states depicted by a pusher ($\beta =-7$) for different plausible combinations of the parameters $\mathcal {S}$ and $l_s$ when launched from the same initial state $\theta _0 =187^{\circ }, h_0= 2$. Here, the upstream rheotaxis (damped oscillations) and undamped oscillations have been numerically distinguished by setting a tolerance of 5 % on the amplitude difference between the first and last oscillations. In the absence of a background shear ($\mathcal {S} =0$), the brown markers denote the mirror symmetric damped (for lower $l_s$) and periodic (for higher $l_s$) oscillatory features of the trajectories.

Figure 7. (a) Regime map of a pusher ($\beta =-7$) type microswimmer launched from the initial condition $\theta _0 =187^{\circ }, h_0= 2$. Green circular markers represent damped oscillation of the trajectories in the upstream direction, whereas black square markers and red triangular markers represent escaping and periodic oscillation behaviours, respectively. The brown markers at $\mathcal {S} =0$ denote the mirror symmetric damped or periodic oscillation in the absence of a background shear. Beyond $\mathcal {S}=0.13$, for all slip lengths, trajectories become escaping in nature; therefore, in the case of pushers the critical shear line is marked in between $\mathcal {S}=0.13$ and $0.14$. (b) Three different values of shear rate (marked with square, triangular and circular enclosures in three different domains) are chosen at $l_s=0.1$ for plotting the trajectories in the $xz$ plane. The inset in (b) is a phase space trajectory, showing the attractor for $\mathcal {S}=0.01$, limit cycle for $\mathcal {S}=0.1$ and escaping nature for $\mathcal {S}=0.15$.

For $\mathcal {S}>0$, the figure shows that lower values of slip length result in undamped oscillations in the upstream direction (red markers) when the external flow intensity is in the range $\mathcal {S} \lesssim 0.01$. This trajectory feature gives rise to a limit cycle in the phase portrait. For this low range of $\mathcal {S}$, upon crossing a critical slip length $l_s \approx 0.1$, a transition from undamped to damped oscillations takes place in the upstream direction, thereby giving rise to rheotactic attractors. The phase space trajectory in the inset to figure 7(b) (green line) exemplifies the above-mentioned rheotactic attractor for $\mathcal {S}=0.01$.

Further enhancement of $\mathcal {S}$ beyond $\mathcal {S} = 0.01$ drives the microswimmer to escape (black) from the near-wall region when the no-slip condition applies at the wall. However, a simultaneous enhancement in $l_s$ and $\mathcal {S}$ causes a transition from escape to undamped oscillating trajectories (red). The line separating these two contrasting zones in the regime map thus becomes a piecewise function of both $l_s$ and $\mathcal {S}$ until a critical shear rate $\mathcal {S}_{cr} =0.13$ is encountered. Beyond this limit, all combinations of $l_s$ and $\mathcal {S}$ promote only the escaping trajectories of the microswimmer. An example of escaping trajectory (black) at $l_s=0.1$ and $\mathcal {S}=0.15$ has been shown in figure 7(b). It is also identified from the regime map that the critical slip length $l_{s,cr}$ required for the transition from undamped oscillations (red) to rheotactic states (green) increases as the external flow strengthens. The uneven phase boundaries in the regime map as a function of both slip length and shear strength indicate the rich interplay between the two competing mechanisms in deciding the swimming state of the microswimmer.

The variations of swimming states for different intensities of external flow can be explained by comparing the trajectories in figure 7(a). For a high shear rate $\mathcal {S}=0.15$, a strong CW rotation $\varOmega _y^{(ex)}$ overwhelms the tendency of the squirmer to rotate towards the wall under a CCW $\varOmega _y^{(sq)}$. Thus, an escaping trajectory results. However, at an intermediate strength of the external shear rate $(\mathcal {S}=0.1)$ the $\varOmega _y^{(sq)}$ and $\varOmega _y^{(ex)}$ become competitive in magnitudes, as highlighted in figure 8(b). Although the microswimmer initially rotates with a resultant CW $\varOmega _y$, the enhanced near-wall hydrodynamic interaction attracts the microswimmer towards the wall with $V_z<0$. A subsequent repulsive interaction imparts an upward motion to the microswimmer. Upon overcoming the short-range repulsive interaction, the microswimmer again moves downward. This motion behaviour becomes periodic in nature and undamped oscillations emerge.

Figure 8. Variation of rotational velocity and its components along the trajectory of a pusher $(\beta =-7)$ launched at $h_0=2$ and $\theta _0=187^{\circ }$ under external shear $(\mathcal {S}=0.1 )$ and slip length (a) $l_s=0.1$ and (b) $l_s=0.15$.

The mechanism behind rheotactic migration found exclusively in the presence of hydrodynamic slip can be explained by considering the effects of $l_s$ on the rotational velocity components $\varOmega _{y}^{(sq)}$ and $\varOmega _{y}^{(ex)}$ (see figure 8b) adjacent to the wall where repulsive interaction with the wall plays a prominent role. As the slip length crosses a threshold, the CCW rotation due to self-propulsion $(\varOmega _{y}^{(sq)})$ is perfectly balanced by an opposing shear-induced rotation $\varOmega _{y}^{(ex)}$ in the CW direction. On the other hand, the downward propulsive thrust is balanced by the repulsive forces, leading to a net zero velocity in the vertical direction. Consequently, a steady-state stable swimming results.

4. Conclusions and remarks

We have performed a theoretical investigation to elucidate the combined interplay of a background flow field and hydrodynamic slippage at a nearby plane wall in dictating the swimming characteristics of a model microswimmer. The self-propelling features of the microswimmer have been captured through the squirmer surface velocity condition at the sphere surface, while the Navier slip condition has been employed to track the surface wettability condition at the wall. We have presented an exact solution of the Stokes flow equation by exploiting the eigenfunctions in the bispherical coordinate system. Consequently, all the sub-problems associated with the flow dynamics have been summed up by exploiting the linearity of the flow problem. Furthermore, going beyond the widely employed image-singularity-based singularity analysis (Lopez & Lauga 2014) and the assumption of an asymptotically small slip length, we have provided a comprehensive solution of the flow problem for a diverse range of distances from the wall and arbitrary slip lengths.

The quasi-steady dynamics of a puller in three dimensions explicate that the enhancement of wall slip transforms an escaping or rheotactic trajectory (depending upon the orientation of the director) to a crashing one, as shown in figure 2. Subsequently, we have restricted the swimmer dynamics in the plane of shear and simplified the analysis based on phase portraits of a plane autonomous system. Our analysis aptly predicts near-wall trapping tendency of a puller (Li & Ardekani 2014; Daddi-Moussa-Ider et al. 2018; Poddar et al. 2020) in the absence of an external flow field. Correspondingly, a pair of oppositely directed mirror symmetric sliding states in the phase portrait has been reported. However, the application of a small strength of background flow alters one mirror symmetric sliding state to a rheotactic state, favoured by the orientation of the director (figure 3a,d). It is further reported that the rheotactic attractor appears along with a saddle point without affecting other fixed points in the phase portrait. Although the upstream rheotaxis of squirmers can be observed even in the presence of a no-slip wall (Ishimoto 2017; Walker et al. 2018; Uspal et al. 2015; Sharan et al. 2022), here, we bring out the immense characteristic modulation of rheotactic conditions and new rheotactic states triggered by finite slippage at the wall. Additionally, a minor enhancement of slip length $(l_s)$ from the no-slip condition draws the rheotactic attractor towards the wall (figure 3b,e) due to the action of an additional CW torque, and finally, the attractor annihilates (figure 3c, f) at the minimum separation distance $(\delta =0.01)$. The evolution of velocity components $\varOmega _y$ and $V_z$ (figure 4c–e) reveal that the slip-mediated torque enhances the negative magnitude of the corresponding velocities, resulting in the squirmer to collide before the CCW shear-assisted torque could reverse the motion direction.

We have reported the existence of a critical shear rate $\mathcal {S}_{cr}$, beyond which the slip-induced swimming characteristics bifurcate into two distinct regimes of emergence and annihilation (see regime maps in figure 6a,b). Initially, the strength of the background flow washes out the microswimmer due to insufficient torque triggered by low strengths of wall slip. However, the subsequent increase in slip length opens up a window of parametric combinations of $\mathcal {S}$ and $l_s$ by virtue of a strengthened CW torque acting on the microswimmer. As a consequence, the shear-dominated escaping trajectories transform to upstream rheotaxis for an intermediate zone of slip lengths, whereas crashing states reappear at higher values of $l_s$ due to wall contact.

In line with the earlier investigations (Lintuvuori et al. 2016; Poddar et al. 2020), we have observed both damped and undamped oscillations of pushers near a no-slip repulsive wall. The combined presence of slip and shear modulates the interaction between a repulsive wall and the hydrodynamic forces in such a fashion that many of the escaping and undamped oscillation states are converted to upstream damped oscillations or rheotaxis (figure 7a). Similar to the pullers, above a limiting strength of background shear, the microswimmer is forced to swim away from the wall (i.e. escape) due to a weaker contribution of the opposing slip-mediated torque (figure 7a). Thus, the pushers exhibit a set of rich thigmotaxis-assisted (Ishimoto 2017) motion governed by the strengths of the wall slip and background flow. The genesis of new rheotactic states for both pullers and pushers as an exclusive outcome of hydrodynamic slip is termed as slippery rheotaxis in the present work.

It is important to note that the current results of high wall slip in the micrometre scale should be used cautiously for a nano-engineered surface having inhomogeneous surface asperities of the order of the slip length (Choi & Kim 2006; Joseph et al. 2006). Moreover, a further detailed study on microswimmer movement near superhydrophobic confinements would account for the variable slippery condition along the wall (Ybert et al. 2007; Asmolov et al. 2013; Pimponi et al. 2014; Hu et al. 2015; Nizkaya et al. 2015). However, the present assumption of uniform surface slip remains applicable to the scenarios where the surface is atomically smooth, and ‘intrinsic slippage’ is exhibited due to monolayers of hydrophobic molecules (Huang et al. 2008; Sega et al. 2013; Gentili et al. 2014). Thus, the present study provides fundamental insights into the various complicated motion characteristics near hydrophobic confinements in the presence of a background flow encountered in diverse in vivo and in vitro conditions. Furthermore, the intriguing locomotion characteristics of spherical microswimmers regulated by well-defined combinations of the background flow strength $(\mathcal {S})$ and slip length $(l_s)$ may motivate the choice of parameters in future experiments on the navigation of microswimmers in complex and confined environments, and sorting of microswimmers based on the knowledge of sustained motion against the flow. In addition, it would be interesting to analyse the effects squirming modes of higher orders (Pak & Lauga 2014) and the role of non-spherical body shapes of micro-organisms (Shum et al. 2010; Ishimoto & Gaffney 2013) on the complicated flow physics considered in the present model.