Open access peer-reviewed chapter

The Boundary Element Method for Fluctuating Active Colloids

Written By

William E. Uspal

Submitted: October 3rd, 2018 Reviewed: May 8th, 2019 Published: June 23rd, 2019

DOI: 10.5772/intechopen.86738

Chapter metrics overview

1,005 Chapter Downloads

View Full Metrics


The boundary element method (BEM) is a computational method particularly suited to solution of linear partial differential equations (PDEs), including the Laplace and Stokes equations, in complex geometries. The PDEs are formulated as boundary integral equations over bounding surfaces, which can be discretized for numerical solution. This manuscript reviews application of the BEM for simulation of the dynamics of “active” colloids that can self-propel through liquid solution. We introduce basic concepts and model equations for both catalytically active colloids and the “squirmer” model of a ciliated biological microswimmer. We review the foundations of the BEM for both the Laplace and Stokes equations, including the application to confined geometries, and the extension of the method to include thermal fluctuations of the colloid. Finally, we discuss recent and potential applications to research problems concerning active colloids. The aim of this review is to facilitate development and adoption of boundary element models that capture the interplay of deterministic and stochastic effects in the dynamics of active colloids.


  • active colloids
  • Brownian dynamics
  • boundary element method

1. Introduction

Over the past 15 years, significant effort has been invested in the development of synthetic micro- and nano-sized colloids capable of self-propulsion in liquid solution [1, 2, 3]. These “active colloids” have myriad potential applications in drug delivery [4, 5], sensing [6], microsurgery [7], and programmable materials assembly [8]. Furthermore, they provide well-controlled model systems for study of materials systems maintained out of thermal equilibrium by continuous dissipation of free energy. In this context, and in comparison with driven systems (e.g., sheared suspensions), a unique aspect of active colloids is that energy is injected into the system at the microscopic scale of a single particle, instead of through macroscopic external fields or at the boundaries of the system. As a consequence of this, novel collective behaviors are possible, including motility-induced phase separation [9], mesoscopic “active turbulence” [10], and formation of dynamic “living crystals” and clusters [11, 12]. Furthermore, since living systems can be regarded as self-organized non-equilibrium materials systems, study of active colloids could yield insight into fundamental principles of living systems, and open a path towards development of biomimetic “dissipative materials” capable of homeostasis [13], self-repair [14], goal-directed behavior [15, 16], and other aspects of life.

Paradigmatic examples of synthetic active colloids include bimetallic Janus rods [17] and Janus spheres consisting of a spherical core with a hemispherical coating of a catalytic material [18]. In both cases, self-propulsion is driven by catalytic decomposition of a chemical “fuel” available in the liquid solution. For instance, for gold/platinum Janus rods, both ends of the rod are involved in the electrochemical decomposition of hydrogen peroxide into water and oxygen: hydrogen peroxide is oxidized at the platinum anode and reduced at the gold cathode. In this reaction process, a hydrogen ion gradient is established between the anode and cathode. The resulting gradient in electrical charge creates an electric field in the vicinity of the rod. The electric field exerts a force on the diffuse layer of ions surrounding the colloid surface, resulting in motion of the suspending fluid relative to the colloid surface. Viewed in a stationary reference frame, the final result is “self-electrophoretic” motion of the colloid in direction of the platinum end. For Janus spheres (e.g., platinum on silica or platinum on polystyrene), the mechanism of motion is still a subject of debate. Since the core material is inert and insulating, it was originally thought that these particles move by neutral self-diffusiophoresis in a self-generated oxygen gradient. Diffusiophoresis is similar to electrophoresis in that motion is driven by interfacial molecular forces. Briefly, in diffusiophoresis, the colloid surface and solute molecules interact through some molecular potential. This interaction potential, in conjunction with a gradient of solute concentration along the surface of the colloid, leads to the pressure gradient in a thin film surrounding the colloid, and therefore fluid flow within the film relative to the colloid surface. Following initial studies on chemically active Janus spheres, subsequent studies revealed a dependence of the Janus particle speed on the concentration of added salt [19], suggesting that a self-electrophoretic mechanism may be implicated in motion of the colloid. Golestanian and co-workers proposed that dependence of the rate of catalysis on thickness of the deposited catalyst can lead to different regions of the catalyst acting as anode and cathode [20]. More recently, it was proposed that if one of the redox reactions is reaction-limited and the other is diffusion-limited, the anodic or cathodic character of a point on the catalytic surface will depend on the local curvature of the surface [21]. Regardless of the detailed molecular mechanism of motion, a key point is that interfacial flows drive self-propulsion of chemically active colloids. A second key point is that particles need to have an intrinsic asymmetry (e.g., from the Janus character of their material composition) in order to exhibit directed motion.

These findings have motivated development of theoretical and numerical concepts for modeling the interfacially driven self-propulsion of active colloids. Motivated by classical work on phoresis in thermodynamic gradients [22, 23], an influential continuum framework for modeling neutral self-diffusiophoresis was established in Ref. 24, and will be reviewed below. This basic framework can be modified or extended to account for electrochemical effects [25], multicomponent diffusion [26], reactions in the bulk solution [27], and confinement [28, 29, 30, 31, 32, 33, 34]. An emerging area of study within this framework is autonomous navigation and “taxis” of chemically active colloids in ambient fields and complex geometries, including chemotaxis in chemical gradients [35] and rheotaxis in confined flows [15, 36]. Theoretical research on synthetic active colloids has also found common ground with an older strand of research on locomotion of biological microswimmers. Here, an important point of contact is again the idea of interfacial flow [37]. For a quasi-spherical microswimmer that is “carpeted” with a layer of cilia, the effect of the periodic, time-dependent, metachronal motion of the cilia can be modeled as a period-averaged interfacial flow. This “squirmer” model of locomotion was introduced by Lighthill [38] and refined by Blake [39]. More recent work has explored collective motion of suspensions of squirmers [40] and squirmer motion in confined geometries [41].

These theoretical frameworks are deterministic, and do not directly address the role of thermal fluctuations. For instance, for the model of a chemically active colloid in Ref. 42, diffusion of the chemical reaction product (i.e., the solute) into the surrounding solution is modeled with the Laplace equation, which has a smooth and unique solution for a given set of boundary conditions describing surface catalysis. Implicit in the use of the Laplace equation are the assumptions that, on the timescale of Janus particle motion, the solute diffuses very fast, and that fluctuations of the solute distribution average out to be negligible. Likewise, fluctuations of the surrounding fluid are neglected, i.e., the deterministic Stokes equation is used to model the fluid in lieu of the fluctuating Stokes equation. On the other hand, micron-sized active Janus particles are observed in experiments to exhibit “enhanced diffusion”: directed motion on short timescales t < τ r and random walk behavior on long timescales t τ r . For the latter, the effective diffusion coefficient D eff is enhanced relative to the “bare” diffusion coefficient D 0 of an inactive colloid, i.e., D eff D 0 . The reason for this behavior is that the orientation of the particle is free to fluctuate, and the particle changes its direction of motion by rotational diffusion over the timescale τ r = D r 1 , where D r is the rotational diffusion coefficient of the particle [18]. Therefore, thermal fluctuations qualitatively affect the motion of even a micron-sized catalytic Janus particle in unbounded, uniform solution. For a catalytic Janus particle in an ambient field or in confinement, thermal fluctuations affect whether and for how long the particle can align with the ambient field [42, 43] or stay near confining surfaces [34, 44]. Overall, a full theoretical understanding of the behavior of micron-sized active colloids requires modeling thermal fluctuations.

Moreover, as part of the general drive towards miniaturization, recent experimental efforts have sought to fabricate and characterize nano-sized chemically active colloids [45, 46, 47]. On the theoretical side, new questions arise when the size of the colloid becomes comparable to the size of the various molecules participating in the catalytic reaction. These questions include: When is using a continuum model appropriate [48]? Can a catalytic particle still display (time- and ensemble-averaged) directed motion when the particle and the surrounding chemical field are fluctuating on similar timescales? Relatedly, can a spherical colloid with a catalytic surface of uniform composition exhibit enhanced diffusion when nano-sized [49]? Can a fluctuating, nano-sized Janus particle effectively follow an ambient chemical gradient, i.e., exhibit chemotaxis [35]? These questions also connect with the burgeoning literature on chemotaxis of biological enzymes [50].

