Cross-stream migration of active particles

See allHide authors and affiliations

Science Advances  26 Jan 2018:
Vol. 4, no. 1, eaao1755
DOI: 10.1126/sciadv.aao1755


For natural microswimmers, the interplay of swimming activity and external flow can promote robust directed motion, for example, propulsion against (upstream rheotaxis) or perpendicular to the direction of flow. These effects are generally attributed to their complex body shapes and flagellar beat patterns. Using catalytic Janus particles as a model experimental system, we report on a strong directional response that occurs for spherical active particles in a channel flow. The particles align their propulsion axes to be nearly perpendicular to both the direction of flow and the normal vector of a nearby bounding surface. We develop a deterministic theoretical model of spherical microswimmers near a planar wall that captures the experimental observations. We show how the directional response emerges from the interplay of shear flow and near-surface swimming activity. Finally, adding the effect of thermal noise, we obtain probability distributions for the swimmer orientation that semiquantitatively agree with the experimental distributions.


In microorganisms, such as bacteria, self-propulsion helps them in the efficient exploration of their surroundings to find nutrient-rich areas or to swim away from toxic environments. A common feature among natural microswimmers is an ability to adjust their self-propulsion in response to local stimuli such as chemical gradients, temperature gradients, shear flows, or gravitational fields. Whereas the mechanism of some of these responses, such as chemotaxis, is active, involving the ability to sense gradients and signals, that of others is passive, purely resulting from external forces and torques. For example, gravitaxis in Paramecium is thought to arise solely due to the organism’s fore-rear asymmetry (1). Recently, there has been substantial effort to develop artificial microswimmers, which mimic their natural counterparts in many ways (2). They transduce energy from their local surroundings to engage in self-propulsion and display persistent random walk trajectories, reminiscent of bacterial run-and-tumble behavior (3). They also respond to gravitational fields (4, 5), physical obstacles (611), and local fuel gradients, seeking fuel-rich regions over depleted ones (1214). A number of applications have been envisioned for these biomimetic swimmers, ranging from targeted drug delivery (15, 16) to environmental remediation (17, 18). In many of these applications, it is inevitable that the artificial microswimmers will encounter flow fields while operating in confined spaces. Although their swimming behavior in various complex, static environments has been studied in detail, little is known about their response when they are additionally exposed to flow fields (1922).

In anticipation of this response, one can look to decades of previous work on biological microorganisms in flows [see reviews of Pedley and Kessler (23) and Guasto et al. (24) and the references therein]. The interplay of flow and swimming activity can lead to robust directional response of these microswimmers through a rich variety of mechanisms (2528). For instance, in “gyrotaxis,” bottom-heavy microorganisms, such as algae, adopt a stable, steady orientation and swim against the direction of gravity (2933). Elongated microswimmers, such as sperms, display “rheotaxis,” the ability to orient and swim against the direction of flow (27, 3438). However, much of the work with biological microswimmers has focused on swimming in the bulk, and the detailed role of confining boundaries [for example, in rheotaxis (21, 27, 36, 38, 39)] is an emerging area of research.

Here, using spherical Janus particles as a model system, we demonstrate with experiments and theory that spherical particles swimming near a surface can exhibit surprising and counterintuitive spontaneous transverse orientational order. The particles align almost—but not exactly—in the vorticity direction. As we discuss in detail, this slight misalignment is a key experimental observation that allows us to identify and understand the physical mechanism behind the directional alignment. An inactive, bottom-heavy particle in a flow only rotates as it is carried downstream. However, near-surface swimming activity adds an effective “rotational friction” to the dynamics of the particle because it tends to drive the particle orientation into the plane of the wall. The effective “friction” damps the particle rotation and stabilizes the cross-stream orientation. Finally, the directional alignment, combined with self-propulsion, leads to the cross-stream migration of active particles. The analysis reveals that this mechanism can generically occur for spherical microswimmers in a flow near a surface. Our findings exemplify the complex behavior that can emerge for individual microswimmers from the interplay of confinement and external fields. Moreover, they show that the qualitative character of the emergent behavior is sensitive to the details of the interactions (for example, hydrodynamic interactions) between individual microswimmers and bounding surfaces.


Active particles without external flow

In our experiments, we use silica colloids (1-μm radius, or 2.5-μm radius, where noted) half-coated with a thin layer of Pt (10 nm) as active particles. When the particles are suspended in an aqueous H2O2 solution, the Pt cap catalyzes the degradation of H2O2, whereas the silica half remains inert. The asymmetric distribution of reaction products creates a concentration gradient along the surface of the particle, which induces a phoretic slip velocity, resulting in its propulsion away from the Pt cap. In addition, particle-generated concentration gradients can induce chemiosmotic slip on a nearby bounding surface, giving an additional contribution to particle motility. The details of this mechanism, shown schematically in Fig. 1B, are comprehensively discussed elsewhere (40, 41) and are the subject of ongoing research (4245).

Fig. 1 Schematic illustration of the model system.

(A) A particle (white sphere) with axisymmetric coverage by catalyst (black) is driven by an external shear flow (large gold arrows) in the Embedded Image direction near a planar wall (gray). The particle has a height h above the wall and an orientation vector p, which can be specified by the angles θ and φ. When the particle is active, the cap emits solute molecules (green spheres). (B) Schematic illustrating the self-diffusiophoretic and chemiosmotic mechanisms that drive the motion of a chemical microswimmer. The self-generated solute gradient (green spheres) drives flows localized to thin boundary layers on the particle surface (magenta arrows) and on the nearby wall (blue arrows). The particle is shown in the primed frame (red arrows). This frame corotates with the particle around the Embedded Image axis so that the Embedded Image and the particle orientation vector p is always in the plane spanned by Embedded Image and Embedded Image.

