Cross-stream migration of active particles

TL;DR

This study investigates how spherical active particles, specifically catalytic Janus particles, align their propulsion axes perpendicular to flow and surfaces in channel flow, combining experiments and theoretical modeling.

Contribution

The paper introduces a deterministic model explaining the directional response of spherical microswimmers near surfaces under shear flow, supported by experimental validation.

Findings

Particles align nearly perpendicular to flow and surfaces.

Theoretical model captures experimental orientation distributions.

Thermal noise effects improve model-experiment agreement.

Abstract

For natural microswimmers, the interplay of swimming activity and external flow can promote robust motion, e.g. 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. Here, 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 semi-quantitatively agree with the experimental distributions.

Equations

$${\dot{\mathbf{p}}}^{(f)}\equiv {\mathbf{\Omega }}^{(f)}\times \mathbf{p}=\frac{1}{2}\dot{\gamma }f(h/R)\left[p_z\hat{x}-p_x\hat{z}\right]$$

$${\dot{p}}_x$$

$${\dot{p}}_y$$

$${\dot{p}}_z$$

$${\dot{h}}_{}$$

$$U_x$$

$$U_y$$

$$p_x^{\ast }=\cos ({\varphi }^{\ast })=\frac{2{\Omega }_{x^{\prime }}^{(ax)}(p_z^{\ast }=0,h^{\ast }/R)}{\dot{\gamma }f(h^{\ast }/R)},$$

$$\ddot{\tilde{p}}_{z}=-\left(\frac{1}{2}\dot{\gamma}f(H/R)p^{*}_{y}\right)^{2}\tilde{p}_{z}+\frac{\partial\Omega_{x^{\prime}}^{{\color[rgb]{0,0,0}(ax)}}}{\partial p_{z}}\biggr{\rvert}_{p_{z}^{*}=0,H}\dot{\tilde{p}}_{z}.$$

$$\frac{\partial\Omega_{x^{\prime}}^{(g)}}{\partial p_{z}}\biggr{\rvert}_{p_{z}^{*}=0,H}{\color[rgb]{0,0,0}=\frac{\partial\dot{\theta}^{(g)}}{\partial\theta}\biggr{\rvert}_{\theta^{*}=90^{\circ},H}}=0,$$

$$\frac{\partial\Omega_{x^{\prime}}^{(s)}}{\partial p_{z}}\biggr{\rvert}_{p_{z}^{*}=0,H}{\color[rgb]{0,0,0}=\frac{\partial\dot{\theta}^{(s)}}{\partial\theta}\biggr{\rvert}_{\theta^{*}=90^{\circ},H}}<0.$$