In this chapter, we review the boundary element approach to modeling the motion of active colloids. This is a “hydrodynamic” approach that resolves the detailed geometry and surface chemistry of the colloids, the velocity of the surrounding solution, and the distribution of chemical species within the solution [30, 40, 51, 52, 53, 54, 55, 56, 57]. The advantage of such an approach—in comparison with, for instance, the active Brownian particle model—is that it can resolve the detailed microscopic physics of how a colloid couples to ambient fields and other features of the surrounding micro-environment. In addition, we discuss how thermal fluctuations can be included within the approach. The aim of this review is to facilitate development and adoption of models that capture the interplay of deterministic and stochastic effects within an integrated framework.


2. Theory

As a starting point, we review the basic deterministic theoretical framework for understanding the motion of active colloids [24]. This is a continuum approach that coarse-grains the interfacial flow that drives colloid motion, discussed above, as a “slip velocity” boundary condition for the velocity of the suspending fluid.

We consider a suspension of N active colloids in an unbounded liquid solution. The position of each colloid α , with α 1 2 N , is described in a stationary reference frame by a vector x α . The solution is modeled as an incompressible Newtonian fluid with dynamic viscosity μ . The solution is governed by the Stokes equation,

P + μ 2 u = 0 , E1

where P x is the pressure at a position x in the solution, and u x is the velocity of the solution. The velocity obeys the incompressibility condition,

u = 0 , E2

and the boundary condition

u x s = U α + Ω α × x s x α + v s , α x s , x s S α , E3

where S α is the surface of colloid α , x s S α is a position on S α , and U α and Ω α are the translational and rotational velocities, respectively, of colloid α . The quantity v s , α x s is the slip velocity on the surface of colloid α , which is either prescribed (for a squirmer) or determined by the distribution of chemical species in solution (for a chemically active colloid). The form of v s , α x s for the two types of particles will be discussed in detail below. Additionally, far away from the N particles, the fluid velocity vanishes:

u x = 0 . E4

In order to close this system of equations, we require 6 N more equations, corresponding to the 6 N unknown components of U α and Ω α . The net force and torque on each colloid vanishes:

S α σ n ̂ dS + F ext , α = 0 , E5
S α x s x α × σ n ̂ dS + T ext , α = 0 , E6

where the integrals are performed over the surface S α of each colloid α , and F ext , α and T ext , α are, respectively, the net external force and net external torque on the colloid. The stress tensor is given by

σ = P I + μ u + u T , E7

where the pressure P x is determined by the incompressibility condition.

Practitioners of Stokesian Dynamics may notice some similarity between Eq. 3 and the boundary condition for an inert or passive sphere in an ambient flow field. If v s , α x s could be expressed as an effective ambient flow field at the position of particle α , the tools of Stokesian dynamics could be straightforwardly applied to simulation of active suspensions. This analogy will be developed in the Appendix.

2.1 The squirmer model: prescribed surface slip

The “squirmer” model was originally introduced by Lighthill to describe the time-averaged motion of ciliated quasi-spherical micro-organisms [38]. Lighthill’s formulation was subsequently corrected and extended by Blake [39]. The basic motivating idea of the squirmer model is that the periodic, metachronal motion of the carpet of cilia on the surface of the micro-organism drives, over the course of one period and in the vicinity of the microswimmer surface, net flow from the “forward” or “leading” pole of the micro-organism to the “rear” pole (see Figure 1, left). This interfacial flow drives flow in the surrounding bulk fluid, leading to directed motion of the micro-organism towards the forward end. The squirmer model captures some essential features of the self-propulsion of micro-organisms, including the hydrodynamic interactions between micro-organisms, and between an individual micro-organism and confining surfaces.

Figure 1.

The two types of microswimmer considered in this chapter. In the spherical “squirmer” model (left), the slip velocity on the surface of the particle is specified by fiat and fixed (in a frame co-moving and co-rotating with the particle) for all time. For an axisymmetric distribution of surface slip, the particle moves in the direction of the green arrow, i.e., opposite to the (surface-averaged) direction of the slip velocity. For a Janus particle (right), a fraction of the particle surface (black) catalyzes a reaction involving various molecular species diffusing in the surrounding solution. The resulting anisotropic distribution of product molecules or “solute” (green spheres) drives a phoretic slip velocity (purple) in an interfacial layer surrounding the particle. For a repulsive interaction between the solute and the particle surface, the slip is towards high concentration of solute, and the particle moves in the direction of the green arrow.

The slip velocity on the surface of a spherical squirmer α is specified by fiat and does not depend on the configuration of the suspension. It is given as [41]:

v s , α x s = n B n , α V n cos θ e ̂ θ p , α , E8


V n x 2 1 x 2 n n + 1 dP n x dx . E9

The unit vector e ̂ θ p , α is defined in following manner (see Figure 1, left). The squirmer has an axis of symmetry, and a propulsion direction d ̂ α oriented along this axis. We define the body-frame polar angle θ p , α at a point x s on the surface of squirmer α as the angle between d ̂ α and a vector x s x α from the center of the squirmer x α to x s . The unit vector e ̂ θ p , α is oriented in the direction of increasing θ p , α , i.e., locally tangent to the squirmer surface along a longitudinal line.

The squirming mode amplitudes B n , α , which can potentially vary from squirmer to squirmer, are fixed a priori and do not depend on the configuration of the suspension. The set of amplitudes determine the detailed form of the flow field in vicinity of the particle. Furthermore, the lowest order squirming mode B 1 determines the velocity of an isolated squirmer in unbounded solution: U fs , α = 2 / 3 B 1 , α . According to our definition of d α and θ p , α , we require that B 1 , α > 0 . Simulations of squirmers typically truncate Eq. 8 to n 2 or n 3 . The justification for this is that the contributions of the higher order squirming modes to the flow around the squirmer decay rapidly with distance from the squirmer.

2.2 Chemically active colloids: diffusiophoretic slip from chemical gradients

For chemically active colloids, the slip velocity on the surface of a colloid is driven by interfacial molecular forces. The molecular physics of phoresis and self-phoresis is reviewed in detail elsewhere [2, 23, 58]; here, we provide a brief summary. Consider a “Janus” colloid with a surface composed of two different materials. In the presence of molecular “fuel” diffusing in the surrounding solution, one of the two Janus particle materials catalyzes the decomposition of the fuel into molecular reaction products. A paradigmatic example of this reaction is the decomposition of hydrogen peroxide by platinum into water and oxygen:

H 2 O 2 Pt H 2 O + 1 2 O 2 . E10

(This equation is a severe simplification of the actual reaction scheme, which most likely involves charged and complex intermediates [20, 27]; nevertheless, proceeding from it, we can capture some essential features of self-phoresis.) If the reaction is reaction-limited—i.e., hydrogen peroxide is plentifully available in solution, and diffuses quickly relatively to the reaction rate—then we can approximate the production of oxygen with zero order kinetics:

D c n ̂ x = x s = κ x s , E11

where D is the diffusion coefficient of oxygen, c x is the number density of oxygen, and κ x s is the rate of oxygen production on the surface of the particle. (The validity of assumption of reaction-limited kinetics is quantified by the Damköhler number Da = κ 0 R / D , where κ 0 is a characteristic reaction rate; we assume Da 1 .) Furthermore, we assume that the Péclet number Pe U 0 R / D is very small, where U 0 is a characteristic particle velocity and R is the particle radius. Accordingly, we can make a quasi-steady approximation for the diffusion of oxygen in the solution:

2 c = 0 . E12

Finally, we assume that

c x = c , E13

where c is a constant. Eqs. 11, 12, and 13 specify a boundary value problem (BVP) for the distribution of oxygen in the fluid domain containing the N active particles. This problem can be solved numerically, e.g., by the boundary element method, as will be outlined in a later section.