When the silica-Pt particles are suspended in water, they are quickly sedimented to the bottom surface because they are density-mismatched (Embedded Image = 2.196 g/cm3). They orient with their caps down (θ = 0°, where θ is defined in Fig. 1A) due to the bottom heaviness induced by the Pt layer (ρPt = 21.45 g/cm3). Addition of H2O2 introduces activity into the system and changes the orientation distribution of the particles. The particles assume an orientation with the propulsion axis parallel to the bottom surface (θ = 90°; see fig. S1). We have previously reported on the dynamics leading to this change in orientation of the active particles (11). Briefly, we could show that, whereas the hydrodynamic interactions and the bottom heaviness of the particles tend to drive the caps of the particles toward the surface, the wall-induced asymmetry of the distribution of the chemical product and the chemiosmotic flow along the substrate (see Fig. 1B) tend to have the opposite effect, leading to a stable orientation at θ ≈ 90°. Once parallel to the surface, the particles are confined to a single plane of motion, where they propel away from their Pt caps (Fig. 2, C and D) with a typical speed Vp = 6 μm/s. Because of the contrast between the dark Pt hemisphere and the transparent silica, we can measure the angular orientations of these particles (see Materials and Methods and fig. S2 for details of the tracking process). Within the two-dimensional (2D) plane, the particles have no preferred directionality and are diffusive on long time scales.

Fig. 2 Effect of imposed flow on passive and active colloidal particles.

(A) Passive silica-Pt particles in a flow (R = 1 μm). The plot shows the dependence of the rotation time on flow velocity. V* is the translational velocity of the passive particle and is used to characterize the flow rate. The red dashed line is a theoretical scaling derived in section S4 and fitted to the data. Inset: Time-lapse images of a passive particle rolling in a flow. (B) Tracked trajectories of passive particles in a flow (V* = 14 μm/s). (C) Optical microscopy image capturing the distribution of orientations of active particles (Vp ≈ 6 μm/s) without an imposed flow. (D) Tracked trajectories of active particles without a flow. (E) Same system as (C) with an external flow imposed, V* = 14 μm/s. (F) Tracked trajectories of active particles in a flow (V* = 14 μm/s). (G) Angular probability distributions of active colloidal particles in the absence of an imposed flow, V* = 0 μm/s (green), at V* = 14 μm/s (blue) and at V* = 24 μm/s (red). Inset shows the angular evolution of two different active colloids at imposed flow rates of V* = 14 μm/s (blue) and V* = 24 μm/s (red).

Passive particles in external flow

Now, we seek to characterize the behavior of these silica-Pt particles in an imposed flow. Initially, a suspension of the particles in water is introduced in a square glass capillary (1 mm; VitroCom) connected to a computer-controlled microfluidic pump (MFCS-EZ, Fluigent). We allow the particles to sediment to the bottom surface before we impose any external flows. The desired flow rate in the capillary is maintained by using a flow rate monitor (Flowboard, Fluigent), which is in a feedback loop with the microfluidic pump. We begin by imposing a flow of water (no activity) in the x direction. Close to the nonslipping capillary surface, the flow velocity varies linearly as Embedded Image, and the particles, which are sedimented near the surface, experience a shear flow (see section S3 for a calculation of the flow profile) (46). In terms of their translational behavior, we observe that the particles act as “tracers” and translate in nearly straight lines along the direction of flow (Fig. 2B). The translational velocity of these particles is proportional to the imposed flow rate. We use the translational velocity V* of these inactive “tracer” particles to characterize the flow rate (39). Before we start the flow, the particles are all in the cap-down orientation (θ = 0°). The shear flow induces a torque on the particles, and they rotate around the axis of flow vorticity Embedded Image as they translate in a flow (see movie S1). The rotation speed of the particles is also dependent on the flow rate (Fig. 2A), with higher flow rates leading to faster rotation. Through a simple model for particle rolling developed in section S4, we predict the rotational period Embedded Image, where a and b are fitting parameters. This relation shows excellent fit to the data (Fig. 2A, red curve), and from the fitted a, we extract Embedded Image.

Active particles in external flow

We then start a flow of H2O2 to introduce activity into the system. We notice that the dynamics of particle behavior in a flow is markedly influenced in the presence of activity. First, the particles stop rolling and reach a stable orientation parallel to the bottom surface (θ = 90°). This is similar to a previously observed behavior for the particles in the absence of flow (11). More surprisingly, the particles also evolve to a stable orientation that is nearly perpendicular to the direction of imposed flow (φ ≈ 90° or −90°; Fig. 2E). While the particles continue to translate in x due to the imposed flow, they also have the self-propulsion velocity Vp away from their Pt caps. A combination of these effects results in the cross-streamline migration of self-propelled particles, that is, migration of particles in the y direction, perpendicular to the flow along the x direction. Typical trajectories of cross-stream migrating particles are presented in Fig. 2F. Further, we observe that the stability of cross-stream migration (due to the stability of the steady orientation angle φ* perpendicular to the direction of flow) is dependent on the flow rate, with higher flow rates resulting in a stronger alignment effect. The inset of Fig. 2G shows the angular evolution of two self-propelled particles in imposed flows of V* = 14 and 24 μm/s. The deviations away from the φ ≈ 90° positions occur more frequently and at larger amplitude than for particles in lower flow rates.

