# Rheotaxis of spherical active particles near a planar wall

## Abstract

For active particles the interplay between the self-generated hydrodynamic flow and an external shear flow, especially near bounding surfaces, can result in a rich behavior of the particles not easily foreseen from the consideration of the active and external driving mechanisms in isolation. For instance, under certain conditions, the particles exhibit “rheotaxis,” i.e., they align their direction of motion with the plane of shear spanned by the direction of the flow and the normal of the bounding surface and move with or against the flow. To date, studies of rheotaxis have focused on elongated particles (e.g., spermatozoa), for which rheotaxis can be understood intuitively in terms of a “weather vane” mechanism. Here we investigate the possibility that spherical active particles, for which the “weather vane” mechanism is excluded due to the symmetry of the shape, may nevertheless exhibit rheotaxis. Combining analytical and numerical calculations, we show that, for a broad class of spherical active particles, rheotactic behavior may emerge via a mechanism which involves “self-trapping” near a hard wall owing to the active propulsion of the particles, combined with their rotation, alignment, and “locking” of the direction of motion into the shear plane. In this state, the particles move solely up- or downstream at a steady height and orientation.

###### pacs:

47.63.Gd, 47.63.mf, 64.75.Xc, 82.70Dd, 47.57.-s## I Introduction

Flows are ubiquitous in microconfined fluid systems. Examples include the flow of blood
through capillaries, the flow of groundwater through pore spaces in soil, the flow of oil through
porous rock, and the flow of analyte through the channels of a lab-on-a-chip device. Flow
presents both challenges and opportunities for natural and engineered active particles,
especially near surfaces. For instance, fluid shear can trap bacteria near surfaces, inhibiting
chemotaxis, but promoting surface attachment and biofilm formation Rusconi et al. (2014).
Active particles can harness shear flow for sensing and navigation. A motile
spermatozoan in shear flow near a surface orients its body axis along the flow direction
– a phenomenon known as rheotaxis^{1}

Experimental studies of the rheotaxis of sperm near surfaces have a long history. Bretherton and Rothschild investigated live and dead sperm under flow in various device geometries Bretherton and Rothschild (1961). They observed that live sperm always exhibited motion upstream, i.e., upstream rheotaxis. Subsequently, Rothschild observed that motile sperm tend to accumulate near surfaces. He attributed this to an attraction involving the sperm heads Rothschild (1963). If heads are attracted to surfaces, then reorientation in flow can be qualitatively understood as due to the greater drag on the tail, which points into the bulk, than on the head, which is near the surface. With the head acting as a pivot point, the flow drags the tail downstream, like a weather vane Winet et al. (1984); Miki and Clapham (2013); Kantsler et al. (2014). Accumulation of sperm and other micro-organisms at surfaces is evidently driven by hydrodynamic interactions Fauci and McDonald (1995); Lauga et al. (2006); Berke et al. (2008); Smith et al. (2009); Elgeti et al. (2010); Spangolie and Lauga (2012) and steric repulsion Li and Tang (2009); Elgeti et al. (2010); Kantsler et al. (2012, 2014), although the relative significance of these contributions is a subject of debate Li and Tang (2009); Kantsler et al. (2012, 2014).

Recently, the experimental efforts have been extended to the study of bacterial rheotaxis. Like sperm, bacteria in quiescent fluid are observed to accumulate near surfaces Lauga et al. (2006); Berke et al. (2008); Li and Tang (2009). Hill et al. studied E. coli under flow in a microfluidic channel Hill et al. (2007). They observed that the bacteria in the middle of the channel tend to orient, relative to the flow direction, with a certain acute angle. They drift to the sides of the channel, migrating across fluid streamlines. Upon reaching the side walls, the bacteria orient against the flow and swim upstream. As with sperm, these observations could be rationalized with a “weather vane” mechanism. More recently, Kaya and Koser observed that, within an intermediate band of shear rates, E. coli can directly align against the flow and swim upstream, without having to migrate to side walls Kaya and Koser (2012).

Furthermore, we note that for micro-organisms in the bulk (far from bounding surfaces), the interplay of shear, swimmer geometry, and swimming strategy can drive rotation of the body axis towards the vorticity direction (which is perpendicular to the shear plane), allowing the organism to migrate across streamlines. In shear flow, bacteria driven by coiled flagella, such as B. subtillis, experience a lift force in the vorticity direction, owing to their chirality. Since this chirality is localized at the flagella, a torque is exerted on the cell body, driving rotation towards the vorticity direction Marcos et al. (2012). Micro-algae also rotate towards the vorticity direction. The mechanism for this orientation is not yet known, but it may be rooted in the fact that their two flagella beat at different rates, producing unequal propulsive forces and, consequently, a torque on the cell body Chengala et al. (2013). This attraction of the orientation vector to the vorticity direction has also been labeled rheotaxis Marcos et al. (2012), but throughout the present study we restrict this notion to denote attraction of the orientation vector to the shear plane.