Accordingly, each Janus particle will be surrounded by an anisotropic “cloud” of oxygen molecules (“solute”), with the oxygen concentration highest near the catalytic cap (see Figure 1, right). Now we suppose that the oxygen molecules interact with the surface of the colloid through some molecular interaction potential with range δ R [23]. Each colloid is surrounded by an interfacial layer of thickness δ in which the molecular interaction of the solute and the colloid is significant. Outside of this layer, the solute freely diffuses in the solution. We can regard c x as the bulk concentration, i.e., the concentration outside the interfacial layer. Near a location x s on the surface of the colloid, the interfacial layer concentration is enhanced or depleted, according to the attractive or repulsive character of the molecular interaction, relative to c x s + . Here, the plus sign emphasizes that c x s + is evaluated outside the interfacial layer. Moreover, since δ R , the interfacial layer concentration can locally, in the direction locally normal to the colloid surface, relax to a Boltzmann (i.e., equilibrium) distribution governed by the molecular interaction potential Φ . (The timescale for this local relaxation is much faster than the characteristic timescale for colloid motion R / U 0 .) Accordingly, the local pressure P x s η can be calculated from Φ and c x s + , where η is a coordinate defined at x s that is locally normal to the colloid surface.

These notions can be made mathematically rigorous through the theory of matched asymptotics. However, for the purpose of this discussion, the essential idea is that the bulk concentration c x determines the pressure in the interfacial layer in the vicinity of a point x s on the colloid surface. Moreover, c x varies over the length scale R of the colloid. Accordingly, within the interfacial layer, the pressure varies over the size of the colloid, driving flow within the interfacial layer. From the perspective of the outer solution for the flow field, this interfacial flow looks like a slip velocity:

v s , α x s = b x s c . E14

Here, the surface gradient operator is defined as I n ̂ n ̂ . The material-dependent parameter b x s encapsulates the details of the molecular interaction between the solute and the surface material, and can be calculated from the molecular potential Φ [23]. Since the surface of the Janus colloid comprises different materials, b depends on the location on the colloid surface. In fact, a spatial variation of b over the surface of colloid is a necessary condition to obtain phoretic rotation of a colloid near a wall [30] or chemotactic alignment with a gradient of “fuel” molecules [35].

2.3 Lorentz reciprocal theorem

The Lorentz reciprocal theorem relates the fluid stresses σ σ and velocity fields u u of two solutions to the Stokes equation within the same domain V :

S u σ n ̂ dS = S u σ n ̂ dS , E15

where S is the boundary of V . For the N active particles in unbounded solution, S = α = 1 N S α .

This theorem can be used to simplify the problem specified above for the velocities of N active particles. We designate that problem as the “unprimed” problem. Additionally, we specify that F ext , α = 0 and τ ext , α = 0 for all α . (Since the Stokes equation is linear, the contributions of the external forces and torques to the velocities of the particles can be calculated separately, using standard methods, and superposed with the contributions from activity to obtain the complete velocities.) We consider 6 N “primed” problems, indexed by j = 1 , 2 , , 6 N . For problem j = α , particle α is exposed to an external force with magnitude F ext in the x ̂ direction. Each particle has unknown translational and rotational velocities U α j and Ω α j , and the velocity field u j is subject to no-slip boundary conditions on each particle:

u j x s = U α j + Ω α j × x s x α , x s S α , E16

Additionally, the flow field vanishes far away from the particles, i.e., u j x = 0 . The unprimed problem and primed problem α are schematically illustrated in Figure 2. Similarly, for problems j = α + 1 and j = α + 2 , particle α is exposed to an external force with magnitude in the y ̂ and z ̂ directions, respectively, with no-slip boundary conditions likewise holding on each particle, and the flow field vanishing far away from the particles. For problems j = α + 3 , j = α + 4 , and j = α + 5 , particle α is exposed to a torque with magnitude τ ext in the x ̂ , y ̂ , and z ̂ directions, respectively, with the same boundary conditions. Each “primed” problem j can be solved for 6 N unknown velocity components U α j and Ω α j .

Figure 2.

Illustration of the “unprimed” problem for the velocities of N active particles, and the “primed” problem α for the velocities of N inert particles when particle α is exposed to a force with magnitude F ext in the x ̂ direction. Similarly, in primed problems α + 1 and α + 2 , particle α is exposed to a force with magnitude F ext in the y ̂ direction and the z ̂ direction, respectively. Moreover, in primed problems α + 3 , α + 4 , and α + 5 , particle α is exposed to torques with magnitude τ ext in the x ̂ , y ̂ , and z ̂ directions, respectively. Note that the primed and the 6 N unprimed problems all have the same geometry, i.e., the same configuration of N spheres.

For problem j , we substitute Eqs. 3 and 16 into Eq. 15 for u and u = u ' j to obtain:

α S α U α + Ω α × x s x α + v s , α x s σ j n ̂ dS = E17
α S α U α j + Ω α j × x s x α σ n ̂ dS . E18

It can be shown that the right hand side of this equation vanishes. Consider the term involving U α j . For each integral over S α , U α j is a constant and can be moved out of the integral,

U α j S α σ n ̂ dS , E19

but the integral is simply the force F α on particle α . Since the particles are free of external forces, F α = 0 . Likewise, the term involving Ω α j can be rearranged as

Ω α j S α x s x α × σ n ̂ dS , E20

but the integral is the torque τ α = 0 on particle α .

Rearranging the left hand side of Eq. 17, we obtain a set of 6 N equations j :

α U α F α j + Ω α τ α j = α S α v s σ j n ̂ dS , j = 1 , , 6 N . E21

These 6 N equations can be written in matrix form:

R V = b , E22

where b is a vector containing the 6 N integrals associated with the right hand side of Eq. 21, V is vector of all 6 N velocity components U α Ω α T , and R is the grand resistance matrix for N spheres at positions x α . Note that the arbitrary magnitudes F ext and τ ext have been divided out of Eq. 22.

The advantage of the reciprocal theorem approach is that if we solve the “primed” problem for a given set of particle positions x α , we can easily compute the set V of 6 N velocities for any set of slip velocities v s , α . This, for instance, facilitates studying how various choices for b x s or B n , α affect particle motion. Additionally, the “primed” problem for the resistance matrix R and stresses σ j in a system of N spheres is a standard problem in microhydrodynamics. An interesting open question is whether this approach is numerically more stable than directly solving for the 6 N particle velocities in the presence of the force-free and torque-free constraints.

2.3.1 Proof of Lorentz reciprocal theorem

We provide a short proof of Eq. 15, following the lines of Ref. 59 because some intermediate results will be useful later in the chapter. We recall that the rate of strain tensor e ij is defined as

e ij = 1 2 u i x j + u j x i , E23

and that, in index notation, the stress tensor is

σ ij = P δ ij + 2 μ e ij E24

We consider the quantity σ ij e ij :

σ ij e ij = P δ ij + 2 μ e ij ' e ij = P e ii + 2 μ e ij e ij = 2 μ e ij e ij , E25

where we have used u = 0 to eliminate e ii , and we assume the Einstein convention for summation over repeated indices. Similarly, we can obtain σ ij e ij = 2 μ e ij e ij , so that

σ ij e ij = σ ij e ij . E26

We can also manipulate σ ij e ij as follows:

σ ij e ij = 1 2 σ ij u i x j + 1 2 σ ij u j x i . E27

Swapping the two indices in the last term,

σ ij e ij = 1 2 σ ij u i x j + 1 2 σ ji u i x j . E28

But σ ij = σ ji , giving

σ ij e ij = σ ij u i x j = x j σ ij u i σ ij x j , E29

so that

x j σ ij u i σ ij x j u i = x j σ ij u i σ ij x j u i . E30

If there are no point forces applied to the fluid in determination of u and u , then σ = 0 and σ = 0 , and we obtain

u σ = u σ . E31

Integrating both sides over the volume V and applying the divergence theorem, we obtain Eq. 15.

2.4 Boundary integral formulation of the Laplace equation

Even with the aid of the Lorentz reciprocal theorem, it is necessary to solve the Stokes and (for self-phoretic particles) Laplace equations in a fluid domain containing the active particles as interior boundaries. For most configurations of the suspension, an analytical solution is intractable, and a numerical approach is required. Many numerical methods (e.g., the Finite Element Method) discretize and solve the governing partial differential equations in the three-dimensional fluid domain. This can be computationally intensive. Moreover, if the domain is unbounded in one or more dimensions, the computational domain must be truncated. Typically, the computational domain must be large in order to accurately approximate an unbounded solution, and significant computational effort must be expended on calculating the flow, pressure, and concentration fields far away from the particles.