$$\text{out}(x,y)=\left\{\begin{array}[]{cc}1&\text{if img}(x,y)\geq T\\ 0&\text{if img}(x,y)<T\end{array}\right.$$

Full text

I Introduction

In micro-organisms 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. While the mechanism of some of these responses, such as chemotaxis, is active, involving the ability to sense gradients and signals, others are passive, purely resulting from external forces and torques. Gravitaxis in Paramecium, for example, is thought to arise solely due to the organism’s fore-rear asymmetry Roberts and Deacon (2002). Recently, there has been substantial effort to develop artificial micro-swimmers which mimic their natural counterparts in many ways Sánchez et al. (2015). They transduce energy from their local surroundings to engage in self-propulsion and display persistent random-walk trajectories, reminiscent of bacterial run and tumble behaviour Howse et al. (2007). They also respond to gravitational fields Campbell and Ebbens (2013); ten Hagen et al. (2014), physical obstacles Volpe et al. (2011); Spangolie and Lauga (2012); Takagi et al. (2016); Uspal et al. (2015a); Das et al. (2015); Simmchen et al. (2016), and local fuel gradients, seeking fuel-rich regions over depleted ones Hong et al. (2007); Baraban et al. (2013); Peng et al. (2015). A number of applications have been envisioned for these bio-mimetic swimmers, ranging from targeted drug delivery Patra et al. (2013); Gao and Wang (2014a) to environmental remediation Soler and Sánchez (2014); Gao and Wang (2014b). In many of these applications, it is inevitable that the artificial microswimmers will encounter flow fields while operating in confined spaces. While their swimming behaviour in various complex, static environments has been studied in detail, little is known about their response when they are additionally exposed to flow fields Tao and Kapral (2010); Zöttl and Stark (2012); Uspal et al. (2015b); Zöttl and Stark (2016).

In anticipation of this response, one can look to decades of previous work on biological micro-organisms in flow (see the review papers Pedley and Kessler (1992); Guasto et al. (2012) 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 Marcos et al. (2012); Chengala et al. (2013); Bukatin et al. (2015); Sokolov and Aranson (2016). For instance, in “gyrotaxis,” bottom-heavy micro-organisms like algae adopt a stable, steady orientation and swim against the direction of gravity Kessler (1985); Pedley and Kessler (1987); Durham et al. (2009); Thorn and Bearon (2010); Durham et al. (2011). Elongated microswimmers such as sperms display “rheotaxis”, the ability to orient and swim against the direction of flow Bretherton and Rothschild (1961); Hill et al. (2007); Kantsler et al. (2014); Figueroa-Morales et al. (2015); Bukatin et al. (2015); Omori and Ishikawa (2016). However, much of the work with biological microswimmers has focused on swimming in the bulk, and the detailed role of confining boundaries (e.g., in rheotaxis Kantsler et al. (2014); Uspal et al. (2015b); Palacci et al. (2015); Bukatin et al. (2015); Omori and Ishikawa (2016)) is an emerging area of research.

In this work, 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 flow only rotates as it is carried downstream. However, near-surface swimming activity adds an effective “rotational friction” to the dynamics of the particle, since it tends to drive the particle orientation into the plane of the wall. The effective “friction” damps the particle rotation and stabilises 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 flow near a surface. Our findings exemplify the complex behaviour that can emerge for individual microswimmers from the interplay of confinement and external fields. Moreover, they show that the qualitative character of the emergent behaviour is sensitive to the details of the interactions (e.g., hydrodynamic interactions) between individual microswimmers and bounding surfaces.

II Results

Active particles without external flow

In our experiments, we use silica colloids (1 $\mu m$ radius, or 2.5 $\mu 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 catalyses the degradation of $H_{2}O_{2}$ while 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. Additionally, particle-generated concentration gradients can induce chemi-osmotic slip on a nearby bounding surface, giving an additional contribution to particle motility. The details of this mechanism, shown schematically in Fig. 1(b), are comprehensively discussed elsewhere Anderson (1989); Golestanian et al. (2007), and are the subject of ongoing research Brown and Poon (2014); Ebbens et al. (2014); Brown et al. (2016, 2017).

When the silica-Pt particles are suspended in water, they quickly sediment to the bottom surface as they are density mismatched ($\rho_{SiO_{2}}$ = 2.196 $g/cm^{3}$). They orient with their caps down ($\theta=0^{\circ}$, where $\theta$ is defined in Fig. 1(a)) due to the bottom heaviness induced by the Pt layer ($\rho_{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 ($\theta=90^{\circ}$, see Fig. S1). We have previously reported on the dynamics leading to this change in orientation of the active particles Simmchen et al. (2016). Briefly, we could show that while the hydrodynamic interactions and the bottom heaviness of the particles tend to drive the particles towards the surface, the wall induced asymmetry of the distribution of the chemical product and the chemi-osmotic flow along the substrate (see Fig. 1(b)) tend to have the opposite effect, leading to a stable orientation at $\theta\approx 90^{\circ}$. 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\;{\mu m/s}$. Due to the contrast between the dark Pt hemisphere and the transparent silica, we can measure the angular orientations of these particles (see methods and Fig. S2 for details of the tracking process). Within the 2D plane, the particles have no preferred directionality and are diffusive on long time scales.

Passive particles in external flow

Now we seek to characterise the behaviour 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 non-slipping capillary surface, the flow velocity varies linearly as $v_{flow}=\dot{\gamma}z$ and the particles which are sedimented near the surface experience a shear flow (see Section S3 for a calculation of the flow profile) Bruus (2008). In terms of their translational behaviour, we observe that the particles act as “tracers” and translate in nearly straight lines along the direction of flow (Fig. 2(b)). 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 Palacci et al. (2015). Before we start the flow, the particles are all in the cap down orientation ($\theta=0^{\circ}$). The shear flow induces a torque on the particles and they rotate around the axis of flow vorticity $\hat{y}$ as they translate in flow (See SM Video 1). The rotation speed of the particles is also dependent on the flow rate (Fig. 2(a)), with higher flow rates leading to faster rotation. Via a simple model for particle rolling developed in Section S4 of the SM, we predict that the rotational period $\tau=2\pi/\sqrt{aV^{*,2}+b}$, 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 $\dot{\gamma}\approx 1.22\,V^{*}/R$.

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 behaviour in flow is drastically influenced in the presence of activity. Firstly, the particles stop rolling and reach a stable orientation parallel to the bottom surface ($\theta=90^{\circ}$). This is similar to previously observed behaviour for the particles in the absence of flow Simmchen et al. (2016). More surprisingly, the particles also evolve to a stable orientation that is nearly perpendicular to the direction of imposed flow ($\phi\approx 90^{\circ}$ or $-90^{\circ}$) (Fig. 2(e)). 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, i.e, 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. 2(f). Further, we observe that the stability of cross-stream migration (due to the stability of the steady orientation angle $\phi^{*}$ 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. 2(g) shows the angular evolution of two self-propelled particles in imposed flows of $V^{*}=14\,\mu m/s$ and $V^{*}=24\,\mu m/s$. The deviations away from the $\phi\approx 90^{\circ}$ positions occur more frequently and at larger amplitude than for particles in lower flow rates.

In order 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 $\phi$ in the system. These are plotted in Fig. 2(g) 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 $\phi$, and at long time scales, the particle behaviour is purely diffusive. However, in the case of imposed flow, we observe distinct peaks that appear at $\phi\approx 90^{\circ}$ and $-90^{\circ}$. These correspond to the particles exhibiting the cross-stream behaviour. These peaks also become sharper when we increase the flow rate, as can be seen in the distributions for $V^{*}=14\,\mu m/s$ and $V^{*}=24\,\mu 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. Since the particles are subject to Brownian fluctuations, they can occasionally also orient with or against the flow ($\phi\approx 0^{\circ}$ or $180^{\circ}$). This leads to an intermittent state where the particles “tumble” in the direction of flow before recovering their $\phi\approx\pm 90^{\circ}$ orientation and the cross-streamline behaviour. The recovered orientation of these particles can be different from their initial states, changing the direction of particle migration (e.g., from $+y$ to $-y$) (see Fig. S5).

Apart from the imposed flow rate $V^{*}$, we find that the stability of the $\phi^{*}$ is also dependent on the particle radius. Using particles of larger radius ($2.5\;{\mu m}$), we show that the distribution of $\phi$ is narrower around $\phi=\pm 90^{\circ}$ as compared to the distribution for $R=1\;{\mu m}$ particles for identical $V_{p}$ and $V^{*}$ (Fig. S6). Since the effect of Brownian noise is significantly lower on 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 behaviour of particle migration by tuning the self-propulsion velocity $V_{p}$ of the particles. Firstly, we find that higher propulsion velocities dampen the fluctuations around $\phi^{*}$ due to higher activity (See Fig. S8). Secondly, the propulsion velocity also controls the “slope” of the cross-stream migration. In order to quantify this, we define $\alpha$ to be the offset between the orientation vector $\mathbf{p}$ and the tracked velocity vector $\mathbf{v}$ of the particle. For a self-propelled particle in the absence of flow, the particles translate in the direction of the orientation vector $\mathbf{p}$ as the particles propel away from their Pt caps (i.e., $\alpha=0^{\circ}$), as shown in Fig. 3(a and b). However, with an imposed flow, the translational direction differs from the orientation vector, and the offset $\alpha$, for a given flow rate, is determined by the $V_{p}$ of the particles. For particles with $V_{p}=6\,\mu m/s$ the probability distribution function of $\alpha$ has a peak around $-64^{\circ}$, whereas for a particle with $V_{p}=3\,\mu m/s$, $\alpha$ is peaked at $-83^{\circ}$ (Fig. 3(a, c and d)). The offset $\alpha$ can be rationalized as a contribution of the transverse propulsion induced by the cross-stream orientation and the longitudinal advection by flow. The $\alpha$ should then simply be given by $\alpha=\arctan(V^{*}/V_{p})$. Substituting the experimental values for $V_{p}$ and $V^{*}$, we get $\alpha=66.9^{\circ}$ and $\alpha=78.3^{\circ}$, which agree reasonably well with the observed values. This clear separation in peaks of $\alpha$ 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 Uspal et al. (2015b). 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 (e.g., self-diffusiophoresis, self-electrophoresis, or mechanical propulsion by motion of surface cilia); it is, in that sense, generic. Accordingly, our major theoretical findings are, first, that there is a physical mechanism that can produce the surprising transverse orientational order observed in experiments, and, secondly, 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 the Methods section, the flow in the suspending fluid is characterised by a low Reynolds number, and hence governed by the Stokes equations. Since these equations are linear, the contributions of external flow, gravity, and swimming to the particle translational and angular velocities can be calculated independently and superposed: $\mathbf{U}=\mathbf{U}^{(f)}+\mathbf{U}^{(s)}+\mathbf{U}^{(g)}$ and $\mathbf{\Omega}=\mathbf{\Omega}^{(f)}+\mathbf{\Omega}^{(s)}+\mathbf{\Omega}^{(g)}$. The velocity of the orientation vector is determined by $\dot{\mathbf{p}}=\mathbf{\Omega}\times\mathbf{p}$. We will generally describe $\mathbf{p}$ in Cartesian coordinates, but will sometimes find it useful to use the spherical coordinates $\mathbf{p}=(\sin(\theta)\cos(\phi),\sin(\theta)\sin(\phi),\cos(\theta))$. Note that $|\mathbf{p}|=1$.

In order 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 $\hat{y}$, as shown in the left panel of Fig. 4. This is because the flow is faster near the upper surface of the particle (i.e., the surface farther away from the wall) than near the bottom surface. Accordingly, the particle has angular velocity $\mathbf{\Omega}^{(f)}=\frac{1}{2}\dot{\gamma}f(h/R)\hat{y}$. Additionally, the particle is carried downstream with velocity $\mathbf{U}^{(f)}=\dot{\gamma}hg(h/R)\hat{x}$. The functions $f(h/R)$ and $g(h/R)$ represent the influence of hydrodynamic friction from the wall Goldman et al. (1967); they are provided in SM Sect. S8. Due to the spinning motion of the particle, the tip of the vector $\mathbf{p}$ traces a circle in the shear plane:

[TABLE]

The component $p_{y}$ is a constant determined by the initial orientation of the particle. The radius of the circular orbit is $\sqrt{1-p_{y}^{2}}$, and the speed of the tip $|\dot{\mathbf{p}}^{(f)}|$ 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 geometry and surface activity profile; accordingly, we define $\mathbf{U}^{(ax)}\equiv\mathbf{U}^{(s)}+\mathbf{U}^{(g)}$ and $\mathbf{\Omega}^{(ax)}\equiv\mathbf{\Omega}^{(s)}+\mathbf{\Omega}^{(g)}$. This system, depicted in the middle panel of Fig. 4, has a plane of mirror symmetry defined by $\hat{z}$ and $\mathbf{p}$. Accordingly, translation of the particle is restricted to this plane, and $\mathbf{p}$ is restricted to remain within this plane (i.e., $\dot{\phi}=0$), although $\mathbf{p}$ can rotate towards or away from the wall ($\dot{\theta}\neq 0$). We write $\dot{\mathbf{p}}^{(ax)}=-\Omega_{x^{\prime}}^{(ax)}(\theta,h/R)\,\hat{\theta}$, where we have defined a new, “primed” frame that has $\hat{y}^{\prime}$ and $\hat{z}^{\prime}$ in the plane of mirror symmetry, with $\hat{z}^{\prime}=\hat{z}$. This frame is convenient for calculations because $\mathbf{\Omega}^{(ax)}$ is strictly in $\hat{x}^{\prime}$ (see Fig. 1(b)). The function $\Omega_{x^{\prime}}^{(ax)}(\theta,h/R)$ incorporates the effect of bottom-heaviness, as well as interactions with the wall (e.g., hydrodynamic interactions) that originate in swimming activity. Both $\Omega_{x^{\prime}}^{(ax)}(\theta,h/R)$ and $\mathbf{U}^{(ax)}(\theta,h/R)$ depend on the particle design and the model for the propulsion mechanism. In order to show that our analysis is generic in the sense discussed above, we will leave them unspecified until later in the text. Since the particle is axisymmetric and the wall is uniform, these functions have no $\phi$ dependence.

Now we superpose all contributions and obtain equations for $\dot{\mathbf{p}}=\mathbf{\Omega}\times\mathbf{p}$ and $\dot{h}=U_{z}$ in the following form:

[TABLE]

Note that the axisymmetric contributions are evaluated in the “primed” frame, which is defined to co-rotate with the particle, and we used $\hat{\theta}=\cos(\theta)\cos(\phi)\hat{x}+\cos(\theta)\sin(\phi)\hat{y}-\sin(\theta)\hat{z}$. The components of the translational velocity $\mathbf{U}$ in the $x$ and $y$ directions are determined by $\mathbf{p}$ and $h$ as:

[TABLE]

Steady state solutions

Eqs. 2-7 fully describe a deterministic active sphere in external shear flow near a planar surface. Now we look for the fixed points $(\mathbf{p}^{*},h^{*})$ of Eqs. 2-5. A fixed point is a particle configuration in which the particle translates along the wall with a steady height and orientation, i.e., $\dot{\mathbf{p}}=0$ and $\dot{h}=0$ Uspal et al. (2015b). 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 ($p_{z}^{*}=0$), and has non-zero components in both the flow and vorticity directions. These criteria are not satisfied by the other two fixed points. Slight misalignment from the vorticity axis is therefore a key experimental observation that discriminates between the steady states predicted by theory. For planar alignment, the streamwise component $p_{x}^{*}$ is

[TABLE]

which can be either downstream ($p_{x}^{*}>0$), in agreement with our experimental observations, or upstream ($p_{x}^{*}<0$), depending on the sign of \Omega_{x^{\prime}}^{{\color[rgb]{0,0,0}(ax)}}(p_{z}^{*}=0,h^{*}). Due to the mirror symmetry of Eqs. 2-5 with respect to $p_{y}=0$, this fixed point always occurs in pairs $(p_{x},\pm p_{y}^{*},0,h^{*})$. Further mathematical details are provided in the SM.

In order 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 $\mathbf{p}^{*}$. We recall that the axisymmetric contributions to $\dot{\mathbf{p}}$ are in the $\hat{\theta}$ direction. For most orientations $\theta$ and $\phi$, $\hat{\theta}$ will have a component in the vorticity direction $\hat{y}$. However, from Eq. 1, we see that the shear contribution never has a $\hat{y}$ component. As a way out of this dilemma, we see from the definition of $\hat{\theta}$ that when $\theta=90^{\circ}$ (i.e., $p_{z}=0$), the $\hat{x}$ and $\hat{y}$ components of $\hat{\theta}$ both vanish, and $\dot{\mathbf{p}}^{(ax)}$ is strictly in the $\hat{z}$ direction. Likewise, when $\theta=90^{\circ}$, $\dot{\mathbf{p}}^{(f)}$ is strictly in $\hat{z}$. This suggests the possibility that all contributions can cancel out for $\theta^{*}=90^{\circ}$ and some unknown $\phi^{*}$. Now we consider the role of $\phi$. The axisymmetric contributions $\dot{\mathbf{p}}^{(ax)}$ have no dependence on $\phi$, as discussed previously. The contribution from shear does depend on $\phi$, as can be seen by substituting $p_{z}=0$ and $p_{x}=\cos(\phi)$ into Eq. 1. Therefore, the sign and magnitude of the shear contribution can be “tuned” via the angle $\phi$ to cancel out the axisymmetric contributions to $\dot{\mathbf{p}}$ (provided these are not too large). Finally, we consider the height of the particle. Since shear does not contribute to vertical motion, the particle must obtain a steady height ($\dot{h}=0$) 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, Eq. 8 predicts that the steady state orientation approaches the vorticity axis as the flow rate is increased, i.e., $\cos(\phi^{*})\sim 1/V^{*}$. Accordingly, we perform further experiments with $5\;{\mu 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\;{\mu m/s}$ to $V^{*}=75\;{\mu m/s}$, for particles with $V_{p}=6\;{\mu m/s}$. In Fig. S11, we plot $\cos|\phi^{*}|$ against $V^{*}$ and find that the data recovers 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. Since, experimentally, the active particles spontaneously adopt a steady cross-stream orientation, the fixed point should be a stable attractor for active particles. Secondly, since 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 a certain condition on the particle/wall interaction is satisfied.

We define the generalized configuration vector $\mathbf{f}\equiv(\mathbf{p},h/R)$, and $\tilde{\mathbf{f}}$ as a small perturbation away from the fixed point: $\mathbf{f}=\mathbf{f}^{*}+\tilde{\mathbf{f}}$. We obtain the linearized governing equations $\dot{\tilde{{\mathbf{f}}}}=\mathbf{J}\tilde{\mathbf{f}}$, where the Jacobian matrix J_{ij}=\frac{\partial\dot{f}_{i}}{\partial f_{j}}\biggr{\rvert}_{{\color[rgb]{0,0,0}\mathbf{f}^{*}}} is given in detail in the SM. As a simplifying approximation, we take the particle to have a constant height $h=H$, so that $\dot{\tilde{h}}=0$. 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 SM):

[TABLE]

The motion of the orientation vector resembles a damped harmonic oscillator with intrinsic frequency $\omega_{0}=\frac{1}{2}\dot{\gamma}f(H/R)p^{*}_{x}$. Recalling that inactive, bottom-heavy particle were observed to continuously rotate in the experiments, we consider the case in which \Omega_{x^{\prime}}^{{\color[rgb]{0,0,0}(ax)}}=\Omega_{x^{\prime}}^{(g)}. We note that

[TABLE]

where we have used $\dot{\theta}^{(ax)}=-\Omega_{x^{\prime}}^{(ax)}$ and $\dot{\theta}^{(g)}\sim\sin(\theta)$. Therefore, the dissipative term in Eq. 9 is zero, and the orientation $\mathbf{p}$ has a continuous family of closed orbits centered on $\mathbf{p}^{*}$. Now we consider active, bottom-heavy particles, for which \Omega_{x^{\prime}}^{{\color[rgb]{0,0,0}(ax)}}=\Omega_{x^{\prime}}^{(s)}+\Omega_{x^{\prime}}^{(g)}. The fixed point is a stable attractor if

[TABLE]

This condition has the following interpretation: if the particle orientation is perturbed out of the plane of the wall, the active contribution to rotation responds (increases or decreases) so as to oppose the perturbation. If this condition is satisfied, activity induces an effective “friction” that damps oscillations of the orientation $\mathbf{p}$, driving attraction of $\mathbf{p}$ to $\mathbf{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, e.g., the comparison of pushers and pullers in Ref. Spangolie and Lauga (2012).)

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-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-5 by using a simple, well-established model of neutral self-diffusiophoresis in confinement Golestanian et al. (2005, 2007); Uspal et al. (2015a); Simmchen et al. (2016); Mozaffari et al. (2016); Ibrahim and Liverpool (2016). In the model, the particle has a hemispherical catalytic cap (Fig. 1(b), black), and the orientation vector $\mathbf{p}$ points from its catalytic pole to its inert pole. The particle emits solute molecules (i.e., 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. 1(b), 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 (e.g., ionic or electrokinetic) specific to the detailed self-phoretic mechanism Brown and Poon (2014); Ebbens et al. (2014); Brown et al. (2016, 2017). Further details are provided in Methods.

We consider some illustrative examples, using dimensionless parameters comparable to those in the experiments (see Methods for definition and estimation of parameters). We fix the particle height as $H/R=1.2$. In Fig. 5(a), we show a particle trajectory in a shear flow with strength $\dot{\gamma}R/U_{0}=0.1$ and initial orientation $\theta_{0}=30^{\circ}$ and $\phi_{0}=315^{\circ}$. The particle rotates so that its inert face points largely in the $-\hat{y}$ 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. 5(b)), rotating as it does so.

Since $h=H$ is taken to be constant, the instantaneous value of $\mathbf{p}$ completely determines the instantaneous velocity of the particle. In Fig. 6(a)-(c), phase space trajectories are shown on the unit sphere $|\mathbf{p}|=1$. The initial orientation in Fig. 5(a) and (b) is indicated by magenta circles in Fig. 6(a) and Fig. 6(b), respectively. We see that an inactive particle (Fig. 6(a)) has a continuous family of closed orbits in $\mathbf{p}$. When the particle is active, these oscillations are damped, and trajectories in phase space are attracted to the fixed point (Fig. 6(b)-(c)). In (c), 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. 6(e), we plot the component $p_{z}$ (black line) for the inactive particle in Fig. 5(a). In Fig. 6(f), we show $p_{z}$ for the active particle in Fig. 5(b). We can also obtain analytical solutions to the linearized equations, characterized by a few input parameters that are evaluated numerically (see SM). These solutions are shown as dashed red lines in Figs. 6(e) and (f), and agree well with the numerical data. Additionally, we note that if the particle height is allowed to change, one can still obtain the phenomenology studied here of attraction to a planar alignment steady state; an example is given in SI Sect. S8.

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. 6(a). How does adding activity transform this portrait into that shown in Fig. 6(b)? The contribution of activity to $\dot{\mathbf{p}}$ is shown as the vector field in Fig. 6(d). Since this vector field is small but non-zero on the equator $p_{z}=0$, adding it to the portrait in Fig. 6(a) will shift the center of oscillation (cyan circle) slightly towards the axis $p_{x}=0$, producing the green circle in Fig. 6(b). More significantly, the addition of this vector field will destabilize continuous oscillatory motion in the following way. Consider the closed trajectories in Fig. 6(a) that are closest to the cyan circle. In the neighborhood of the equator, activity always “pushes” the vector $\mathbf{p}$ towards 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 Methods). In Fig. 5(c), we show the probability distribution functions for the $\phi$ obtained for a Janus particle driven by shear flow at dimensionless inverse temperature $Pe_{p}=500$ (see Methods for definition) and shear rate $\dot{\gamma}R/U_{0}=0.1$. The distribution is symmetric, and has two peaks near the steady angle \phi^{*}=\pm{\color[rgb]{0,0,0}80.9}^{\circ} predicted by the deterministic model (Fig. 5(c)). For a higher shear rate $\dot{\gamma}R/U_{0}=0.5$, we find that the peaks in $\phi$ are sharper (Fig. 5(d)), in qualitative agreement with the experiments, with the peaks shifted near the deterministic prediction for this shear rate, $\phi^{*}=\pm{88.2}^{\circ}$. Additional probability distributions for the components of $\mathbf{p}$ are given in Fig. S14.

III Discussion

In this paper, 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 behaviour 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 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 (e.g., 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 (e.g. hydrodynamic interactions) can drive a rich directional response of spherical particles to external flows. This is in contrast with previous works on natural microswimmers (e.g., bacteria) where complex body shapes and flagellar beat patterns were implicated in directional response Marcos et al. (2012); Bukatin et al. (2015). Additionally, we have obtained semi-quantitative 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.

IV 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\,\mu 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\,\mu m$ silica particles (Sigma Aldrich) were prepared by drop casting of 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\,\mu 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 is first cleaned of noise by using blurring techniques, which substitute each pixel with an average of its surroundings. The particles are 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 results in a binary image out(x,y) given by:

[TABLE]

The gradient method is used for images with irregular brightness or when the particles are hard to distinguish from the background. In the first step of this method the gradient of the image $\nabla img=(\frac{\partial img}{\partial x},\frac{\partial img}{\partial y})$, is approximated by convolving the original frame with a Sobel operator. This results in two images, one which is the derivative along the X axis and another 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 is approximated as the center of mass of the contours obtained after segmentation. Particle trajectories are calculated using Bayesian decision-making, linking every particle center with the previous closest one. Intermediate missing positions, if any, are interpolated using cubic or linear splines.

In order to calculate the $\mathbf{p}$ vector, a line of predefined radius (approximately equal to the particle radius) is drawn through the center at a test angle. The standard deviation of the pixel values along the test-line is calculated and stored. The process is repeated multiple times and the standard deviations are compared. The test-angle with the least standard deviation is 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 present our model for calculation of the particle velocities $\mathbf{U}$ and $\mathbf{\Omega}$ as a function of $h$ and $\mathbf{p}$. We take the instantaneous position of the particle to be $\mathbf{x}_{p}=(x_{p},y_{p},h)$ in a stationary reference frame. The catalytic cap emits a product molecule at a constant and uniform rate $\kappa$. We take the solute number density $c(\mathbf{x})$ to be quasi-static, i.e., it obeys the equation $\nabla^{2}c=0$ with boundary conditions $-D\hat{n}\cdot\nabla c=\kappa$ on the cap, $\hat{n}\cdot\nabla c=0$ on the inert region of the particle surface, and $\hat{n}\cdot\nabla c=0$ on the wall. Here, $D$ is the diffusion coefficient of the solute molecule, $\mathbf{x}$ is a location in the fluid, $\hat{n}$ and is the normal vector pointing from a surface into the liquid. For each instantaneous configuration ($h$, $\mathbf{p})$, this set of equations can be solved for $c(\mathbf{x})$, e.g., numerically by using the boundary element method Pozrikidis (2002); Uspal et al. (2015b).

We take the velocity $\mathbf{u}(\mathbf{x})$ in the fluid solution to obey the Stokes equation $-\nabla P+\eta\nabla^{2}\mathbf{u}=0$, where $\eta$ is the solution viscosity and $P(\mathbf{x})$ is the pressure in the solution. Additionally, the fluid is incompressible, so that $\nabla\cdot\mathbf{u}=0$. The Stokes equation is a linear equation. Therefore, the contributions to $\mathbf{U}$ and $\mathbf{\Omega}$ from various boundary conditions for $\mathbf{u}$, as well as from various external forces and torques, can be calculated individually as the solution to separate subproblems and then superposed. Hence, we write \mathbf{U}=\mathbf{U}^{(f)}+\mathbf{U}^{{\color[rgb]{0,0,0}(ax)}} and \mathbf{\Omega}=\mathbf{\Omega}^{(f)}+\mathbf{\Omega}^{{\color[rgb]{0,0,0}(ax)}}, where $(f)$ indicates the contributions of the external flow, and the axisymmetric contributions are from gravity and swimming activity: \mathbf{U}^{{\color[rgb]{0,0,0}(ax)}}=\mathbf{U}^{(g)}+\mathbf{U}^{{\color[rgb]{0,0,0}(s)}} and \mathbf{\Omega}^{{\color[rgb]{0,0,0}(ax)}}=\mathbf{\Omega}^{(g)}+\mathbf{\Omega}^{{\color[rgb]{0,0,0}(s)}}. For each subproblem, the fluid is governed by the Stokes equation and incompressibility condition.

We first consider the subproblem for the contribution of the external flow. The fluid velocity is subject to no-slip boundary conditions $\mathbf{u}=0$ on the planar wall and $\mathbf{u}=\mathbf{u}_{ext}+\mathbf{U}^{(f)}+\mathbf{\Omega}^{(f)}\times(\mathbf{x}-\mathbf{x}_{p})$ on the particle. Here, $\mathbf{u}_{ext}$ is the external flow velocity, $\mathbf{u}_{ext}=\dot{\gamma}z\hat{x}$. Additionally, the particle is free of external forces and torques, closing the system of equations for $\mathbf{U}^{(f)}$ and $\mathbf{\Omega}^{(f)}$. The solution of this subproblem is well-known; see, for instance, Goldman et al. Goldman et al. (1967)

The two axisymmetric subproblems are calculated in the “primed” frame co-rotating with the particle (Fig. 1(b)). For the subproblem associated with particle activity, we employ the classical framework of neutral self-diffusiophoresis Anderson (1989); Golestanian et al. (2005, 2007); Uspal et al. (2015a); Simmchen et al. (2016). In this subproblem, the self-generated solute gradients drive surface flows on the wall and the particle surface, $\mathbf{v}_{s}=-b(\mathbf{x}_{s}^{\prime})\nabla_{||}c(\mathbf{x}^{\prime})$, where $\nabla_{||}\equiv(\mathbf{1}-\hat{n}\hat{n})\cdot\nabla$ and $\mathbf{x}_{s}^{\prime}$ is a location on a surface (in the primed frame). The “surface mobility” $b(\mathbf{x}_{s}^{\prime})$ encapsulates the details of the molecular interaction between the solute and the bounding surfaces. We write the boundary conditions \mathbf{u}=\mathbf{U}^{{\color[rgb]{0,0,0}(s)}}+\mathbf{\Omega}^{{\color[rgb]{0,0,0}(s)}}\times(\mathbf{x}^{\prime}-\mathbf{x}_{p}^{\prime})+\mathbf{v}_{s}(\mathbf{x}_{s}^{\prime}) on the particle, and $\mathbf{u}=\mathbf{v}_{s}(\mathbf{x}_{s}^{\prime})$ on the wall. Again, specifying that the particle is force and torque free closes the system of equations for \mathbf{U}^{{\color[rgb]{0,0,0}(s)}} and \mathbf{\Omega}^{{\color[rgb]{0,0,0}(s)}}. This subproblem can be solved numerically using the boundary element method Pozrikidis (2002); Uspal et al. (2015b). For the subproblem associated with gravity, we use the “eggshell” model of Campbell and Ebbens 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 Campbell and Ebbens (2013). Details concerning this subproblem, which is solved using standard methods Simmchen et al. (2016), are provided in Supplementary Materials.

We now specify the parameters characterizing the system. We choose 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 is far away from bounding surfaces, moves in the $\mathbf{p}$ direction (i.e., away from its cap) with a velocity $U_{fs}=13/80\;U_{0}$, where $U_{0}=|b_{cap}|\kappa/D$ Michelin and Lauga (2014). The wall is characterized by a surface mobility $b_{w}$. We choose $b_{w}/b_{cap}=-0.35$. These surface mobility ratios are chosen to be similar to those used in previous work (Ref. Simmchen et al. (2016), where we had $b_{inert}/b_{cap}=0.3$ and $b_{w}/b_{cap}=-0.2$), and to give a slightly downstream steady orientation. We non-dimensionalize length with $R$, velocity with $U_{0}$, and time with $T_{0}=R/U_{0}$. In order to non-dimensionalize the gravitational and shear contributions, we must estimate $U_{0}$ in real, dimensional units. Rather than calculate $U_{0}$ directly, which requires estimates for $\kappa$ and $|b_{cap}|$, we use the expression $U_{fs}=13/80\;U_{0}$ Michelin and Lauga (2014). Knowing that, experimentally, the particle characteristically moves at $U_{fs}\approx 5\;{\mu m/s}$ when it is far from surfaces, we obtain $U_{0}\approx 30\;{\mu m/s}$. The parameters describing the heaviness of the particle are given in Supplementary Information. Finally, we consider the shear rate $\dot{\gamma}$. This is not known experimentally, but can be roughly estimated by considering inactive particles to act as passive tracers. From Goldman et al., the velocity of a spherical particle driven by shear flow near a wall is $U^{(f)}_{x}=g(h/R)\dot{\gamma}h$ Goldman et al. (1967). Experimentally, inactive particles in flow are observed to move with $V^{*}\sim 10\,{\mu m/s}$. The particle height $h$ is difficult to observe experimentally, but we take it to be set by the balance of gravity and electrostatic forces. Hence, the particle/wall gap $\delta$ is on the order of a Debye length $\lambda_{D}\sim 0.1R$, so that $h\sim 1.1R$. For $h/R\geq 1.05$, the factor $g(h/R)\approx 1$. We therefore estimate $\dot{\gamma}\approx 10\,{s^{-1}}$, and a typical dimensionless shear rate to be $\dot{\gamma}R/U_{0}\approx 0.8$. As a reminder, our aim is to establish semi-quantitative agreement with experiments, and therefore we seek only 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\equiv U_{fs}R/D$ and $Pe_{\dot{\gamma}}\equiv\dot{\gamma}R^{2}/D$. At room temperature, the diffusion coefficient of oxygen is $D\sim 4\times 10^{-9}m^{2}/s$ Popescu et al. (2010), so that $Pe\approx 0.003$ and $Pe_{\dot{\gamma}}\approx 0.015$. Furthermore, in taking the fluid velocity $\mathbf{u}$ to be governed by the Stokes equation, we neglected fluid inertia. This approximation is justified for low Reynolds number, $Re=\rho_{fluid}U_{0}R/\eta$, where $\eta$ is the dynamic viscosity of the solution. Using $\eta\sim 10^{-3}Pa\,s$ for water, we obtain $Re\approx 10^{-4}$.

Effect of thermal noise

The particle Peclet number $Pe_{p}=U_{fs}R/D_{0}^{t}$ characterizes the relative strengths of self-propulsion and translational diffusion of the particle. Here, $D_{0}^{t}=k_{B}T/6\pi\eta R$ is the translational diffusion coefficient of the particle in free space. For a particle with $R=2.5\;{\mu m}$ in water at room temperature ($k_{B}T_{room}\sim 4\times 10^{-21}J$) with U_{fs}\approx 5\;{\mu m{\color[rgb]{0,0,0}/s}}, we estimate $Pe_{p}\approx 880$. In order to perform Brownian dynamics simulations, we adapt the Euler-Maruyama integration scheme introduced by Jones and Alavi Jones and Alavi (1992), and later presented by Lisicki et al. Lisicki et al. (2014), which explicitly includes the effects of the wall on diffusion. Our principal modification to the method is inclusion of deterministic contributions to $\dot{\mathbf{p}}$ (Eqs. 2-4). Further details are provided in Supplementary Information and Ref. Lisicki et al. (2014).

Acknowledgements.

The authors thank M. N. Popescu for helpful discussions. W.E.U. acknowledges financial support from the DFG, 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 MINECO 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 311529.