To study this effect at a population scale, we flow a suspension of self-propelled particles with H2O2 and record the angular orientations and positions of every particle in each frame, which allows us to determine the probability distribution of φ in the system. These are plotted in Fig. 2G for two different flow rates and compared to the system of self-propelled particles without any imposed flow. In the absence of flow, the distribution is nearly flat, indicating the lack of preference for any orientation φ, and at long time scales, the particle behavior is purely diffusive. However, in the case of imposed flow, we observe distinct peaks that appear at φ ≈ 90° and − 90°. These correspond to the particles exhibiting the cross-stream behavior. These peaks also become sharper when we increase the flow rate, as can be seen in the distributions for V* = 14 and 24 μm/s. In both cases, closer observation of the angular probability distributions reveals a small bias of orientations in the direction of flow for particles migrating across flow streamlines. Because the particles are subject to Brownian fluctuations, they can also occasionally orient with or against the flow (φ ≈ 0° or 180°). This leads to an intermittent state where the particles “tumble” in the direction of flow before recovering their φ ≈ ±90° orientation and the cross-streamline behavior. The recovered orientation of these particles can be different from their initial states, changing the direction of particle migration (for example, from +y to −y; see fig. S5).

Apart from the imposed flow rate V*, we find that the stability of the φ* is also dependent on the particle radius. Using particles of larger radius (2.5 μm), we show that the distribution of φ is narrower around φ = ±90°, as compared to the distribution for R = 1 μm particles for identical Vp and V* (fig. S6). Because the effect of Brownian noise is significantly lower for the larger particles, we also observe that they seldom switch their migration direction within the width of our capillary (fig. S7).

We can further control the behavior of particle migration by tuning the self-propulsion velocity Vp of the particles. First, we find that higher propulsion velocities damp the fluctuations around φ* due to higher activity (see fig. S8). Second, the propulsion velocity also controls the “slope” of the cross-stream migration. To quantify this, we define α to be the offset between the orientation vector p and the tracked velocity vector v of the particle. For a self-propelled particle in the absence of flow, the particles translate in the direction of the orientation vector p as the particles propel away from their Pt caps (that is, α = 0°), as shown in Fig. 3 (A and B). However, with an imposed flow, the translational direction differs from the orientation vector, and the offset α, for a given flow rate, is determined by the Vp of the particles. For particles with Vp = 6 μm/s, the probability distribution function of α has a peak around −64°, whereas for a particle with Vp = 3 μm/s, α is peaked at −83° (Fig. 3, A, C, and D). The offset α can be rationalized as a contribution of the transverse propulsion induced by the cross-stream orientation and the longitudinal advection by flow. The α should then simply be given by α = arctan(V*/Vp). Substituting the experimental values for Vp and V*, we get α = 66. 9° and 78. 3°, which agree reasonably well with the observed values. This clear separation in peaks of α for different propulsion velocities could eventually be used for the separation of particles based on activity.

Fig. 3 The dependence of cross-stream migration on Vp.

(A) Probability distributions of α for a particle in the absence of a flow with Vp = 6 μm/s and with an imposed flow corresponding to V* = 14 μm/s for particles of Vp = 6 and 3 μm/s. (B to D) Tracked orientation vectors (line segments) and instantaneous velocities (arrows) for trajectories at different flow and self-propulsion speeds.

Construction of theoretical model

The experimental observations can be qualitatively captured and understood within a generic mathematical model of swimming near a surface in external flow (21). In this model, we construct dynamical equations for the height and orientation of a heavy spherical microswimmer by exploiting physical symmetries and the mathematical linearity of Stokes flow. As detailed below, mathematical analysis of these equations reveals that qualitatively distinct steady-state behaviors, including cross-stream migration, can emerge from the interplay of external shear flow, near-surface swimming, and gravity. We note that the analysis does not depend on a particular mechanism of self-propulsion (for example, self-diffusiophoresis, self-electrophoresis, or mechanical propulsion by motion of surface cilia); it is, in that sense, generic. Accordingly, our major theoretical findings are (i) that there is a physical mechanism that can produce the surprising transverse orientational order observed in experiments and (ii) that this mechanism can generically occur for a spherical, axisymmetric microswimmer exposed to flow near a bounding surface.

We now proceed to construction and analysis of the dynamical equations. As discussed in Materials and Methods, the flow in the suspending fluid is characterized by a low Reynolds number and hence governed by the Stokes equations. Because these equations are linear, the contributions of external flow (f), gravity (g), and swimming (s) to the particle translational and angular velocities can be calculated independently and superposed: U = U(f) + U(s) + U(g) and Ω = Ω(f) + Ω(s) + Ω(g). The velocity of the orientation vector is determined by Embedded Image. We will generally describe p in Cartesian coordinates but will sometimes find it useful to use the spherical coordinates p = (sin(θ) cos(φ), sin(θ) sin(φ), cos(θ)). Note that |p| = 1.

To calculate the contribution of the external flow to the particle velocity, we consider a neutrally buoyant and inactive sphere of radius R in shear flow at a height h above a uniform planar wall. Shear spins the particle around the vorticity axis Embedded Image, as shown in the left panel of Fig. 4. This is because the flow is faster near the upper surface of the particle (that is, the surface farther away from the wall) than near the bottom surface. Accordingly, the particle has angular velocity Embedded Image. In addition, the particle is carried downstream with velocity Embedded Image. The functions f(h/R) and g(h/R) represent the influence of hydrodynamic friction from the wall (47); they are provided in section S8. Because of the spinning motion of the particle, the tip of the vector p traces a circle in the shear planeEmbedded Image(1)

Fig. 4 Schematic of the contributions to planar alignment steady state.