An alternative approach proceeds from the following insight: a linear boundary value problems can be reformulated as a boundary integral equation (BIE) on the domain boundaries [51, 60]. Furthermore, the boundary integral equation can be discretized for numerical solution, yielding a dense linear system of coupled boundary element equations in the form of A q = b . One significant advantage of the boundary element method is that it requires discretization of only bounding surfaces; for instance, to represent an unbounded system of N particles, one need only mesh the surfaces of the N spheres. As a disadvantage, the coefficient matrix A is typically fully populated and non-symmetric; therefore, for a system of N elm elements, the required computer memory scales as O N elm 2 , and the required computation time scales as O N elm 3 .

In order to obtain the BIE for the Laplace equation, we first consider the divergence theorem:

V A dV = S A n ̂ dS , E32

where the volume integral on the left hand side is carried out over the entire solution domain V , and the surface integral on the right hand side is carried out over all boundaries S . We include a negative sign on the right hand side of the equation because we define n ̂ to point into the solution domain (see Figure 3). If A = ϕ ψ , where ϕ x and ψ x are scalar fields, then the divergence term A = ϕ 2 ψ + ϕ ψ , and we obtain Green’s first identity:

Figure 3.

Schematic illustration of the geometry for development of the boundary integral equations for the Laplace and stokes equations. The fluid domain is denoted by V , the solid domain by V p , and the interior of particle α by V p , α , where V p = α = 1 N V p , α . The solid and fluid domains are separated by the particle surfaces S α , with S = α = 1 N S α . The observation point x 0 can occur in V , in V , or on S ; we show x 0 in V .

V ϕ 2 ψ + ϕ ψ dV = S ϕ ψ n ̂ dS . E33

We can also write Green’s first identity for A = ψ ϕ :

V ψ 2 ϕ + ψ ϕ dV = S ψ ϕ n ̂ dS . E34

Subtracting Eq. 34 from Eq. 33, we obtain Green’s second identity:

V ϕ 2 ψ ψ 2 ϕ dV = S ϕ ψ ψ ϕ n ̂ dS . E35

Now, we let ϕ = c x , with 2 c = 0 . Furthermore, we let ψ = G x x 0 , where the Green’s function G x x 0 satisfies Poisson’s equation:

2 G x x 0 + δ x x 0 = 0 . E36

We obtain:

V c x 2 G x x 0 dV = S c x G x x 0 G x x 0 c x n ̂ dS . E37

We have not yet specified the location of the pole x 0 . If x 0 is located in the domain V , then, using the properties of the Dirac delta function, we obtain an integral representation of the concentration field c x 0 at a point x 0 V :

c x 0 = S c x G x x 0 G x x 0 c x n ̂ dS . E38

Using the divergence theorem, can show that 2 1 x x 0 = 4 πδ x x 0 , so that

G x x 0 = 1 4 π x x 0 . E39

We recall from electrostatics that G x x 0 represents the electrostatic potential at x from a point charge of unit strength located at x 0 . Within the context of the steady diffusion equation 2 c = 0 , it has a different physical interpretation: it can be regarded as the steady concentration c x at a point x due to a point-like, steady source of concentration, continuously injected into the system, located at x 0 and with unit strength (i.e., unit number density flux per unit time). One can take derivatives of G x x 0 with respect to x 0 to obtain higher order multipole singularities. For instance, we can obtain the Green’s function G dp x x 0 for a source/sink dipole located at x 0 as

G dp x x 0 x 0 G x x 0 = x x 0 4 π x x 0 3 . E40

As G dp x x 0 is a vector quantity, we obtain a scalar contribution to c x by multiplying with a dipole vector d ; the magnitude and direction of d specify the strength and orientation of the dipole.

By inspection, the Green’s function obeys the symmetry property G x x 0 = G x 0 x , so we can rewrite Eq. 38 as

c x 0 = S c x G x 0 x G x 0 x c x n ̂ dS . E41

Interestingly, we have obtained an expression for c x 0 in the solution domain in terms of the values of c x and c x n ̂ on the domain boundaries. Note that this is not a solution to a boundary value problem for c x 0 , since a BVP specifies only one of the quantities c x and c x n ̂ on the domain boundaries. Specifically, for the problem of a system of catalytic particles outlined above, we only know c x n ̂ a priori. Eq. 41 has an interesting physical interpretation: c x 0 can be regarded as the concentration due to a distribution of monopoles (i.e., point sources of mass flux) with strength c x n ̂ on the boundaries, plus a distribution of point dipoles (i.e., infinitesimally separated pairs of mass sources and sinks) with strength c x and orientation n ̂ on the boundaries.

We still have two other options for where to place x 0 : inside the boundary S (i.e., outside the solution domain V ) or somewhere on the boundary S . If we place x 0 inside S , then the integral on the left hand side of Eq. 37 vanishes:

S c x G x x 0 G x x 0 c x n ̂ dS = 0 . E42

Placing x 0 on the boundary requires some care in how to handle the Dirac delta function on the left hand side of Eq. 37. If we regard the Dirac delta as the limit of a sequence of distributions, then it is clear that a factor of one half should arise when we integrate over V :

1 2 c x 0 = S c x G x x 0 G x x 0 c x n ̂ dS . E43

This is a boundary integral equation (BIE) because the left hand side is the concentration c x 0 at a point on the boundary, while the right hand side is an integral of c x and c x n ̂ over the same boundary. A boundary value problem typically specifies one of these two quantities on the boundary; the other can be obtained with Eq. 43.

In the boundary element method, the boundary integral equation is discretized for numerical solution. Here, we briefly summarize the method, and direct the reader to consult the useful and comprehensive book of Pozrikidis for further information [51]. Each particle is represented as a meshed, closed surface. The meshing only needs to be done once; for a dynamical simulation, no remeshing during the simulation is required, even if the particles move relative to each other. The concentration c and its normal derivative c n ̂ are assumed to be uniform over element i . For a point x 0 on the surface, we obtain the boundary integral equation:

1 2 c x 0 = i = 1 N elm c i S i G x x 0 n ̂ dS c n ̂ i S i G x x 0 dS . E44

Choosing x 0 as the midpoint x j of element j , we can write N elm equations:

1 2 c j = i = 1 N elm c i S i G x x j n ̂ dS c n ̂ i S i G x x j dS . E45

The N elm equations can be written in matrix form:

A ij 1 2 δ ij c j = B ij c n ̂ j , E46


A ij S i G x x j n ̂ dS E47


B ij S i G x x j dS . E48

Given either a specification of either c j (Dirichlet boundary conditions) or c n ̂ j (Neumann boundary conditions), the algebraic system in Eq. 46 can be solved numerically with standard methods.

A certain difficulty becomes apparent when we consider the element B ii : the evaluation point x i lies within the element of integration, and therefore the integral contains the singularity of Eq. 39. We are saved from a potentially disastrous situation by the fact that the integral is carried out over an area. Nevertheless, this singular integral has to be handled with care. Further technical information, as well as a wealth of practical details concerning implementation of the BEM, is available in Ref. [51].

As a further note, issues with singular integrals have motivated development of regularized boundary element methods, which use a regularized Green’s function, i.e., a Green’s function with the singularity “smeared out” over a finite size ε [52, 54, 61].

2.5 Boundary integral formulation of the Stokes equation

A similar approach can be taken for the Stokes equation [51, 59]. Recall that the Stokes equation is:

σ = P + μ 2 u = 0 . E49

We can define a Green’s function G x x 0 as the solution u x G x x 0 F to the Stokes equation with a body point force F located at x 0 :

P + μ 2 u + F δ x x 0 = 0 , E50


σ = F δ x x 0 . E51

It can be shown that the Green’s function is

G ij x x 0 = 1 8 πμr δ ij + x ˜ i x ˜ j r 2 , E52

where r x x 0 and where x ˜ j x j x 0 , j . Eq. 52 is commonly called the Oseen tensor or “Stokeslet”. The fluid pressure in response to the point force is given by P x P x x 0 F , where

P j = x ˜ j 4 π r 3 . E53

The stress in the fluid is given by σ Σ x x 0 F , where

Σ ijk = 3 x ˜ i x ˜ j x ˜ k 4 π r 5 . E54

Now we wish to apply Eq. 30. We specify the “primed” fields u and σ as the fields due to a point force F at x 0 in unbounded fluid. For the “unprimed” fields u and σ , we specify that they are fields of interest in the domain V bounded by S (see Figure 3). Furthermore, V is free of body forces, so that σ = 0 . We obtain:

x j u i Σ ijk F k u i x j Σ ijk F k = x j σ ij G ik F k E55
F k x j u i Σ ijk + u i F i δ x x 0 = F k x j G ik σ ij E56