Theoretical studies have sought to reproduce and isolate generic aspects of swimming in flow, including rheotaxis. Zöttl and Stark considered a model spherical microswimmer driven by Poiseuille flow in narrow tubes and slits Zöttl and Stark (2012). They found that hydrodynamic interactions with the confining walls can stabilize upstream swimming. So-called “pullers” migrated to the centerline and swam upstream, while “pushers” were attracted to trajectories in which they swim upstream while oscillating between the walls. In numerical simulations of elongated run-and-tumble particles (posing as model bacteria) driven by Poiseuille flow through slit-like channels, the particles accumulated near walls due to steric repulsion between them and the wall and, due to the subsequent alignment of their bodies by the flow, swam upstream Nash et al. (2010); Costanzo et al. (2012). Chilukuri et al. modeled swimmers as Brownian “dumb-bells” (i.e., two beads connected by a spring) Chilukuri et al. (2014). Likewise, when driven by Poiseuille flow through a wide slit, a fraction of the particles accumulated near surfaces and swam upstream. This accumulation was enhanced upon including hydrodynamic interactions between the dumb-bells and the walls, indicating that, for this system, both steric repulsion and hydrodynamic interactions promote rheotaxis. Additionally, for a Brownian self-propelled particle (i.e., one exposed to stochastic forces), ten Hagen et al. found that an unbounded shear flow can modify the scaling with time of the mean squared displacement of the particle ten Hagen et al. (2011).

In the last decade, significant progress has been made in developing synthetic self-propelled particles. These particles are envisioned as key components in future lab-on-a-chip, drug delivery, and pollution remediation systems Ebbens and Howse (2010); Patra et al. (2013). In these and other applications, active particles will undoubtedly encounter shear flows near surfaces and will need to autonomously sense and respond to flow. Many generic aspects of artificial microswimmers in flow are likely captured by the theoretical and numerical studies discussed above. However, achieving rheotaxis in artificial systems, for which active stresses can be generated through, e.g., local changes in the chemistry of the environment, requires an in-depth analysis of the interplay between the external flow, the presence of bounding surfaces, and the mechanism of propulsion (i.e., activity). For example, for a large class of active particles, a reaction involving a “fuel” molecule in the surrounding solution is catalyzed by a region of the particle surface, leading to propulsion either via the formation and expulsion of bubbles Sanchez et al. (2011), or via “self-phoresis,” i.e., the generation of a tangential surface flow powered by the gradient of the product and/or reactant molecules Anderson (1989); Golestanian et al. (2005). For the latter case, we have recently shown that the complex behavior emerging when the particle motion occurs near confining boundaries can be understood only if the coupling loop between the distributions of chemical species and hydrodynamics is fully accounted for Uspal et al. (2015).

There have been comparatively few studies of artificial microswimmers in flow, and the issue of rheotactic behavior was rarely raised. Sanchez et al. have shown that bubble-powered, rod-shaped catalytic “microjets” can move upstream in a channel when guided by an external magnetic field Sanchez et al. (2011). Frankel and Khair studied theoretically a model of self-phoretic Janus particles and have shown that drift across streamlines in unbounded shear flow may occur at finite Péclet number Frankel and Khair (2014). Tao and Kapral have used numerical simulations in order to investigate the motion of a self-phoretic dimeric nanomotor moving upstream against a Poiseuille flow in a channel and to examine the effect of confined flow on the propulsion velocity of the motor Tao and Kapral (2010). However, they used a potential which confined the motor to the channel centerline, and therefore did not probe the possibility of rheotaxis. Recently, Palacci et al. employed microfluidic chips and self-propelled dimers made out of polymer spheres and hematite cubes partially embedded in the polymer matrix in order to provide the first experimental demonstration of rheotaxis, via the “weather vane” mechanism, for an artificial swimmer Palacci et al. (2015). The heavy dimers sediment to the bottom surface and, once being exposed to blue light and to hydrogen peroxide “fuel”, the hematite becomes catalytically active and promotes the decomposition of the hydrogen peroxide. This leads to an attraction of the hematite end to the surface and thus, while the whole dimer tends to move due to the concentration gradients, the hematite end serves as an “anchor” for the dimer to orient with the ambient, externally controlled flow within the microfluidic chip.

In the present study we address the issue of whether artificial spherical active particles can exhibit rheotaxis when moving in the vicinity of a planar surface. Without a distinction between major and minor axes, the particle cannot align via the intuitive “weather vane” effect. Thus it is not a priori clear if rheotaxis is possible for such systems; for rheotatic behavior to emerge nonetheless a different physical mechanism must be at work. Considering the overwhelming importance of spherical particles for experiments and applications of artificial swimmers (a spherical particle is naturally simpler and easier to fabricate than an elongated body), this question out of basic research turns out to be one of significant practical interest, too. Furthermore, the answer to this question may also shed light on the motility in confined flows of micro-organisms with spherical or nearly spherical shape.

We focus our theoretical and numerical study on the case of spherical particles with an axisymmetric propulsion mechanism which can be described via an effective surface slip velocity. This framework captures a number of important classes of microswimmers, including self-phoretic catalytically active Janus particles, as well as the classical “squirmer” model, introduced by Lighthill Lighthill (1952) and subsequently refined by Blake Blake (1971), which is known to account very well for the essential features of the self-propulsion of ciliated micro-organisms. We show that if the swimming activity of the particle can lead to a “self-trapping” near a bounding surface rheotactic behavior may emerge. We establish the conditions under which this “self-trapping” is accompanied by a mechanism for resisting rotation by the flow and maintaining a steady orientation. In contrast with the “weather vane” mechanism, rotation of the orientation vector (i.e., the symmetry axis of the particle) into the shear plane is now driven not by the external flow, but rather by the near-surface swimming activity of the particle. Building on these results, as well as on our recent study of the self-propulsion of a catalytic Janus particle near a surface in a quiescent fluid Uspal et al. (2015), we show how the surface chemistry of a catalytically active Janus particle can be designed to achieve rheotaxis. The wide applicability of our results is further emphasized by demonstrating that rheotactic behavior also occurs for squirmers which in quiescent fluids are attracted by surfaces Ishimoto and Gaffney (2013); Li and Ardekani (2014).