(Left) A neutrally buoyant, inactive sphere driven by an external shear flow (gold arrows) near a planar wall. Because the external flow spins the particle around the vorticity axis (that is, the normal to the shear plane), the tip of the particle orientation vector p traces a circular path (magenta). (Middle) An active, heavy sphere in quiescent fluid. Because the wall is uniform and the geometry and activity profile of the particle are axisymmetric, the tip of the orientation vector can only rotate in the Embedded Image direction, that is, directly toward or away from the wall. (Right) Planar alignment steady state, where Embedded Image, but the particle orientation vector has nonzero components Embedded Image and Embedded Image in the flow and vorticity directions, respectively. For planar alignment, the component Embedded Image in the flow direction can be either upstream (Embedded Image) or downstream (Embedded Image), as determined by the function Embedded Image; the downstream case is shown. All contributions to Embedded Image are in the Embedded Image direction (see Eqs. 2 to 5). At a certain angle φ*, all contributions to Embedded Image balance, as shown by the arrows, so that Embedded Image. Note that this fixed point always occurs in pairs related by mirror symmetry across the shear plane; we show the state with py > 0.

The component py is a constant determined by the initial orientation of the particle. The radius of the circular orbit is Embedded Image, and the speed of the tip Embedded Image is proportional to the circle radius.

Now, we consider a heavy, active particle moving in quiescent fluid (no external flow) near the same surface. The particle has an axisymmetric (ax) geometry and surface activity profile; accordingly, we define U(ax)U(s) + U(g) and Ω(ax) ≡ Ω(s) + Ω(g). This system, depicted in the middle panel of Fig. 4, has a plane of mirror symmetry defined by Embedded Image and p. Accordingly, translation of the particle is restricted to this plane, and p is restricted to remain within this plane (that is, Embedded Image), although p can rotate toward or away from the wall (Embedded Image). We write Embedded Image, where we have defined a new, “primed” frame that has Embedded Image and Embedded Image in the plane of mirror symmetry, with Embedded Image. This frame is convenient for calculations because Ω(ax) is strictly in Embedded Image (see Fig. 1B). The function Embedded Image incorporates the effect of bottom heaviness and interactions with the wall (for example, hydrodynamic interactions) that originate in swimming activity. Both Embedded Image and U(ax)(θ, h/R) depend on the particle design and the model for the propulsion mechanism. To show that our analysis is generic in the sense discussed above, we will leave them unspecified until later in the text. Because the particle is axisymmetric and the wall is uniform, these functions have no φ dependence.

Now, we superpose all contributions and obtain equations for Embedded Image and Embedded Image in the following formEmbedded Image(2)Embedded Image(3)Embedded Image(4)Embedded Image(5)Note that the axisymmetric contributions are evaluated in the primed frame, which is defined to corotate with the particle, and we used Embedded Image. The components of the translational velocity U in the x and y directions are determined by p and h asEmbedded Image(6)Embedded Image(7)

Steady-state solutions

Equations 2 to 7 fully describe a deterministic active sphere in external shear flow near a planar surface. Now, we look for the fixed points (p*, h*) of Eqs. 2 to 5. A fixed point is a particle configuration in which the particle translates along the wall with a steady height and orientation, that is, Embedded Image and Embedded Image (21). We find that the system has three fixed points, shown schematically in fig. S10. Of these three, planar alignment shows excellent qualitative agreement with the experiment observations: The particle orientation is within the plane of the wall (Embedded Image) and has nonzero components in both the flow and vorticity directions. These criteria are not satisfied by the other two fixed points. Therefore, slight misalignment from the vorticity axis is a key experimental observation that discriminates between the steady states predicted by theory. For planar alignment, the streamwise component Embedded Image isEmbedded Image(8)which can be either downstream (Embedded Image), in agreement with our experimental observations, or upstream (Embedded Image), depending on the sign of Embedded Image. Because of the mirror symmetry of Eqs. 2 to 5 with respect to py = 0, this fixed point always occurs in pairs Embedded Image. Further mathematical details are provided in the Supplementary Materials.

To understand why planar alignment is a steady state, we consider the conditions for the contributions of shear flow, activity, and gravity to the angular velocity to cancel out at some p*. We recall that the axisymmetric contributions to Embedded Image are in the Embedded Image direction. For most orientations θ and φ, Embedded Image will have a component in the vorticity direction Embedded Image. However, from Eq. 1, we see that the shear contribution never has a Embedded Image component. As a way out of this dilemma, we see from the definition of Embedded Image that, when θ = 90° (that is, pz = 0), the Embedded Image and Embedded Image components of θ both vanish, and Embedded Image is strictly in the Embedded Image direction. Likewise, when θ = 90°, Embedded Image is strictly in Embedded Image. This suggests the possibility that all contributions can cancel out for θ* = 90° and some unknown φ*. Now, we consider the role of φ. The axisymmetric contributions Embedded Image have no dependence on φ, as discussed previously. The contribution from shear does depend on φ, as can be seen by substituting pz = 0 and px = cos(φ) into Eq. 1. Therefore, the sign and magnitude of the shear contribution can be “tuned” via the angle φ to cancel out the axisymmetric contributions to Embedded Image (provided that these are not too large). Finally, we consider the height of the particle. Because shear does not contribute to vertical motion, the particle must obtain a steady height (Embedded Image) through the combined effects of near-surface swimming activity and gravity.

As an additional note, if we assume that the steady height of the particle is not significantly affected by flow rate, then Eq. 8 predicts that the steady-state orientation approaches the vorticity axis as the flow rate is increased, that is, cos(φ*) ~ 1/V*. Accordingly, we perform further experiments with 5-μm-diameter particles, which allow for excellent resolution of their orientation (see fig. S2). We performed the experiments at six different flow rates in the range accessible with the current experimental setup, V* = 5 to 75 μm/s, for particles with Vp = 6 μm/s. In fig. S11, we plot cos|φ*| against V* and find that the data recover the predicted asymptotic behavior.

Linear stability analysis

It is not enough to find a fixed point that corresponds to experimental observations; we must also consider its stability. Because, experimentally, the active particles spontaneously adopt a steady cross-stream orientation, the fixed point should be a stable attractor for active particles. Second, because the inactive particles are observed to continuously rotate in the experiments, the fixed point should be associated with closed orbits for inactive particles. In the following, we carry out a linear stability analysis and show that our model meets both criteria, provided that a certain condition on the particle/wall interaction is satisfied.