We integrate both sides over the domain V :

V u k F k δ x x 0 dV = F k x j G ik σ ij x j u i Σ ijk dV E57

Now we apply the divergence theorem:

F k V u k δ x x 0 dV = F k S G ik x x 0 σ ij u i Σ ijk x x 0 n j dS , E58

where the negative sign appears because of our convention that n ̂ points into V . We note that G ik x x 0 = G ik x 0 x and Σ ijk x x 0 = Σ ijk x 0 x . Additionally, the choice of F k was arbitrary. We can therefore write:

V u k δ x x 0 dV = S G ik x 0 x σ ij + u i Σ ijk ( x 0 x ) n j dS . E59

If we choose to place x 0 in V , we obtain a boundary integral representation for u k x 0 :

u k x 0 = S G ik x 0 x σ ij x + Σ ijk ( x 0 x ) u i x n j dS . E60

As with Eq. 41, the boundary integral representation for the flow field has an interesting physical interpretation. The first term on the right hand side of Eq. 60 can be regarded as a “single layer potential” due to a distribution of point forces with strength σ n ̂ over the surface of the particle. The second term on the right hand side is the “double layer potential.” Detailed examination of this term reveals that it can be decomposed into the superposition of the flow due to a distribution u n ̂ of point sources and sinks of fluid mass, plus the flow to a distribution of point force dipoles [59].

If we place x 0 outside V , i.e., inside the space V p enclosed by the particles, we obtain

S G ik x 0 x σ ij x + Σ ijk ( x 0 x ) u i x n j dS = 0 . E61

Finally, if we place x 0 on the boundary S , we obtain a boundary integral equation:

1 2 u k x 0 = S G ik x 0 x σ ij x + Σ ijk ( x 0 x ) u i x n j dS . E62

For rigid body motion, including the 6 N “primed” problems for the Lorentz reciprocal theorem, the double layer can be eliminated from the boundary integral equation as follows [59]. Consider extending the volume filled by “fluid” to V p . Within V p , the flow field u is simply the flow field u RBM for rigid body motion, since it must obey no-slip on the surface S . Now we apply Eq. 59 for the field u RBM inside V p and x 0 V , noting that we must use a normal n ̂ = n ̂ pointing into V p :

S G ik x 0 x σ ij RBM x + Σ ijk ( x 0 x ) u i RBM x n j dS = 0 . E63

For rigid body motion, there is no shear stress and the pressure is uniform, i.e., σ ij RBM x n ̂ = p 0 n ̂ . The first term is simply the integral of G ik x 0 x n ̂ over the surface S for x 0 V , which vanishes identically by incompressibility. This leaves:

S Σ ijk x 0 x u i RBM x n j dS = 0 . E64

Examining Eq. 62, we note that u ( x ), i.e., the flow velocity in V , is equal to u RBM on S . Therefore, we conclude:

u k x 0 = S G ik x 0 x σ ij x n j dS , x 0 V . E65

In order to obtain a single-layer boundary integral equation for x 0 S , note that the jump discontinuity responsible for the factor of 1 / 2 in Eq. 62 is strictly from the double-layer potential [59]. The contribution of the single layer potential to the velocity field is continuous across S . We obtain:

u k x 0 = S G ik x 0 x σ ij x n j dS , x 0 S . E66

This single layer boundary integral equation can be discretized and solved numerically in a similar manner as the Laplace equation; Ref. 51 provides a comprehensive account.

2.6 Active suspensions in confined geometries

In the preceding, we considered a suspension of N particles in an unbounded three-dimensional geometry. However, the presence of confining boundaries can have a significant effect on the dynamics of a suspension. It is possible to include a solid surface by explicitly meshing it and including it as a “fixed” or immobile particle in the calculations [53]. This approach is typically necessary for solid surfaces with corners or complex topography. One disadvantage of this approach is that an infinite surface (e.g., an infinite planar wall) must be truncated and included as a finite size object. Care must be taken that the mesh is sufficiently fine near the particles, so that, for instance, the concentration and flow fields do not “leak” through a solid wall, but also that the mesh is sufficiently coarse far away from the particles, so that computation time is tractable.

A second, “mesh-free” approach is suitable for confining geometries with high symmetry, such as an infinite planar wall [39], an interface between two fluids with different viscosities [62], a fluid domain bounded by a solid wall and a free interface [63], or even two infinite planar walls. Additionally, it can be suitable if the domain is periodic in two or three dimensions. In this approach, the Green’s functions for the Laplace and Stokes equations are replaced with Green’s functions that obey the desired boundary conditions on the bounding surfaces. The Green’s function in the confined geometry can often be constructed by the method of images.

2.7 Thermal fluctuations

So far, we have considered the deterministic contributions to the 6 N components of velocity for a suspension of N particles. However, as outlined in the Introduction, the interplay of these deterministic contributions and the stochastic Brownian forces on the particles is important—and in some problems, such as the long-time behavior of an active colloid, it is absolutely essential.

One approach to include Brownian forces on an active particle, the hybrid boundary element/Brownian dynamics method, simply calculates them separately and superposes them with the deterministic contributions. Using the Itô convention for stochastic differential equations, this superposition is expressed by the overdamped Langevin equation for the generalized, 6 N -component coordinate q :

d q dt = V + k B T M + 2 k B T B W , E67

where V is the deterministic contribution of activity to the generalized velocity U α Ω α T , i.e., the solution to the problem outlined above; M is the grand mobility matrix M = R 1 ; B satisfies B B T = M ; and W is a collection of independent Wiener processes. Discretizing time in steps of Δ t , one can write a generalized displacement Δ q as the following Euler-Maruyama equation [34, 64]:

Δ q = V Δ t + k B T M Δ t + 2 k B T B Δ w , E68

where Δ w is a stochastic variable with Δ w = 0 and Δ w i Δ w j = δ ij Δ t . The stochastic drift term M is a consequence of having a configuration-dependent mobility tensor in the framework of the Itô interpretation.

The update of the orientation of each particle α should respect the constraint that d α = 1 and avoid any errors arising from application of (non-commuting) rotation matrices in arbitrary order to d α . There are robust algorithms for rigid body motion that represent the orientations of the particles with quaternions [65], Euler angles [66], or rotation matrices that transform between body-fixed and global reference frames [67].

The stochastic drift term in Eq. 68 can present some difficulty for numerical calculations [66]. For some simple situations, such as a single spherical colloid near a planar wall [34, 42], solutions for the configuration dependence of the mobility tensor are available in the literature [68, 69]. Alternatively, Eq. 67 can be discretized and solved via Fixman’s midpoint method to avoid calculation of the drift term [70].

This approach assumes that that the colloid and the fluid are not fluctuating on the same timescale, i.e., the fluid velocity is integrated out as a fast variable. Additionally, for self-phoretic particles, this approach necessarily neglects fluctuations of the chemical field c x in the fluid domain V .

A recently developed variation of the boundary element method for Stokes flow, the fluctuating boundary element method, does not make this post hoc superposition of deterministic and Brownian contributions to particle motion. Rather, fluctuations are directly incorporated into boundary integral equation via a random velocity field on the boundary S [71].


3. Discussion and conclusions

The boundary element method is emerging as a powerful and important method for numerical simulation in the field of synthetic active colloids [30, 52, 54, 55, 56, 57]. This new area of application follows many years of fruitful application to modeling biological microswimmers, including with the squirmer model [40, 53]. For active colloids, a major advantage of the boundary element approach is that it can resolve the microscopic details of phoretic self-propulsion, including the chemical and flow fields generated by an active colloid, the surface chemistry and shape of the colloid, and the microscopic physics of how the colloid can couple to ambient fields and confining surfaces.