## Ii Theory

### ii.1 Problem formulation

We consider a spherical, axisymmetric, active particle suspended in a linear shear flow in the vicinity of a planar wall (Fig. 1(a)). We adopt a coordinate frame (the “lab frame”) in which the wall is stationary. The wall is located at and has a normal vector in the direction. The particle has radius and its center is located at ; denotes the height of the particle center above the wall. We describe the particle orientation in terms of the “director” , , which is the direction of active motion if the particle would be suspended in an unconfined, quiescent fluid. Due to axial symmetry, is parallel to the axis of symmetry of the particle. (For instance, for a catalytic Janus particle which moves “away” from its catalytic cap, the direction of is given by the vector from the pole of the catalytic cap to the pole of the inert region.) The height and orientation completely specify the instantaneous particle configuration. We shall find it convenient to alternatively describe the particle orientation in terms of spherical coordinates and , where is the angle between and , and is the angle between and the projection of onto the plane (see Fig. 1(a)). We seek to calculate the translational velocity and angular velocity of the particle as a function of the configuration (, ).

The active particles we consider here are such that the self-propulsion mechanism can be accounted for via an effective surface slip velocity . The slip velocity can encode, for example, the surface flow generated by a “squirming” mechanical deformation of the particle surface (e.g., the time averaged motion of the cilia of a micro-organism), or the surface flow generated by the interactions between reactant and product molecules with the particle (as in the case of a catalytically active Janus colloid). Here we take the slip velocity as a known function of the position along the particle surface and of the configuration (, ). We postpone for Section II.E the discussion of how to obtain this slip velocity for a given swimmer model.

The velocity in the suspending fluid is , where is the external background flow, is the disturbance flow created by the particle, and denotes the spatially and temporally constant shear rate. The fluid velocity is governed by the Stokes equation , where is the pressure field and is the viscosity of the solution, as well as the incompressibility condition . On the planar wall, the velocity satisfies the no-slip boundary condition . On the surface of the particle, accounting for the effective slip implies that the fluid velocity satisfies . Finally, we specify that the active particle is force and torque free, thus closing the system of equations for the six unknowns and .

### ii.2 Solution by linear superposition

We exploit the linearity of the Stokes equation in order to split the complete problem into two subproblems. In subproblem I we consider the motion of an inert (non-active) sphere with configuration driven by an external shear flow. We note that although, for consistency of notations, is included as a configuration variable, in the case considered in this subproblem the velocity of the particle actually does not depend on (but, obviously, will evolve in time due to rotation of the particle.) In subproblem II, a spherical active particle with instantaneous configuration moves through a quiescent fluid. The particle velocity for the complete problem is then obtained by linear superposition: and , where the velocities and are due to flow (subproblem I), and the velocities and are due to the activity-induced motion (subproblem II).

Subproblem I was solved analytically by Goldman et al. using bispherical coordinates Goldman et al. (1967). The translational and rotational velocities of the particle are given by

(1) |

Note that the external flow does not cause motion of the particle in the direction. The functions and encode the hydrodynamic interaction with the no-slip wall, which retards both the translation and the rotation of the particle. These functions increase monotonically with and have the properties , , , and . They are tabulated for selected values of in Ref. Goldman et al. (1967).

We now turn to subproblem II. Due to the symmetry of the system and the absence of thermal fluctuations, the motion of the particle is confined to the plane spanned by the orientation vector and the wall normal . Therefore it is convenient to introduce a “primed” reference frame to solve this subproblem (Fig. 1(b)). In the primed frame, the wall is stationary and located at . The vector is normal to the wall, the vector lies in the plane containing the particle symmetry axis and (which means that the projection of onto the plane points into the direction of ), and is orthogonal to and . In general, this subproblem must be solved numerically. We discuss the numerical solution in Section III. Here we assume the solution to be known and we use the most general form compatible with the symmetry constraints discussed above, i.e., a translational velocity and an angular velocity , where , , and are functions of and . Transforming back to the original reference frame, we obtain

(2) |

### ii.3 Particle motion and fixed points for the dynamics

In the absence of thermal fluctuations (which we assume to be negligible), the time evolution of the particle configuration ) is determined, in the overdamped limit, by the translational and rotational velocities derived above. By noting that the motion of the director obeys and that (see Fig. 1(a))

we arrive at the following dynamical equations (using , , ):

(3) | |||||

(4) | |||||

(6) |

where the dot over a quantity denotes its time derivative. Since , the dynamics of the director must obey the constraint ; this is obviously satisfied by Eqs. (3)-(II.3). Since, accordingly, the three components of are not independent, the dynamical system is de facto three-dimensional.