We define the generalized configuration vector f ≡ (p, h/R) and Embedded Image as a small perturbation away from the fixed point: Embedded Image. We obtain the linearized governing equations Embedded Image, where the Jacobian matrix Embedded Image is given in detail in the Supplementary Materials. As a simplifying approximation, we take the particle to have a constant height h = H so that Embedded Image. This approximation is motivated by the experimental observation that the particles never leave the microscope focal plane. Having made this approximation, we can find an intuitive analogy between our system and a damped harmonic oscillator (see the Supplementary Materials)Embedded Image(9)

The motion of the orientation vector resembles a damped harmonic oscillator with intrinsic frequency Embedded Image. Recalling that inactive, bottom-heavy particles were observed to continuously rotate in the experiments, we consider the case in which Embedded Image. We note thatEmbedded Image(10)where we have used Embedded Image and Embedded Image. Therefore, the dissipative term in Eq. 9 is zero, and the orientation p has a continuous family of closed orbits centered on p*. Now, we consider active, bottom-heavy particles, for which Embedded Image. The fixed point is a stable attractor ifEmbedded Image(11)

This condition has the following interpretation: If the particle orientation is perturbed out of the plane of the wall, then the active contribution to rotation responds (increases or decreases) so as to oppose the perturbation. If this condition is satisfied, then the activity induces an effective friction that damps oscillations of the orientation p, driving attraction of p to p*. Satisfaction of the condition depends on the details of the interactions between the particle and the wall that originate in swimming activity. In turn, these details depend on the character of the self-propulsion mechanism [see, for example, the comparison of pushers and pullers in the study of Spangolie and Lauga (7)].

Application to a self-phoretic particle

Up to this point, we have proceeded without specifying a model for the particle composition and self-propulsion, and we derived and analyzed Eqs. 2 to 5 in general terms. Now, we seek to obtain illustrative particle trajectories by numerical integration. Accordingly, hereafter, we calculate the various terms in Eqs. 2 to 5 by using a simple, well-established model of neutral self-diffusiophoresis in confinement (9, 11, 41, 4850). In the model, the particle has a hemispherical catalytic cap (Fig. 1B, black), and the orientation vector p points from its catalytic pole to its inert pole. The particle emits solute molecules (that is, oxygen) at a constant, uniform rate from its cap, leading to self-generated solute gradients in the surrounding solution. These gradients drive surface flows on the particle (Fig. 1B, magenta arrows) and on the wall (blue arrows), leading to directed motion of the particle. We note that our generic model does not account for effects (for example, ionic or electrokinetic) specific to the detailed self-phoretic mechanism (4245). Further details are provided in Materials and Methods.

We consider some illustrative examples using dimensionless parameters comparable to those in the experiments (see Materials and Methods for the definition and estimation of parameters). We fix the particle height as H/R = 1.2. In Fig. 5A, we show a particle trajectory in a shear flow with strength Embedded Image and initial orientation θ0 = 30° and φ0 = 315°. The particle rotates so that its inert face points largely in the Embedded Image direction but with a small downstream orientation (px > 0). Notably, this slight downstream orientation agrees with experimental observations. If the particle is inactive (but still bottom-heavy), then from the same initial orientation, the particle simply translates in the flow direction (Fig. 5B), rotating as it does so.

Fig. 5 Numerically obtained trajectories of inactive and active particles in flow.

(A) Numerically computed trajectory of a passive bottom-heavy particle near a wall (gray) and driven by an external shear flow (gold arrows) with dimensionless strength Embedded Image (concerning dimensionless parameters, see Materials and Methods). The initial condition of the particle is θ0 = 30° and φ0 = 315°, and the particle height is fixed as h/R = 1.2. The particle rotates and is carried downstream by the external flow with no cross-streamline migration. (B) Numerically computed trajectory of an active Janus particle with the external flow strength, particle materials, and initial conditions as in (A). The particle rotates so that its inert face points largely in the Embedded Image direction, with a slight downstream orientation (px > 0). With this steady orientation, the particle swims across flow streamlines as it moves downstream. (C) Probability density function (pdf) for φ, for a catalytic Janus particle in a shear flow with Embedded Image. (D) Probability density function for the same particle with Embedded Image.

Because h = H is taken to be constant, the instantaneous value of p completely determines the instantaneous velocity of the particle. In Fig. 6 (A to C), phase space trajectories are shown on the unit sphere |p| = 1. The initial orientation in Fig. 5 (A and B) is indicated by magenta circles in Fig. 6 (A and B, respectively). We see that an inactive particle (Fig. 6A) has a continuous family of closed orbits in p. When the particle is active, these oscillations are damped, and trajectories in the phase space are attracted to the fixed point (Fig. 6, B and C). In Fig. 6C, we show the effect on the structure of trajectories as the flow strength is increased. For stronger flows, the approach to the fixed point is more oscillatory. In Fig. 6E, we plot the component pz (black line) for the inactive particle in Fig. 5A. In Fig. 6F, we show pz for the active particle in Fig. 5B. We can also obtain analytical solutions to the linearized equations, characterized by a few input parameters that are evaluated numerically (see the Supplementary Materials). These solutions are shown as dashed red lines in Fig. 6 (E and F) and agree well with the numerical data. In addition, we note that if the particle height is allowed to change, then one can still obtain the phenomenology studied here of attraction to a planar alignment steady state; an example is given in section S13.

Fig. 6 Phase portraits of inactive and active particles in flow.