A few examples serve to illustrate the utility of the approach. Ref. [30] considers the dynamics of a spherical active Janus colloid near a planar wall. The colloid can “sense” and respond to the wall through self-generated chemical and hydrodynamic fields. Specifically, the wall provides a no-flux boundary condition for the solute concentration, and a no-slip boundary condition for the flow field. By confining the solute, the wall enriches the concentration of solute in the space between the particle and the wall, breaking the axial symmetry of the concentration field. Concerning the flow, the flow created by the particle scatters off the wall and back to the particle. These effects are captured by the boundary element method, including their dependence on the size of the catalytic cap and the spatial variation in the surface mobility b over the surface of the particle. As another example, Ref. [43] considers the dynamics of a photo-active spherical Janus colloid. The catalytic cap of the colloid is only active when exposed to incident light. This self-shadowing effect, in conjunction with the spatial variation of b on the surface of the colloid, leads to phototaxis (rotation of the cap towards the light) or anti-phototaxis (rotation of the cap away from the light.) Notably, this work uses the hybrid BEM/BD method to calculate the distribution of particle orientations as a function of illumination intensity and particle surface chemistry. Concerning the interaction of multiple particles, Ref. 57 uses the regularized BEM to calculate the dynamics of multiple isotropic spherical colloids. Interestingly, a group of N = 5 particles can form a stable cluster with broken rotational symmetry. This broken symmetry allows propulsion of the whole cluster. Finally, concerning shape, the BEM has been used to model toroidal [54] and spherocylindrical [72] self-phoretic particles.

However, some caveats are in order. For the hybrid boundary element/Brownian dynamics method discussed in this work, neither the fluctuations of the suspending fluid nor of the chemical field(s) are explicitly resolved. For self-phoretic particles in the ångstrom to nanometer size range, the particles, the solute, and the solvent fluctuate on similar timescales. Additionally, the validity of the continuum description of the surrounding solution is questionable. Molecular and mesoscopic simulation methods that resolve discrete solute and solvent particles may be more appropriate in this size range [48]. As a second caveat, boundary element methods are most suited to solution of linear governing PDEs, such as the Laplace and Stokes equations. Introducing nonlinearity in the governing equations (e.g., for a solution with nonlinear rheology or nonlinear bulk reaction kinetics) leads to the appearance of volume integrals in the boundary integral formulation. Thirdly and relatedly, the boundary element method is not as easily extensible as other methods (e.g., the finite element method) for inclusion of more complicated multiphysics. Finally, there is a caveat specific to active colloids. Much remains unknown about the reaction kinetics for self-phoretic particles. The boundary element method can have many free microscopic parameters (e.g., the values of the surface mobility b on different surfaces); this raises the danger of overfitting to experimental results.

As a potential direction of research, we suggest developing a hybrid computational method combining the advantages of BEM and Stokesian Dynamics (SD). Stokesian dynamics is a method for simulating the dynamics of colloidal suspensions [73, 74, 75, 76]. Far-field hydrodynamic interactions are included in SD, truncated at the level of the stresslet (i.e., the first moment of the stress on the surface of a particle, which produces a hydrodynamic disturbance decaying as 1 / r 2 .) Near-field hydrodynamic interactions are typically included via lubrication forces acting between particle pairs. Due to these approximations, Stokesian dynamics is computationally much cheaper than BEM, allowing access to collective dynamics, the rheology of dense suspensions, etc. On the other hand, SD does not typically resolve the microscopic details of individual particles, such as shape or heterogeneous surface chemistry. A hybrid BEM-SD method could combine the detailed microscopic resolution of BEM for near-field interactions with the ability of SD to capture many-body phenomena driven by far-field interactions. (This hybrid approach would bear some similarity to the fast multipole method.) In the Appendix at the end of this chapter, we develop a starting point for including interfacial flows v s x within the standard SD formalism for spherical particles.

As a second potential research direction, one could consider deformable active particles using the BEM. The boundary element method for Stokes flow has been coupled to methods to model particle elasticity, including the finite element method, in order to study the deformation of fluid-filled capsules [77] and elastic particles in shear flow [78], as well as the deformation of blood cells squeezing through constrictions [79].

The boundary element method could also be used to investigate questions touching upon fundamental nonequilibrium statistical mechanics. For instance, do nonequilibrium steady states of squirmers or self-phoretic particles (e.g., stable clusters of catalytic particles [57]) minimize the rate of entropy production [80]? When do hydrodynamic interactions suppress or enhance motility-induced phase separation and other nonequilibrium phase transitions? Does the pressure of an active suspension on a boundary obey an equation of state when hydrodynamic and phoretic interactions with the boundary are considered [76, 81]?

In any case, we anticipate that the boundary element method will continue to find successful application in the microswimmers field. A few potential problems include: modeling the collision dynamics and scattering of two or more non-spherical active colloids [72, 82]; the interaction of an active colloid and a passive colloid, possibly including the formation of dimeric bound states for cargo transport; and further exploration of motion near bounding surfaces and interfaces, especially fluid/fluid interfaces.


Consider an inert (non-active) sphere of radius R in an ambient flow field u x . The sphere has translational velocity U and rotational velocity Ω . The flow field u can be written as u x = u x + u D x , where u D x is the velocity disturbance created by the presence of the sphere. The boundary condition for u D x on the sphere surface S is

u D x s = U + Ω × x s x 0 u x s , x s S . E69

Additionally, u D x 0 . Taking the sphere position to be x 0 = 0 , we can expand the ambient flow field around the sphere center as:

u i x = u i 0 + u i x j x = 0 x j + 1 2 2 u i x j x k x = 0 x j x k + E70

Now we recall the definitions of the (symmetric) rate of strain tensor e ij ,

e ij = 1 2 u i x j + u j x i , E71

and the (anti-symmetric) vorticity tensor

W ij = 1 2 u i x j u j x i , E72

The vorticity tensor can be related to the vorticity vector w = × u by

W = 1 2 ε w . E73

Here, ε is the Levi-Civita tensor. The first derivative in Eq. 70 can be decomposed into symmetric and anti-symmetric contributions:

u i x j = e ij + W ij . E74

Using the Lorentz reciprocal theorem, one can obtain Faxén’s law for the drag force on the sphere (see Ref. [59] for details):

F drag = 6 πμR 1 + R 2 6 2 u x 0 U . E75

(In our shorthand notation, the Laplacian is first applied to u x and then evaluated at x 0 .) One can also obtain Faxén law for the drag torque:

τ drag = 8 πμ R 3 ω x 0 Ω , E76

where the angular velocity of the fluid ω 1 2 w . Finally, there is a Faxén law for the stresslet [59, 83, 84]

S = 20 3 πμ R 3 1 + R 2 10 2 e x 0 , E77

where S is defined as an integral over the particle surface

S ij = 1 2 x j σ ik n k + x i σ jk n k 1 3 x k σ kl n l δ ij μ u i n j + u j n i dS . E78

So far we have only presented standard results, but now we raise the following question. Consider an active sphere with a slip velocity v s x . Comparing the boundary conditions in Eq. 3 and Eq. 69, can we construct an ambient linear flow field

u x = u 0 + e x + ω × x E79

with u x s = v s x s on S ? Constructing an effective flow field would allow us to obtain F drag , τ drag , and S by Faxén laws, without having to solve the complete hydrodynamic problem posed in Section II. Moreover, an understanding of how v s x determines F drag , τ drag , and S would pave the way towards development of a hybrid BEM-Stokesian Dynamics scheme, since these quantities are central to SD.

As our starting point, we write the Taylor expansion of u x :

u i x = u i 0 + u i x j x = 0 x j . E80

To obtain u i 0 , we integrate both sides of Eq. 80 over the surface of the sphere:

u i x dS = u i 0 dS + u i x j x = 0 x j dS . E81

We identify u i x on the surface of the sphere as v s x s . The second integral on the right hand side of Eq. 81 vanishes, giving

u 0 = 1 4 π R 2 v s x s dS . E82

Using Eq. 75, we obtain

F drag = 6 πμR 1 4 π R 2 v s x s dS U . E83

If we consider a force-free swimmer, F drag = 0 , giving the result:

U = 1 4 π R 2 v s x s dS . E84

This equation is one of the major results obtained in Ref. 37 by use of the Lorentz reciprocal theorem. However, our rederivation and interpretation in terms of an effective ambient flow field u is (to our knowledge) novel. To obtain the vorticity associated with u , we multiply Eq. 80 by ε lmi x m and integrate over the sphere surface:

ε lmi x m u i x dS = u i 0 ε lmi x m dS + u i x j x = 0 ε lmi x m x j dS . E85

The first integral on the right hand side of Eq. 85 vanishes. For the second integral on the right hand side, we use the identity

x m x j dS = 4 π R 4 3 δ jm . E86

We obtain:

ε lmi x m u i x dS = 4 π R 4 3 ε lji u i x j x = 0 E87
x × v s x s dS = 4 π R 4 3 × u x = 0 = 8 π R 4 3 ω 0 , E88

so that