Defining the generalized velocity vector , we search for fixed points of the dynamics, which are determined by . Such configurations correspond to the particle translating along the wall with a fixed height and orientation, i.e., steady states with properties which are compatible with (although not all of them are necessary for) rheotaxis. In order to obtain , we start by considering the three possibilities for to be satisfied. For completeness, we discuss here all three cases and, for each of them, derive the additional conditions following from , , and . In anticipation, we note that for the types of swimmer considered here, it turns out that only the third case, which corresponds to , is compatible with rheotaxis. Briefly, in the first two cases, the combined effects of shear and activity produce a fixed point with , i.e., . Most types of swimmer do not satisfy the restrictive conditions required to produce such a fixed point. Consequently, the rest of the study in this paper will focus on the third case.

(i) The case . This is a particle
configuration for which in the absence of external flow the particle only translates. In
order to also satisfy , , and , it is
then required that , , and . The first two of the latter relations imply that the director is along
the -axis, which coincides with the vorticity axis of the external shear
flow^{2}

(ii) The case . In this case the director lies within a plane parallel to the wall (see Fig. 1(a)). The equation is then automatically satisfied (Eq. (4)). The condition implies that the curve defined by passes through (i.e., ) at a certain height . The condition then implies . The corresponding physical picture is as follows. As the director is within a plane parallel to the wall, the contributions to from both shear and activity are entirely along the direction: for Eqs. (3) - (II.3) imply . The magnitude of the contribution from the shear flow ( in Eq. (II.3)) depends on the angle . For instance, if the director is oriented along the -axis (, so that ), the contribution of shear to is zero (see Eq. (II.3)). At (see above, i.e., for a certain angle ), due to the definition of the contributions from shear and activity to sum to zero. For this fixed point to occur, due to it is required that . This requirement, together with that of the vanishing of at a certain height for , places weaker demands on subproblem II than the requirements in case (i). Nevertheless, as in case (i), it turns out that none of the types of active particles which will be considered in Sec. III satisfies the fixed point requirement on , and thus none of them can fit case (ii) either.

(iii) The case . In this case the director lies within the shear plane (i.e., the plane in Fig. 1(a)). If a fixed point exists, in that state the particle translates either only upstream or only downstream because once having reached the fixed point of the dynamics the orientation cannot vary in time. If the fixed point is stable (this issue will be addressed analytically in the next subsection and numerically in Sec. III), this state corresponds to rheotaxis. Since in this case and thus , Eqs. (4) and (II.3) are identical (as expected, because only two of the Eqs. (3) - (II.3) are independent). Therefore, the steady height and orientation are obtained as the solutions of

(7) |

The first equation simply expresses that at the steady state the sum of the contributions to the angular velocities from shear and activity must add up to zero. The term appears because shear tends to rotate the director away from the wall when (i.e., ) and towards the wall when (i.e., ) (see Fig. 1(a) and recall that in the present case lies in the plane), whereas, owing to the symmetry of subproblem II, there is no dependence on for the direction of rotation by the swimming activity. In contrast to the cases (i) and (ii), it turns out that the conditions in Eqs. (II.3) can be satisfied by all three classes of spherical active particles considered here (see Sec. III). Therefore we proceed further with the analysis of the case .

### ii.4 Linear stability analysis of the fixed point

In order to establish the conditions under which the fixed point (case (iii) above) is stable, so that the corresponding steady state exhibits rheotaxis, we proceed by carrying out a linear stability analysis. Since for the director lies within the shear plane and, due to the symmetry of the problem, cannot rotate out of this plane, in the three-dimensional phase space defined by , , and a trajectory with initial condition in the plane will be confined to that plane for all times . Upon linearization of the equations of motion near the fixed point, a small perturbation decouples from the other variables:

(8) |

After using Eq. (II.3) for substituting

, we obtain

(9) |

Therefore decays or grows exponentially,

(10) |

depending on the sign of the quantity multiplying in the argument of the exponential.

Since and , the stability of the fixed point is determined by the sign of the ratio . For the particle to swim oriented upstream (i.e., to exhibit “positive” rheotaxis), it must align against the flow, so that ; stability then requires that , i.e., the director points towards the wall. Additionally, if the swimming velocity is large enough to overcome the shear flow, i.e., for , where (see in Eq. (II.2) and in Eq. (II.2) with ), net upstream swimming emerges. The other possibility, that an attractor exists for which and are both positive, is in general unlikely, because the particle would be pointing away from the wall (see Fig. 1(a)). Furthermore, it would correspond to the less interesting case of negative rheotaxis, in which the particle motion is in the flow direction. Therefore, it will not be discussed further.

The analysis of the linear stability shows that we need to examine particle dynamics only in the symmetry plane in order to determine whether rheotaxis occurs. In the following, we shall examine the two-dimensional dynamics in this plane in detail. Additionally, we shall perform fully three-dimensional numerical calculations of trajectories with initial conditions in order to explore the limits of validity of the linear analysis, as well as to understand the details of the evolution towards alignment in the shear plane. We introduce a new angle , shown in Fig. 2(a), which will turn out to be convenient for describing the director orientation if . The requirement that the attractor occurs with and translates into .