(A to C) Phase portraits on the sphere |p| = 1 for an inactive, bottom-heavy particle (A) in shear flow with Embedded Image and for an active, bottom-heavy particle (B and C) in a flow with Embedded Image (B) and Embedded Image (C). The red dashed lines indicate the plane of mirror symmetry py = 0. The cyan circle in (A) indicates the center of oscillatory motion. Green circles in (B) and (C) indicate stable fixed points (“attractors”). For (A) and (B), the magenta circles show the initial conditions for the trajectory in Fig. 5A (A) and that in Fig. 5B (B). The parameters characterizing the particle heaviness and activity are given in Materials and Methods. (D) Vector field on the unit sphere representing the contribution of activity to the motion in (B) and (C). The green dashed line indicates pz = 0. (E) Oscillation in pz with time for the inactive particle in Fig. 5A, also corresponding to the magenta circle in (A). The black dashed line gives the result obtained by numerical integration. The red dashed line gives the analytical solution for the linearized equations. The slight disagreement between them is due to the large amplitude of the oscillation, for which the effect of the nonlinearity in Eqs. 2 to 4 is important. (F) For the active particle in Fig. 5B, with initial conditions corresponding to the magenta circle in (B), the oscillation in pz decays with time, and the orientation p eventually approaches the stable fixed point p*. The analytical solution is shown for the later part of the decay, for which pz has a small amplitude and is well described by linear theory.

The phase portraits also provide an intuitive way to understand the stability condition expressed in Eq. 11. Consider the bottom-heavy, inactive particle with the phase portrait shown in Fig. 6A. How does adding activity transform this portrait into that shown in Fig. 6B? The contribution of activity to Embedded Image is shown as the vector field in Fig. 6D. Because this vector field is small but nonzero on the equator pz = 0, adding it to the portrait in Fig. 6A will shift the center of oscillation (cyan circle) slightly toward the axis px = 0, producing the green circle in Fig. 6B. More significantly, the addition of this vector field will destabilize continuous oscillatory motion in the following way. Consider the closed trajectories in Fig. 6A that are closest to the cyan circle. In the neighborhood of the equator, activity always “pushes” the vector p toward the equator (here, neglecting the small value on the equator, the main effect of which is to shift the fixed point). Therefore, over each period of oscillation, the orientation vector will get slightly closer to the fixed point, transforming the closed circular orbits into decaying spiral orbits.

Effect of thermal noise

All of the preceding analysis assumed deterministic motion. As a further exploration, we consider the effect of thermal noise on the particle orientation by performing Brownian dynamics simulations (see Materials and Methods). In Fig. 5C, we show the probability distribution functions for the φ obtained for a Janus particle driven by shear flow at dimensionless inverse temperature Pep = 500 (see Materials and Methods for the definition) and shear rate Embedded Image. The distribution is symmetric and has two peaks near the steady angle φ* = ±80.9° predicted by the deterministic model (Fig. 5C). For a higher shear rate Embedded Image, we find that the peaks in φ are sharper (Fig. 5D), in qualitative agreement with the experiments, with the peaks shifted near the deterministic prediction for this shear rate, φ* = ±88.2°. Additional probability distributions for the components of p are given in fig. S14.


Here, we use catalytic Janus particles as a well-controlled model experimental system to study spherical active particles in confined flows. We demonstrate that spherical active particles near surfaces, when exposed to external flows, can exhibit robust alignment and motion along the cross-stream direction. Our model reveals how this behavior arises from the interplay of shear flow and swimming in confinement. The steady orientation is determined by a balance of contributions to the angular velocity of the particle from shear flow, bottom heaviness, and swimming near a planar substrate. Near-surface swimming introduces an effective friction opposing the rotation of the particle away from the preferred orientation. The mechanism is generic in the sense that it can occur for any spherical microswimmer with axisymmetric actuation and is not specific to a particular mechanism of propulsion (for example, chemical or mechanical). As a consequence of the alignment, the particles migrate across the streamlines of the external flow as they are carried downstream.

To the best of our knowledge, our results are the first to demonstrate that swimmer/surface interactions (for example, hydrodynamic interactions) can drive a rich directional response of spherical particles to external flows. This is in contrast with previous works on natural microswimmers (for example, bacteria), where complex body shapes and flagellar beat patterns were implicated in directional response (25, 27). In addition, we have obtained semiquantitative agreement between particle orientational statistics obtained from the experiments and from the theoretical model.

For lab-on-a-chip devices that use continuous flows and artificial microswimmers, our findings imply that the microswimmers would have a tendency to migrate to the confining side walls of the device. Our findings additionally raise the possibility that in dense suspensions of microswimmers, for which fluid flows are self-generated, the collective behavior of the suspension may be sensitive to the detailed interactions between individual microswimmers and bounding surfaces.


Sample preparation

The Janus particles were obtained by electron-beam deposition of a Pt layer on a monolayer of silica colloids. The monolayers were prepared either by a Langmuir-Blodgett (LB) method or a drop-casting method. For the LB method, silica colloids (R = 1 μm; Sigma-Aldrich) were first surface-treated with allyltrimethoxysilane to make them amphiphilic. A suspension of these particles in chloroform-ethanol mixture [80%:20% (v/v)] was carefully dropped onto the LB trough and compressed to create a closely packed monolayer. The monolayer was then transferred onto a silicon wafer at a surface pressure of 20 mN/m. The silicon wafer was then shifted to a vacuum system for the electron-beam deposition of a thin layer of Pt (10 nm) at 10−6 mmHg. The Janus particles were released into deionized water using short ultrasound pulses. The suspension of Janus particles in water was stored at room temperature. The monolayers of R = 2.5 μm silica particles (Sigma-Aldrich) were prepared by drop-casting the suspension of colloids onto an oxygen plasma–treated glass slide. The plasma treatment was used to make the glass slide hydrophilic and ensure uniform spreading of the particle suspension. The solvent was subsequently removed by slow evaporation. The Pt deposition step was identical to the one used for R = 1 μm particles.