ω 0 = 3 8 π R 4 x × v s x s dS . E89

Using Eq. 76, we obtain:

τ drag = 8 πμ R 3 3 8 π R 4 x × v s x s dS Ω . E90

For a torque-free swimmer, τ drag = 0 , and we obtain a second major result from Ref. [37]:

Ω = 3 8 π R 4 x × v s x s dS . E91

Finally, we consider how to obtain the stresslet S . We multiply Eq. 80 by x m and integrate over the surface of the sphere:

u i x x m dS = u i 0 x m dS + u i x j x = 0 x m x j dS . E92

The first integral on the right hand vanishes, giving

u i x x m dS = 4 π R 4 3 u i x m x = 0 . E93

Swapping the indices i and m , we can also write:

u m x x i dS = 4 π R 4 3 u m x i x = 0 . E94

Adding these two equations and dividing by two, we obtain

1 2 u i x x m + u m x x i dS = 4 π R 4 3 e im 0 . E95


e im 0 = 3 8 π R 4 v s , i x x m + v s , m x x i dS . E96

Using the Faxén Law in Eq. 77, we obtain:

S = 5 μ 2 R v s x s x + x v s x s dS . E97

This is the major result obtained in Ref. 84 via the Lorentz reciprocal theorem. As before, this manuscript provides a novel alternative derivation and interpretation of Eq. 97 in terms of an effective ambient flow field. (Note that, due to the linearity of the Stokes equation, our approach is easily extended to model active particles in a real ambient flow field.)


  1. 1. Sánchez S, Soler L, Katuri J. Chemically powered micro- and nanomotors. Angewandte Chemie, International Edition. 2015;54:1414
  2. 2. Moran JL, Posner JD. Phoretic self-propulsion. Annual Review of Fluid Mechanics. 2016;49:511
  3. 3. Aubret A, Ramananarivo S, Palacci J. Eppur si muove, and yet it moves: Patchy (phoretic) swimmers. Current Opinion in Colloid & Interface Science. 2017;30:8189
  4. 4. Sundararajan S, Lammert PE, Zudans AW, Crespi VH, Sen A. Catalytic motors for transport of colloidal cargo. Nano Letters. 2008;8:1271
  5. 5. Baraban L, Tasinkevych M, Popescu MN, Sánchez S, Dietrich S, Schmidt OG. Transport of cargo by catalytic Janus micro-motors. Soft Matter. 2012;8:48
  6. 6. Wu J, Balasubramanian S, Kagan D, Manesh KM, Campuzano S, Wang J. Motion-based DNA detection using catalytic nanomotors. Nature Communications. 2010;1:36
  7. 7. Solovev AA, Xi W, Gracias DH, Harazim SM, Deneke C, Sanchez S, et al. Self-propelled nanotools. ACS Nano. 2012;6:1751
  8. 8. Aubret A, Youssef M, Sacanna S, Palacci J. Targeted assembly and synchronization of self-spinning microgears. Nature Physics. 2018;14:1114-1118
  9. 9. Cates ME, Tailleur J. Motility-induced phase separation. Annual Review of Condensed Matter Physics. 2015;6:219
  10. 10. Wensink HH, Dunkel J, Heidenreich S, Drescher K, Goldstein RE, Lowen H, et al. Meso-scale turbulence in living fluids. Proceedings of the National Academy of Sciences. 2012;109:14308
  11. 11. Palacci J, Sacanna S, Steinberg AS, Pine DJ, Chaikin PM. Living crystals of light-activated colloidal surfers. Science. 2013;339:936
  12. 12. Buttinoni I, Bialké J, Kümmel F, Löwen H, Bechinger C, Speck T. Dynamical clustering and phase separation in suspensions of self-propelled colloidal particles. Physical Review Letters. 2013;110(1):238301
  13. 13. He X, Aizenberg M, Kuksenok O, Zarzar LD, Shastri A, Balazs AC, et al. Synthetic homeostatic materials with chemo-mechano-chemical self-regulation. Nature. 2012;487:214
  14. 14. Li J, Shklyaev OE, Li T, Liu W, Shum H, Rozen I, et al. Self-propelled nanomotors autonomously seek and repair cracks. Nano Letters. 2015;15:7077
  15. 15. Uspal WE, Popescu MN, Dietrich S, Tasinkevych M. Rheotaxis of spherical active particles near a planar wall. Soft Matter. 2015;11:6613
  16. 16. Bechinger C, Di Leonardo R, Löwen H, Reichhardt C, Volpe G, Volpe G. Active particles in complex and crowded environments. Reviews of Modern Physics. 2016;88(1):045006
  17. 17. Paxton WF, Kistler KC, Olmeda CC, Sen A, St SK, Angelo YY, et al. Catalytic nanomotors: Autonomous movement of striped nanorods. Journal of the American Chemical Society. 2004;126:13424
  18. 18. Howse JR, Jones RAL, Ryan AJ, Gough T, Vafabakhsh R, Golestanian R. Self-motile colloidal particles: From directed propulsion to random walk. Physical Review Letters. 2007;99(1):048102
  19. 19. Brown A, Poon W. Ionic effects in self-propelled Pt-coated Janus swimmers. Soft Matter. 2014;10:4016
  20. 20. Ebbens S, Gregory DA, Dunderdale G, Howse JR, Ibrahim Y, Liverpool TB, et al. Electrokinetic effects in catalytic pt-insulator janus swimmers. EPL. 2014;106(1):58003
  21. 21. Brooks AM, Tasinkevych M, Sabrina S, Velegol D, Sen A. Shape-directed rotation of homogeneous micromotors via catalytic self-electrophoresis. Nature Communications. 2019;10:495
  22. 22. Derjaguin BV, Yalamov YI, Storozhilova AI. Diffusiophoresis of large aerosol particles. Journal of Colloid and Interface Science. 1966;22:117
  23. 23. Anderson JL. Colloid transport by interfacial forces. Annual Review of Fluid Mechanics. 1989;21:61
  24. 24. Golestanian R, Liverpool T, Ajdari A. Designing phoretic micro- and nano-swimmers. New Journal of Physics. 2007;9(1):126
  25. 25. Moran JL, Posner JD. Electrokinetic locomotion due to reaction-induced charge auto-electrophoresis. Journal of Fluid Mechanics. 2011;680:31
  26. 26. de Graaf J, Rempfer G, Holm C. Diffusiophoretic self-propulsion for partially catalytic spherical colloids. IEEE Transactions on Nanobioscience. 2015;14:272
  27. 27. Brown AT, Poon WCK, Holm C, de Graaf J. Ionic screening and dissociation are crucial for understanding chemical self-propulsion in polar solvents. Soft Matter. 2017;13:1200
  28. 28. Popescu MN, Dietrich S, Oshanin G. Confinement effects on diffusiophoretic self-propellers. The Journal of Chemical Physics. 2009;130(1):194702
  29. 29. Crowdy DG. Wall effects on self-diffusiophoretic Janus particles: A theoretical study. Journal of Fluid Mechanics. 2013;735:473
  30. 30. Uspal WE, Popescu MN, Dietrich S, Tasinkevych M. Self-propulsion of a catalytically active particle near a planar wall: From reflection to sliding and hovering. Soft Matter. 2015;11:434
  31. 31. Das S, Garg A, Campbell AI, Howse JR, Sen A, Velegol D, et al. Boundaries can steer active Janus spheres. Nature Communications. 2015;6(1):8999
  32. 32. Simmchen J, Katuri J, Uspal WE, Popescu MN, Tasinkevych M, Sánchez S. Topographical pathways guide chemical microswimmers. Nature Communications. 2016;7(1):10598
  33. 33. Liu C, Zhou C, Wang W, Zhang HP. Bimetallic microswimmers speed up in confining channels. Physical Review Letters. 2016;117(1):198001
  34. 34. Mozaffari A, Sharifi-Mood N, Koplik J, Maldarelli C. Self-propelled colloidal particle near a planar wall: A Brownian dynamics study. Physical Review Fluids. 2018;3:014104
  35. 35. Popescu MN, Uspal WE, Bechinger C, Fischer P. Chemotaxis of active janus nanoparticles. Nano Letters. 2018;18:5345
  36. 36. Palacci J, Sacanna S, Abramian A, Barral J, Hanson K, Grosberg AY, et al. Artificial rheotaxis. Science Advances. 2015;1(1):e1400214
  37. 37. Stone HA, Samuel ADT. Propulsion of microorganisms by surface distortions. Physical Review Letters. 1996;77:4102
  38. 38. Lighthill MJ. On the squirming motion of nearly spherical deformable bodies through liquids at very small Reynolds numbers. Communications on Pure and Applied Mathematics. 1952;5:109
  39. 39. Blake JR. A spherical envelope approach to ciliary propulsion. Journal of Fluid Mechanics. 1971;46:199
  40. 40. Ishikawa T, Locsei JT, Pedley TJ. Development of coherent structures in concentrated suspensions of swimming model micro-organisms. Journal of Fluid Mechanics. 2008;615:401
  41. 41. Ishimoto K, Gaffney EA. Squirmer dynamics near a boundary. Physical Review E. 2013;88:062702
  42. 42. Katuri J, Uspal WE, Simmchen J, López AM, Sánchez S. Cross-stream migration of active particles. Science Advances. 2018;4:eaao1755
  43. 43. Uspal WE. Theory of light-activated catalytic Janus particles. The Journal of Chemical Physics. 2019;150:114903
  44. 44. Schaar K, Zöttl A, Stark H. Detention times of microswimmers close to surfaces: Influence of hydrodynamic interactions and noise. Physical Review Letters. 2015;115:038101
  45. 45. Mirkovic T, Zacharia NS, Scholes GD, Ozin GA. Nanolocomotion - Catalytic Nanomotors and Nanorotors. Small. 2010;6:159
  46. 46. Ebbens SJ, Howse JR. In pursuit of propulsion at the nanoscale. Soft Matter. 2010;6:726
  47. 47. Lee T, Alarc’on-Correa M, Miksch C, Hahn K, Gibbs JG, Fischer P. Self-Propelling Nanomotors in the Presence of Strong Brownian Forces. Nano Letters. 2014;14:2407-2412
  48. 48. Reigh SY, Huang M, Schofield J, Kapral R. Microscopic and continuum descriptions of Janus motor fluid flow fields. Philosophical Transactions of the Royal Society A. 2016;374
  49. 49. de Buyl P, Mikhailov AS, Kapral R. Self-propulsion through symmetry breaking. EPL. 2013;103
  50. 50. Dey KK, Das S, Poyton MF, Sengupta S, Butler PJ, Cremer PS, et al. Chemotactic Separation of Enzymes. ACS Nano. 2014;8
  51. 51. Pozrikidis C. A Practical Guide to Boundary Element Methods with the Software Library BEMLIB. Boca Raton: CRC Press; 2002
  52. 52. Spagnolie SE, Lauga E. Hydrodynamics of self-propulsion near a boundary: Predictions and accuracy of far-field approximations. Journal of Fluid Mechanics. 2012;700:105
  53. 53. Shum H, Gaffney EA. Hydrodynamic analysis of flagellated bacteria swimming in corners of rectangular channels. Physical Review E. 2015;92:063016
  54. 54. Montenegro-Johnson TD, Michelin S, Lauga E. A regularised singularity approach to phoretic problems. European Physical Journal E: Soft Matter and Biological Physics. 2015;38:139
  55. 55. Schmiedling L, Lauga E, Montenegro-Johnson TD. Autophoretic flow on a torus. Physical Review Fluids. 2017;2:034201
  56. 56. Singh DP,Uspal WE, Popescu MN, Wilson LG and Fischer P. Photo-gravitactic microswimmers. Advanced Functional Materials. 2018;28:1706660
  57. 57. Varma A, Montenegro-Johnson TD, Michelin S. Clustering-induced self-propulsion of isotropic autophoretic particles. Soft Matter. 2018;14:7155
  58. 58. Popescu MN, Uspal WE, Dietrich S. Self-diffusiophoresis of chemically active colloids. The European Physical Journal Special Topics. 2016;225:2189
  59. 59. Kim S, Karrila SJ. Microhydrodynamics: Principles and Selected Applications. Mineola, NY: Dover Publications; 2005
  60. 60. Jackson JD. Classical Electrodynamics. 3rd ed. Hoboken, NJ, USA: Wiley; 1999
  61. 61. Cortez R. The method of regularized Stokeslets. SIAM Journal on Scientific Computing. 2001;23:1024
  62. 62. Aderogba K, Blake J. Action of a force near the planar surface between two semi-infinite immiscible liquids at very low Reynolds numbers. Bulletin of the Australian Mathematical Society. 1978;18:345
  63. 63. Mathijssen AJTM, Doostmohammadi A, Yeomans JM, Shendruk TN. Hotspots of boundary accumulation: Dynamics and statistics of micro-swimmers in flowing films. Journal of the Royal Society Interface. 2016;13:20150936
  64. 64. Delong S, Usabiaga FB, Delgado-Buscalioni R, Griffith BE, Donev A. Brownian dynamics without Green’s functions. The Journal of Chemical Physics. 2014;140:134110
  65. 65. Rapaport DC. Molecular dynamics simulation using quaternions. Journal of Computational Physics. 1985;60:306
  66. 66. Wittkowski R, Löwen H. Self-propelled Brownian spinning top: Dynamics of a biaxial swimmer at low Reynolds numbers. Physical Review E. 2012;85:021406
  67. 67. Beard DA, Schlick T. Unbiased rotational moves for rigid-body dynamics. Biophysical Journal. 2003;85:2973
  68. 68. Jones RB, Alavi FN. Rotational diffusion of a tracer colloid particle: IV. Brownian dynamics with wall effects. Physica A. 1992;187:436-455
  69. 69. Lisicki M, Cichocki B, Rogers SA, Dhont JKG, Lang PR. Translational and rotational near-wall diffusion of spherical colloids studied by evanescent wave scattering. Soft Matter. 2014;10:4312
  70. 70. Fixman M. Implicit algorithm for brownian dynamics of polymers. Macromolecules. 1986;19:1195
  71. 71. Bao Y, Rachh M, Keaveny EE, Greengard L, Donev A. A fluctuating boundary integral method for Brownian suspensions. Journal of Computational Physics. 2018;334
  72. 72. Uspal WE, Popescu MN, Tasinkevych M, Dietrich S. Shape-dependent guidance of active Janus particles by chemically patterned surfaces. New Journal of Physics. 2018;20:015013
  73. 73. Brady JS, Bossis G. Stokesian dynamics. Annual Review of Fluid Mechanics. 1988;20:111
  74. 74. Kim AS, Stolzenbach KD. The Permeability of Synthetic Fractal Aggregates with Realistic Three-Dimensional Structure. Journal of Colloid and Interface Science. 2002;253:315328
  75. 75. Kim AS, Stolzenbach KD. Aggregate formation and collision efficiency in differential settling. Journal of Colloid and Interface Science. 2004;271:110
  76. 76. Yan W, Brady JF. The force on a boundary in active matter. Journal of Fluid Mechanics. 2015;785:1
  77. 77. Walter J, Salsac A-V, Barthés-Biesel D, Le Tallec P. Coupling of finite element and boundary integral methods for a capsule in a Stokes flow. International Journal for Numerical Methods in Engineering. 2010;83:829
  78. 78. Gao T, Hu HH. Deformation of elastic particles in viscous shear flow. Journal of Computational Physics. 2009;228:2132
  79. 79. Zhao H, Isfahani AHG, Olson LN, Freund JB. A spectral boundary integral method for flowing blood cells. Journal of Computational Physics. 2010;229:3726
  80. 80. Würger ACR. Is Soret equilibrium a non-equilibrium effect?. Mécanique. 2013;341:438
  81. 81. Nikola N, Solon AP, Kafri Y, Kardar M, Tailleur J, Voituriez R. Active particles with soft and curved walls: equation of state, ratchets, and instabilities. Physical Review Letters. 2016;117:098001
  82. 82. Montenegro-Johnson TD. Microtransformers: Controlled microscale navigation with flexible robots. Physical Review Fluids. 2018;3:062201(R)
  83. 83. Batchelor GK, Green JT. The hydrodynamic interaction of two small freely-moving spheres in a linear flow field. Journal of Fluid Mechanics. 1972;56:375
  84. 84. Lauga E, Michelin S. Stresslets induced by active swimmers. Physical Review Letters. 2016;117:148001

Written By

William E. Uspal

Submitted: October 3rd, 2018 Reviewed: May 8th, 2019 Published: June 23rd, 2019