Before proceeding with the analysis, we discuss the physical picture behind the stability criteria in Eq. (10). While we have provided mathematical conditions for rheotaxis, the corresponding physical interpretation can provide further insight. Equation (3) shows that particle activity has an effect on , but shear does not. The contribution of shear to the angular velocity of the particle is rigid rotation around the axis. In the absence of activity, shear rotates the tip of the director in a coaxial circle around the axis, keeping constant. Now we consider activity. Due to the mirror symmetry with respect to the plane containing , the contribution of activity can only drive rotation of the director with arbitrary orientation towards or away from the planar wall. On the spherical surface defined by , this corresponds to rotation of the director along a line of longitude, where and define the poles [Fig. 2(c)]. The question, then, is whether, for a small perturbation , activity will tend to drive the director towards the equator of or towards the nearest pole. From Fig. 2(c) we note that the effect of the small perturbation is to shift the tip of the director to a neighboring line of longitude. Therefore, if activity tends to rotate the director towards the equator, will increase (see Fig. 2(c)) because the distance between neighboring lines of longitude increases as they approach the equator from the pole. The direction of this activity induced rotation (i.e., towards the pole or towards the equator) is determined by the value of . For instance, consider a positive rheotactic state, for which . In order to balance rotation by shear, activity must tend to rotate the director towards the pole , as shown in Fig. 2(b) and (c). Hence, this dynamical fixed point is stable against small perturbations in .

### ii.5 Calculation of the slip velocity

Finally, we briefly outline the calculation (if necessary) of the slip velocity (where denotes the unit vector in the primed coordinate frame corresponding to the polar angle measured from the director, i.e., the symmetry axis of the particle) along the surface of the spherical particle which, for the specific swimmer models considered here, is an input to the problem outlined in Subsec. II.A.

For a “squirmer”, is specified by fiat and does not depend on the configuration of the particle. Following Li and Ardekani Li and Ardekani (2014), we consider a squirmer for which the slip velocity is given by

(11) |

The model is characterized by the squirming mode amplitudes and . In an unconfined quiescent fluid (i.e., free space) such a squirmer moves with velocity .

A self-phoretic Janus particle generates solute molecules over the surface of a catalytic cap. We parameterize the extent of the cap by , where the angle is defined in Fig. 2(a). It is assumed that the interaction between the active particle and the (non-uniformly distributed) solute molecules has a range much smaller than and . Thus it can be accounted for as driving a tangential flow in a thin boundary layer of thickness surrounding the particle. This can be modeled as an effective slip velocity , where is the solute number density field Anderson (1989); Golestanian et al. (2005). The operator is the projection of the gradient onto the active particle surface, with denoting the surface normal oriented into the fluid. The quantity is the so-called “surface mobility,” determined by the details of the interaction between the active particle and the solute molecules. For instance, the sign of expresses whether the active particle is attracted to or repelled from the solute molecules. Throughout the present study, we assume repulsion. Since the solute number density is affected by the presence of the wall, in this model the slip velocity depends on the configuration .

If the effect of advection on can be neglected, one can determine (and with this ) prior to the consideration of the hydrodynamic problem posed in Subsec. II.A (see, e.g., Refs. Anderson (1989); Popescu et al. (2009)). (Within this approach it is disregarded how the flow field couples back to .) In this case, the solute number density is quasi-static, and obeys the Laplace equation , subject to the boundary conditions for the normal derivatives on the wall, over the inert surface of the colloid, and on the catalytic cap. The quantity is the diffusion coefficient of the solute molecules, and is the flux of solute per unit area from the catalytically active cap. Concerning the first and second boundary conditions, the wall and the inert surface of the particle are taken to be completely impenetrable to solute molecules. In the third boundary condition, it is assumed that the rate of the reaction is independent of the local concentration of solute molecules.

It is valid to neglect advective effects if the Péclet numbers and are small, where is the particle velocity in free space. (The Péclet number is the ratio between the time scale for solute to diffuse a distance and the timescale for the particle to move a distance . Similarly, the Péclet number characterizes the advection of the solute molecules by the external flow, for which the velocity scale is , relative to the diffusion of the solute molecules Frankel and Khair (2014).) Additionally, we assume that the Reynolds number (where is the mass density of the suspending fluid and is its dynamic viscosity) is small, permitting the use of the Stokes equations to account for the hydrodynamics. For example, a self-propelled Janus colloid which catalyzes the decomposition of hydrogen peroxide into water and oxygen typically has a diameter of and moves with a velocity through the aqueous solution ( and ). In this case, . Furthermore, estimating the diffusion coefficient of oxygen in water as , we obtain Popescu et al. (2010). A representative shear rate in microfluidic devices is (see, for instance, Refs. Kantsler et al. (2014) and Chengala et al. (2013)), which renders . Thus for typical active particles the above assumptions of small , , and numbers are well justified.

Before proceeding to an analysis of the numerical results, we remark on the physical realizations which would be captured by our model. In many experiments, the catalytic reactions involve more than one product as well as possibly a number of reactants. If the system remains in a reaction-rate limited regime (i.e., the reactants are in abundance and transported sufficiently fast so that there is no noticeable depletion near the catalyst), accounting for more than one product means to modify the expression for the slip velocity in a linear manner: the gradient of each reaction product “i” multiplied by its corresponding surface mobility is simply added in order to obtain the phoretic slip around the particle. In this case, the results presented here can be easily extended by replacing the effective number density by the sum over the densities of all products with an effective such that . On the other hand, if the catalytic reaction is diffusion limited and involves at least two reactants, the source boundary condition at the catalytic-active region involves products of the densities of the reactants. In this case the equations describing the system become nonlinear, and the present model is neither applicable, nor does it lend itself to an obvious extension.

## Iii Numerical Results and Discussion

