## Abstract

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.

## INTRODUCTION

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 (*6*–*11*), and local fuel gradients, seeking fuel-rich regions over depleted ones (*12*–*14*). 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 (*19*–*22*).

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 (*25*–*28*). For instance, in “gyrotaxis,” bottom-heavy microorganisms, such as algae, adopt a stable, steady orientation and swim against the direction of gravity (*29*–*33*). Elongated microswimmers, such as sperms, display “rheotaxis,” the ability to orient and swim against the direction of flow (*27*, *34*–*38*). 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.

## RESULTS

### 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 H_{2}O_{2} solution, the Pt cap catalyzes the degradation of H_{2}O_{2}, 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 (*42*–*45*).

When the silica-Pt particles are suspended in water, they are quickly sedimented to the bottom surface because they are density-mismatched ( = 2.196 g/cm^{3}). 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/cm^{3}). Addition of H_{2}O_{2} 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 *V*_{p} = 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.

### 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 , 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 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 , 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 .

### Active particles in external flow

We then start a flow of H_{2}O_{2} 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 *V*_{p} 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 H_{2}O_{2} 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 *V*_{p} 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 *V*_{p} 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 *V*_{p} of the particles. For particles with *V*_{p} = 6 μm/s, the probability distribution function of α has a peak around −64°, whereas for a particle with *V*_{p} = 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**/*V*_{p}). Substituting the experimental values for *V*_{p} 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.

### 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 . 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 , 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 . In addition, the particle is carried downstream with velocity . 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 plane(1)

The component *p*_{y} is a constant determined by the initial orientation of the particle. The radius of the circular orbit is , and the speed of the tip 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 and **p**. Accordingly, translation of the particle is restricted to this plane, and **p** is restricted to remain within this plane (that is, ), although **p** can rotate toward or away from the wall (). We write , where we have defined a new, “primed” frame that has and in the plane of mirror symmetry, with . This frame is convenient for calculations because Ω^{(ax)} is strictly in (see Fig. 1B). The function incorporates the effect of bottom heaviness and interactions with the wall (for example, hydrodynamic interactions) that originate in swimming activity. Both 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 and in the following form(2)(3)(4)(5)Note that the axisymmetric contributions are evaluated in the primed frame, which is defined to corotate with the particle, and we used . The components of the translational velocity **U** in the *x* and *y* directions are determined by **p** and *h* as(6)(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, and (*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 () 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 is(8)which can be either downstream (), in agreement with our experimental observations, or upstream (), depending on the sign of . Because of the mirror symmetry of Eqs. 2 to 5 with respect to *p*_{y} = 0, this fixed point always occurs in pairs . 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 are in the direction. For most orientations θ and φ, will have a component in the vorticity direction . However, from Eq. 1, we see that the shear contribution never has a component. As a way out of this dilemma, we see from the definition of that, when θ = 90° (that is, *p*_{z} = 0), the and components of θ both vanish, and is strictly in the direction. Likewise, when θ = 90°, is strictly in . This suggests the possibility that all contributions can cancel out for θ* = 90° and some unknown φ*. Now, we consider the role of φ. The axisymmetric contributions have no dependence on φ, as discussed previously. The contribution from shear does depend on φ, as can be seen by substituting *p*_{z} = 0 and *p*_{x} = 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 (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 () 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 *V*_{p} = 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 as a small perturbation away from the fixed point: . We obtain the linearized governing equations , where the Jacobian matrix 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 . 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)(9)

The motion of the orientation vector resembles a damped harmonic oscillator with intrinsic frequency . Recalling that inactive, bottom-heavy particles were observed to continuously rotate in the experiments, we consider the case in which . We note that(10)where we have used and . 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 . The fixed point is a stable attractor if(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*, *48*–*50*). 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 (*42*–*45*). 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 and initial orientation θ_{0} = 30° and φ_{0} = 315°. The particle rotates so that its inert face points largely in the direction but with a small downstream orientation (*p*_{x} > 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.

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 *p*_{z} (black line) for the inactive particle in Fig. 5A. In Fig. 6F, we show *p*_{z} 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.

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 is shown as the vector field in Fig. 6D. Because this vector field is small but nonzero on the equator *p*_{z} = 0, adding it to the portrait in Fig. 6A will shift the center of oscillation (cyan circle) slightly toward the axis *p*_{x} = 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 *Pe*_{p} = 500 (see Materials and Methods for the definition) and shear rate . 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 , 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.

## DISCUSSION

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.

## MATERIALS AND METHODS

### 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.

### Tracking

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 by(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 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 **x**_{p} = (*x*_{p}, *y*_{p}, *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 ∇^{2}*c* = 0 with boundary conditions on the cap, on the inert region of the particle surface, and on the wall. Here, *D* is the diffusion coefficient of the solute molecule, **x** is a location in the fluid, and 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* + η∇^{2}**u** = 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** = **u**_{ext} + **U**^{(f)} + Ω^{(f)} × (**x** − **x**_{p}) on the particle. Here, **u**_{ext} is the external flow velocity, . 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, **v**_{s} = −*b*(**x**_{s}′)∇_{||}*c*(**x**′), where , and **x**_{s}′ is a location on a surface (in the primed frame). The “surface mobility” *b*(**x**_{s}′) encapsulates the details of the molecular interaction between the solute and the bounding surfaces. We wrote the boundary conditions **u** = **U**^{(s)} + Ω^{(s)} × (**x**′ − **x**_{p}′) + **v**_{s}(**x**_{s}′) on the particle and **u** = **v**_{s}(**x**_{s}′) 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, *b*_{inert} and *b*_{cap}, with *b*_{inert}/*b*_{cap} = 0.3 and *b*_{cap} < 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 *U*_{fs} = 13/80 *U*_{0}, where *U*_{0} = |*b*_{cap}|κ/*D* (*52*). The wall was characterized by a surface mobility *b*_{w}. We chose *b*_{w}/*b*_{cap} = −0.35. These surface mobility ratios were chosen to be similar to those used in a previous work [where we had *b*_{inert}/*b*_{cap} = 0.3 and *b*_{w}/*b*_{cap} = −0.2 (*11*)], and to give a slightly downstream steady orientation. We nondimensionalized length with *R*, velocity with *U*_{0}, and time with *T*_{0} = *R*/*U*_{0}. To nondimensionalize the gravitational and shear contributions, we must estimate *U*_{0} in real, dimensional units. Rather than to calculate *U*_{0} directly, which requires estimates for κ and |*b*_{cap}|, we used the expression *U*_{fs} = 13/80 *U*_{0} (*52*). Knowing that, experimentally, the particle characteristically moves at *U*_{fs} ≈ 5 μm/s when it is far from surfaces, we obtained *U*_{0} ≈ 30 μm/s. The parameters describing the heaviness of the particle are given in the Supplementary Materials. Finally, we considered the shear rate . 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 . 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.1*R* so that *h* ~ 1.1*R*. For *h*/*R* ≥ 1.05, the factor *g*(*h*/*R*) ≈ 1. Therefore, we estimated and a typical dimensionless shear rate to be . 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 ≡ *U*_{fs}*R*/*D* and . At room temperature, the diffusion coefficient of oxygen was *D* ~ 4 × 10^{−9}m^{2}/s (*53*) so that Pe ≈ 0.003 and . 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 = ρ_{fluid}*U*_{0}*R*/η, where η is the dynamic viscosity of the solution. Using η ~ 10^{−3}Pa s for water, we obtained Re ≈ 10^{−4}.

### Effect of thermal noise

The particle Peclet number characterizes the relative strengths of self-propulsion and translational diffusion of the particle. Here, is the translational diffusion coefficient of the particle in free space. For a particle with *R* = 2.5 μm in water at room temperature (*k*_{B}*T*_{room} ~ 4 × 10^{−21}J) with *U*_{fs} ≈ 5 μm/s, we estimated Pe_{p} ≈ 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 (Eqs. 2 to 4). Further details are provided in the Supplementary Materials and by Lisicki *et al*. (*55*).

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/4/1/eaao1755/DC1

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 *V*_{p} *≈* 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 *V*_{p} ≈ 6 μm/s and *V** ≈ 24 μm/s.

fig. S8. Fluctuations in φ obtained from five different particles for two values of *V*_{p} 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 *p*_{x} = |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.

## REFERENCES AND NOTES

**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.

- Copyright © 2018 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).