Particle tracking was performed using an automated tracking program developed in-house. The Python-based program uses OpenCV library for image processing and NumPy for data handling. In source videos filmed in grayscale, each frame was first cleaned of noise by using blurring techniques, which substituted each pixel with an average of its surroundings. The particles were then separated from the background by using either of the two segmentation methods: threshold and gradient. In the threshold method, given a grayscale image img(x,y) and a threshold value T, this operation resulted in a binary image out(x,y) given byEmbedded Image(12)

The gradient method was used for images with irregular brightness or when the particles were hard to distinguish from the background. In the first step of this method, the gradient of the image Embedded Image is approximated by convolving the original frame with a Sobel operator. This results in two images, one that is the derivative along the x axis and another that is along the y axis. These images are then thresholded and joined together to obtain the segmented image. The final result has the edges of the detected particles.

The center of each particle was approximated as the center of mass of the contours obtained after segmentation. Particle trajectories were calculated using Bayesian decision-making, linking every particle center with the previous closest one. Intermediate missing positions, if any, were interpolated using cubic or linear splines.

To calculate the p vector, a line of predefined radius (approximately equal to the particle radius) was drawn through the center at a test angle. The SD of the pixel values along the test line was calculated and stored. The process was repeated multiple times, and the SDs were compared. The test angle with the least SD was assumed to correspond to the separator between the silica and Pt halves and the vector orthogonal to it, the orientation vector (see fig. S2).

Theoretical calculation of particle velocity

Here, we presented our model for the calculation of the particle velocities U and Ω as a function of h and p. We took the instantaneous position of the particle to be xp = (xp, yp, h) in a stationary reference frame. The catalytic cap emitted a product molecule at a constant and uniform rate κ. We took the solute number density c(x) to be quasi-static, that is, it obeyed the equation ∇2c = 0 with boundary conditions Embedded Image on the cap, Embedded Image on the inert region of the particle surface, and Embedded Image on the wall. Here, D is the diffusion coefficient of the solute molecule, x is a location in the fluid, and Embedded Image and is the normal vector pointing from a surface into the liquid. For each instantaneous configuration (h, p), this set of equations can be solved for c(x), for example, numerically by using the boundary element method (21, 51).

We took the velocity u(x) in the fluid solution to obey the Stokes equation −∇P + η∇2u = 0, where η is the solution viscosity and P(x) is the pressure in the solution. In addition, the fluid was incompressible so that ∇ ⋅ u = 0. The Stokes equation is a linear equation. Therefore, the contributions to U and Ω from various boundary conditions for u, and from various external forces and torques, can be calculated individually as the solution to separate subproblems and then superposed. Hence, we write U = U(f) + U(ax) and Ω = Ω(f) + Ω(ax), where (f) indicates the contributions of the external flow, and the axisymmetric (ax) contributions are from gravity (g) and swimming (s) activity: U(ax) = U(g) + U(s) and Ω(ax) = Ω(g) + Ω(s). For each subproblem, the fluid is governed by the Stokes equation and incompressibility condition.

We first considered the subproblem for the contribution of the external flow. The fluid velocity was subject to no-slip boundary conditions u = 0 on the planar wall and u = uext + U(f) + Ω(f) × (xxp) on the particle. Here, uext is the external flow velocity, Embedded Image. In addition, the particle was free of external forces and torques, closing the system of equations for U(f) and Ω(f). The solution of this subproblem is well known [see, for instance, the study of Goldman et al. (47)].

The two axisymmetric subproblems were calculated in the primed frame corotating with the particle (Fig. 1B). For the subproblem associated with particle activity, we used the classical framework of neutral self-diffusiophoresis (9, 11, 40, 41, 48). In this subproblem, the self-generated solute gradients drive surface flows on the wall and the particle surface, vs = −b(xs′)∇||c(x′), where Embedded Image, and xs′ is a location on a surface (in the primed frame). The “surface mobility” b(xs′) encapsulates the details of the molecular interaction between the solute and the bounding surfaces. We wrote the boundary conditions u = U(s) + Ω(s) × (x′ − xp′) + vs(xs′) on the particle and u = vs(xs′) on the wall. Again, specifying that the particle was force- and torque-free closed the system of equations for U(s) and Ω(s). This subproblem can be solved numerically using the boundary element method (21, 51). For the subproblem associated with gravity, we used the “eggshell” model of Campbell and Ebbens (4) for the shape of the cap, taking the cap thickness to vary smoothly from zero at the particle “equator” to a maximum thickness of t at the active pole. Details concerning this subproblem, which is solved using standard methods (11), are provided in the Supplementary Materials.

In the following, we specify the parameters characterizing the system. We chose to take the inert and catalytic regions of the particle to have different surface mobilities, binert and bcap, with binert/bcap = 0.3 and bcap < 0. For this parameter, a neutrally buoyant Janus particle, when it was far away from bounding surfaces, moves in the p direction (that is, away from its cap) with a velocity Ufs = 13/80 U0, where U0 = |bcap|κ/D (52). The wall was characterized by a surface mobility bw. We chose bw/bcap = −0.35. These surface mobility ratios were chosen to be similar to those used in a previous work [where we had binert/bcap = 0.3 and bw/bcap = −0.2 (11)], and to give a slightly downstream steady orientation. We nondimensionalized length with R, velocity with U0, and time with T0 = R/U0. To nondimensionalize the gravitational and shear contributions, we must estimate U0 in real, dimensional units. Rather than to calculate U0 directly, which requires estimates for κ and |bcap|, we used the expression Ufs = 13/80 U0 (52). Knowing that, experimentally, the particle characteristically moves at Ufs ≈ 5 μm/s when it is far from surfaces, we obtained U0 ≈ 30 μm/s. The parameters describing the heaviness of the particle are given in the Supplementary Materials. Finally, we considered the shear rate Embedded Image. This was not known experimentally but could be roughly estimated by considering inactive particles to act as passive tracers. From Goldman et al. (47), the velocity of a spherical particle driven by shear flow near a wall is Embedded Image. Experimentally, inactive particles in a flow were observed to move with V* ~ 10 μm/s. The particle height h is difficult to observe experimentally, but we took it to be set by the balance of gravity and electrostatic forces. Hence, the particle/wall gap δ was on the order of a Debye length λD ~ 0.1R so that h ~ 1.1R. For h/R ≥ 1.05, the factor g(h/R) ≈ 1. Therefore, we estimated Embedded Image and a typical dimensionless shear rate to be Embedded Image. As a reminder, our aim was to establish semiquantitative agreement with experiments, and therefore, we sought only an order of magnitude accuracy in the dimensionless parameters.