When considering a Janus particle, we use the boundary element method (BEM) in order to solve the diffusion equation for the solute number density field . Additionally, for both propulsion mechanisms (the squirmer and the catalytically active particle), we use the BEM in order to solve the hydrodynamic subproblems I and II. A detailed introduction to the BEM is provided by Pozrikidis Pozrikidis (2002). We adopted his freely available BEMLIB code and modified it for the present study. For subproblem I, we calculate and for an (inactive) sphere in a shear flow at various heights . We find that the numerically obtained values are in good agreement with the analytically obtained values given by Goldman et al Goldman et al. (1967). In a previous study, we provided a detailed description of how to apply the BEM in order to solve the diffusion equation (for a catalytic Janus particle) and subproblem II Uspal et al. (2015). Within this approach we have obtained and on a grid of and .

In order to determine full particle trajectories for an initial condition , we interpolate , , , and and integrate Eqs. (3) - (6) numerically. In the following, we restrict our numerical calculations to the region . Above , we consider the particle to have “escaped” the surface. Below , various of our physical approximations (e.g., that the effects of the solute distribution around the particle can be accounted for by a phoretic slip calculated within a thin layer, including that part of the surface of the particle which is in close proximity of the wall) are expected to break down. For a catalytic Janus particle, provides a characteristic velocity, and a characteristic timescale. For instance, it was shown analytically that a catalytic Janus particle with uniform surface mobility and half coverage by catalyst has a velocity in free space Popescu et al. (2010). When we take the catalytic cap and the inert region to have unequal surface mobilities, we use . Similarly, for squirmers the amplitude of the first squirming mode provides a characteristic velocity scale, and the time a characteristic timescale. The ratio between the squirming mode amplitudes defines a parameter . In each of the cases which we shall study in this section, the corresponding velocity and time scales defined above and the particle radius will be employed to render the dynamical equations, Eqs. (3) - (6), dimensionless.

First, we shall apply the theoretical results derived in the previous section to show that the surface chemistry of a catalytically active Janus particle can be tailored such that it leads to the occurrence of positive (upstream) rheotaxis. We shall provide two rather distinct examples of such a design of the surface properties, each exploiting a particular pathway to produce the stabilizing wall-induced rotation component discussed in Subsec. II.4. These two examples are depicted schematically in Fig. 3. In order to illustrate the generality of our theoretical results, we shall show that rheotaxis can occur for certain spherical squirmers.

Our approach to design is guided by the idea outlined in Fig. 2(b). For positive rheotaxis, the particle director must point upstream and towards the wall. The external flow contributes a clockwise component, shown by the cyan arrow, to the angular velocity of the particle. For the particle to maintain a steady orientation, the near-wall swimming activity of the particle must contribute a counterclockwise component , shown by the magenta arrow, to . Since the axisymmetric particle does not rotate in free space, the counterclockwise component must be due to the effect of the wall on the fluid velocity field and the solute number density field. Additionally, the particle must be attracted to a steady height through its near-wall swimming activity. In other words, the tendency of the particle to swim away from its cap must be counteracted by the effect of the wall on the fluid velocity and solute density field. In the following two subsections, we shall introduce two particle surface chemistries which fulfill these criteria. By tailoring the surface chemistry, we can turn on or off, and rationally control, various physical mechanisms which contribute to the particle motion.

We note that the following three subsections (corresponding to the three particle designs shown in Fig. 3) have a repetitive structure of arguments and presentation, so that the reader can choose to read the first subsection, skip ahead to the Conclusions, and return to read Subsections B and C at leisure.

### iii.1 Catalytic Janus particles with high coverage and uniform surface mobility

Previously Uspal et al. (2015) we studied the dynamics of a model catalytically active
Janus particle suspended in a quiescent fluid and near a wall. We found that, for
a uniform surface mobility , in the course of time a particle with a very high
coverage by catalyst can be stably attracted to a “hovering” state in which it remains
motionless at a height and angle (see Fig.
2(b)), i.e., with its catalytic cap oriented
vertically and away from the wall Uspal et al. (2015). In this state, depicted in the left panel of Fig. 4, the tendency of the particle to translate away from its cap (i.e., trying to avoid high solute concentrations) is
balanced by the accumulation of a solute “cushion” near the impenetrable wall
(which is due to the confinement of the solute between the particle and the wall).
The stability of this hovering state against perturbations in
the direction can be understood easily: if the particle moves closer to
the wall,
solute accumulation is enhanced, and repulsion from the wall is stronger; if the particle
moves away from the wall, repulsion from the wall is weakened, but the particle still
translates away from its cap. Less obviously, hydrodynamic interactions with the wall
stabilize the particle against perturbations in the angle . Hydrodynamically, the
particle is characterized as a “puller” (see the flow lines in Fig. 3(b) in Ref.
Uspal et al. (2015)). Pullers are known to orient themselves perpendicular to planar surfaces
via hydrodynamic interactions Spangolie and Lauga (2012).^{3}

The phase plane for a hoverer with and uniform surface mobility in a quiescent fluid (i.e., no external flow) is shown in Fig. 5(a). The phase plane indicates the evolution of the particle configuration for any initial condition . In the absence of flow, the system is symmetric for rotations around the axis. This rotational symmetry is evinced by two mirror symmetries in the phase plane: symmetry of the region for reflection across , and symmetry of the region for reflection across . The “hovering” attractor is shown as a filled magenta circle at and . As discussed in Ref. Uspal et al. (2015), many trajectories in the basin of attraction are drawn to a “slow manifold” indicative of two timescale dynamics: rotation is much slower than motion in direction. For initial conditions , the particle escapes the surface.

Now we consider the same “hoverer” in shear flow. The flow will rotate the particle clockwise, decreasing from . On the other hand, the hydrodynamic interaction which stabilized “hovering” will tend to rotate the particle back towards . As depicted in the center and right panels of Fig. 4, these two contributions to rotation can balance for sufficiently large angular displacement from the cap-down orientation. Hence the right panel of Fig. 4 provides a potential realization of the mechanism illustrated in Fig. 2. Therefore we anticipate that, at least for sufficiently slow external flows, there will a stable angle . Moreover, we anticipate that the “cushion” effect will be preserved to produce a steady height . The numerical simulations for, e.g., lead to a phase plane as shown in Fig. 5(b) and confirm these expectations: the attractor (filled green circle) is preserved and it migrates into the region , being now located at and . Trajectories from a large section of phase space are drawn to this attractor. (Note that the phase plane is periodic in the direction.) Additionally, we indicate a saddle point (orange cross). The saddle point and attractor are rather close, as we have chosen a shear rate close to the (numerically estimated) upper critical value . At the critical value, the saddle point and the attractor collide and annihilate each other. Above the critical shear rate, there is no stable rheotaxis, and, for this choice of the surface chemistry, all trajectories escape from the surface. We have chosen near the upper critical shear rate since, as will be discussed in the Conclusions, “strong” external flows have the greatest experimental relevance and accessibility. Additionally, the relaxation time scale in Eq. (10) decreases as the shear rate is increased. Therefore, we expect the approach towards the rheotactic state to occur most rapidly for being close to the upper critical shear rate. Although we leave an exhaustive parametric study to future research, we note that results of additional numerical calculations at lower values of , omitted here, agree with this expectation.

Since the attractor is in the region , it should be stable against perturbations of the director out of the shear plane, as discussed in Sec. II. In order to probe the stability of the attractor and its basin of attraction, we consider an ensemble of trajectories launched from various initial director angles and (see Fig. 1(a)), initial heights , , and , and the initial lateral position . For these three initial heights Figs. 6(a), (b), and (c) show three-dimensional trajectories of the coordinates of the centers of the particles. Each trajectory has been obtained from an initial position and orientation by numerically integrating Eqs. (3)-(6). Green trajectories are “rheotactic.” Particles following these trajectories are attracted to the steady height and angle and move upstream. Particles following blue trajectories escape from the surface (i.e., the trajectories cross from below.) Phase maps, indicating how the particle behavior depends on the initial height and orientation, are shown in Figs. 6(d), (e), and (f), which correspond to , , and , respectively. Clearly, rheotaxis is achieved for a large basin of initial conditions.

An example of a rheotactic trajectory is shown in Fig. 7. Starting from the initial orientation , , and , the particle has nearly attained the steady height and orientation after moving only a few particle diameters. The particle is attracted to a configuration in which it is tilted slightly away from the “hovering” state by the shear flow. As the (blue) cap is oriented slightly downstream (i.e., ), the particle moves upstream. In Fig. 8, the particle starts near the wall, but pointing away from it. Due to this initial orientation, the particle moves a few diameters away from the wall, where the hydrodynamic and chemical influence of the wall is very weak. However, the flow rotates the particle to swim back towards the wall. Upon returning to the vicinity of the wall, the particle rotates into the “tilted hoverer” configuration.

For selected rheotactic trajectories, in Fig. 9 we plot the time evolution of the particle height and of the components of the director. As expected, the component of the director perpendicular to the shear plane decays to zero. The height and the two in-plane components and attain asymptotically non-zero values. In Fig. 9(e), we plot the time evolution of for all rheotactic trajectories studied for Fig. 6. The decay of is clearly exponential, and the timescale for decay closely agrees with the prediction of Eq. (10), which is plotted as the dotted line in Fig. 9(e). Note that in Figs. 7, 8, and 9 time is given in dimensionless units as . In a previous section, has been defined as . As will be discussed in the Conclusions, we estimate to be for catalytic Janus particles as typically used in experiments.

### iii.2 Catalytic Janus particle: half coverage and inhomogeneous surface mobilities

In order to demonstrate the general character of our theoretical findings, we seek to reach the rheotactic state of Fig. 2(b) via the alternative pathway of a different surface chemistry designed to realize a distinct physical mechanism. Specifically, we consider a particle which is half covered by catalyst (), but the surface mobilities of the inert region and of the catalytic cap are taken to differ: . In part, this choice is motivated by the fact that by default experimental studies use particles with half coverage, because these can readily be prepared by vapor deposition. Moreover, since the catalytic and inert surface regions consist of different materials, they are likely to give rise to different surface mobilities, too.

In our previous study of systems without flow, we have isolated two distinct wall-induced contributions to the rotation of a Janus particle Uspal et al. (2015). One contribution is due to the hydrodynamic interaction of the particle with the wall. Disturbance flows created by the motion of the particle are reflected by the wall, coupling back to the particle. We found that, for half coverage (), hydrodynamic interactions always tend to rotate the catalytic cap of the particle towards the wall (and thus the director away from the wall). Therefore, this contribution to particle rotation cannot oppose the rotation by the shear flow (cyan arrow in Fig. 2(b)) for , i.e., it cannot provide the magenta arrow in Fig. 2(b). However, we also found that if , wall-induced chemical gradients can contribute to particle rotation. If , repulsion of the solute (i.e., the reaction product) from the inert region is weaker than repulsion from the catalytic cap. Accordingly, wall-induced chemical gradients (i.e., a higher concentration of solutes on the side of the particle surface closer to the wall) tend to rotate the catalytically active cap away from the wall. (Note that chemical gradients will not drive rotation of the particle in the absence of the wall because the Janus particle is axially symmetric.) Therefore, taking generates a contribution to the rotational velocity of the particle which corresponds to the magenta arrow in Fig. 2(b).