In assuming the concentration field to be quasi-static, we neglected the advective effects on the solute field by the external shear flow and by the finite velocity of the particle. These approximations are valid for small Peclet numbers Pe ≡ UfsR/D and Embedded Image. At room temperature, the diffusion coefficient of oxygen was D ~ 4 × 10−9m2/s (53) so that Pe ≈ 0.003 and Embedded Image. Furthermore, in taking the fluid velocity u to be governed by the Stokes equation, we neglected fluid inertia. This approximation is justified for low Reynolds number, Re = ρfluidU0R/η, where η is the dynamic viscosity of the solution. Using η ~ 10−3Pa s for water, we obtained Re ≈ 10−4.

Effect of thermal noise

The particle Peclet number Embedded Image characterizes the relative strengths of self-propulsion and translational diffusion of the particle. Here, Embedded Image is the translational diffusion coefficient of the particle in free space. For a particle with R = 2.5 μm in water at room temperature (kBTroom ~ 4 × 10−21J) with Ufs ≈ 5 μm/s, we estimated Pep ≈ 880. To perform Brownian dynamics simulations, we adapted the Euler-Maruyama integration scheme introduced by Jones and Alavi (54) and later presented by Lisicki et al. (55), which explicitly included the effects of the wall on diffusion. Our principal modification to the method was the inclusion of deterministic contributions to Embedded Image (Eqs. 2 to 4). Further details are provided in the Supplementary Materials and by Lisicki et al. (55).


Supplementary material for this article is available at

section S1. Activity-induced reorientation of particles

section S2. Method to detect the in-plane angle of Janus particles

section S3. Flow profile in square capillary

section S4. Estimation of shear rate from particle rotation

section S5. Interim states in cross-stream migration of active particles

section S6. Effect of particle size on orientational stability

section S7. Effect of propulsion velocity on orientational stability

section S8. Contributions of shear to the equations of motion

section S9. Fixed points of governing equations

section S10. Linear stability of planar alignment

section S11. Steady angle of a particle as a function of flow rate

section S12. Calculation of gravitational contribution to particle motion

section S13. Motion of particle in three dimensions

section S14. Brownian dynamics simulations

section S15. Legends for movies S1 to S3

fig. S1. Experimentally measured probability distribution of θ for active particles.

fig. S2. Illustration of method to detect the in-plane angle of Janus particles.

fig. S3. Flow profile in square capillary.

fig. S4. Rotation time τ as a function of speed V* for inactive Janus particles in a flow.

fig. S5. Interim states in cross-stream migration of active particles.

fig. S6. Probability distribution function of φ for particles of R = 1 and 2.5 μm at Vp 6 μm/s and V* ≈ 24 μm/s.

fig. S7. Sample trajectories of R = 2.5 μm particles show very little deviation from their preferred cross-stream orientation at Vp ≈ 6 μm/s and V* ≈ 24 μm/s.

fig. S8. Fluctuations in φ obtained from five different particles for two values of Vp at V* ≈ 24 μm/s.

fig. S9. The functions f(h/R) and g(h/R) as obtained by Goldman et al. and with the boundary element method.

fig. S10. Schematic illustration of the three fixed point solutions to Eqs. 2 to 5 in the main text.

fig. S11. The mean px = |cos (φ*)| plotted as a function of V*.

fig. S12. Probability distribution of |φ| plotted for two different flow velocities shows a clear shift of the peak position φ* toward 90° at higher flow rates.

fig. S13. 3D trajectory of a bottom-heavy, chemically active particle in a shear flow.

fig. S14. Probability distribution functions for components of the particle orientation vector p, obtained from stochastic numerical simulations.

table S1. Comparison of f(h/R) and g(h/R) as calculated by Goldman et al. and in this work using the boundary element method.

movie S1. Silica-Pt active Janus particles in the absence of any external shear flow.

movie S2. Silica-Pt Janus inactive particles in a shear flow (no hydrogen peroxide).

movie S3. Silica-Pt Janus active particles in a shear flow.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank M. N. Popescu for the helpful discussions. Funding: W.E.U. acknowledges financial support from the Deutsche Forschungsgemeinschaft grant no. TA 959/1-1. S.S., J.S., and J.K. acknowledge the DFG grant no. S.A 2525/1-1. S.S. acknowledges the Spanish Ministry of Economy, Industry, and Competitiveness for grant CTQ2015-68879-R (MICRODIA). This research also has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement no. 311529. Author contributions: S.S. and J.K. designed the experiments. A.M.-L. wrote the automated tracking program. J.K. and J.S. performed the experiments and analyzed the data. W.E.U. performed the theoretical and numerical calculations. J.K., W.E.U., and S.S. wrote the manuscript. All the authors discussed the results and commented on the manuscript. Competing interests: All authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.
View Abstract

Navigate This Article