Therefore, we consider a particle with half coverage () and
. A close-up of the phase plane for the dynamics is shown in Fig.
10(a) for the case that the director lies in the
plane and that there is no external flow. There is an attractor very close to the wall at
and . At this point, the hydrodynamic
interaction with the wall, which tends to rotate the cap of the particle towards the wall,
is balanced by the effect generated by wall-induced chemical gradients, which
tends to rotate the cap away from the wall. The particle moves along the wall with a
steady height and steady orientation, which we called a “sliding” steady state Uspal et al. (2015). Importantly, there are regions of the phase space in which the rotational effect of wall-induced chemical gradients is stronger than the rotational effect of the hydrodynamic interaction, such that in sum the cap tends to rotate away from the wall. In particular, in the interval in Fig. 10(a), the cap rotates away
from the wall whenever trajectories flow towards larger (see Fig.
2(b)), i.e., to the right of the plot in Fig.
10(a). In this region, rotation away from the
wall by chemical activity is so strong that it can balance the rotational effect of shear
towards the wall, as illustrated in Fig. 2(b). In addition,
whenever the tangent to the trajectories in
Fig. 10(a) is horizontal, i.e., .
Therefore, the swimming activity of the particle can potentially also on its
own produce a steady height . We also note that some trajectories cross
below the
minimum height which we consider. The particles seemingly “crash”
into
the wall. It is likely that many of these trajectories recross the line
and
flow to the attractor if the numerical calculations were extended below ;
however, as
discussed above, the physical approximations inherent in the calculations break down
very close to the wall.^{4}

We now consider the same particle in shear flow. The phase plane with shear rate is shown in Fig. 10(b). As anticipated, there is a rheotactic attractor (filled green circle), which is located at and . Additionally, there is a saddle point (orange cross). As in the case of the “hoverer,” we have chosen a shear rate close to the upper critical value (numerically estimated to be for this surface chemistry), and therefore the saddle point and attractor are very close to each other.

As for Fig. 6, we launch an ensemble of particles from the lateral position of the center of the particles, for various initial director orientations and , and heights , , and . The three-dimensional trajectories of the center of the particles are shown in Figs. 11(a), (b), and (c). The green trajectories are rheotactic. Particles following blue trajectories escape from the surface. Finally, the red trajectories are those which “crash” into the surface, as discussed above. Phase maps indicating the particle behavior as a function of the initial orientation are shown in Figs. 11(d), (e), and (f). For selected rheotactic trajectories, in Fig. 12(a)-(d) we plot the time evolution of the height and of the director components , , and . In Fig. 12(e), we plot for all rheotactic trajectories studied in Fig. 11. As in Fig. 9(e), the decay time for closely agrees with the prediction of Eq. (10), plotted as the dotted line in Fig. 12(e). Finally, a representative rheotactic trajectory is shown in detail in Fig. 13. The radius of curvature of the trajectory in its evolution towards the rheotactic steady state is clearly much larger than for the “hoverer” shown in Fig. 7. As with the hoverer, the dimensionless times in Figs. 12 and 13 can be converted to dimensional, and thus experimentally relevant values via the estimate .

### iii.3 Squirmer

In order to further reveal the general character of our theoretical results and predictions, we consider a spherical “squirmer.” A squirmer has a prescribed surface velocity. It interacts with a bounding wall hydrodynamically, but not chemically, because the surface flow is not driven by a distribution of solute. The prescribed surface velocity is often taken to model the time-averaged motion of cilia on the surface of a micro-organism.

We follow Li and Ardekani Li and Ardekani (2014) in restricting our consideration to the first two squirming modes, given by and in Eq. (11). Li and Ardekani found that stable dynamical attractors exist for various values of and of the Reynolds number Li and Ardekani (2014). (The lowest value of they considered was .) For instance, for and , they found an attractor with and . By construction our numerical method probes the limit . For (and ), we have found an attractor at and . We note that the dynamics of a squirmer near a boundary, including the possibility of moving at a stable height and orientation, was studied also by Ishimoto and Gaffney for Ishimoto and Gaffney (2013).

Here, we consider a squirmer with . We choose this large ratio in order to achieve a strong hydrodynamic interaction with the wall, permitting rheotaxis for a relatively high shear rate , as will be discussed below. (By comparison, for rheotactic states occur only for .) In Fig. 14(a) we show the phase plane for if the director lies in the plane and if there is no external flow. There are two stable dynamical attractors (filled magenta circles) for which the particle moves at a steady height and a steady orientation with its director pointing towards the wall, similar to the “sliding” states we found for catalytic Janus particles Uspal et al. (2015). Due to the rotational symmetry of subproblem II (see Fig. 5 and the corresponding discussion in Sec. III.A), these two attractors are actually the same in the sense that they correspond to the same sliding state, only that the particle slides towards the negative and positive direction for the attractor with and , respectively. In addition, there are two unstable fixed points with