The Boundary Element Method for Fluctuating Active Colloids DOI: http://dx.doi.org/10.5772/intechopen.86738

The squirming mode amplitudes Bn,α, 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<sup>1</sup> determines the velocity of an isolated squirmer in unbounded solution: Ufs,<sup>α</sup> ¼ ð Þ 2=3 B1,α. According to our definition of d<sup>α</sup> and θp,α, we require that B1,α>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

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

H2O2 !Pt

the production of oxygen with zero order kinetics:

H2O þ

(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

where D is the diffusion coefficient of oxygen, cð Þ x is the number density of oxygen, and κð Þ x<sup>s</sup> 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 ¼ κ0R=D, where κ<sup>0</sup> is a characteristic reaction rate; we assume Da ≪ 1.) Furthermore, we assume that the Péclet number Pe � U0R=D is very small, where U<sup>0</sup> is a characteristic particle velocity and R is the particle radius. Accordingly, we can make a quasi-steady approximation for the diffusion of oxygen

∇2

where c<sup>∞</sup> 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

1 2

�D½ �j <sup>∇</sup><sup>c</sup> � <sup>n</sup>^ <sup>x</sup>¼x<sup>s</sup> <sup>¼</sup> <sup>κ</sup>ð Þ <sup>x</sup><sup>s</sup> , (11)

c ¼ 0: (12)

cð Þ¼ jxj ! ∞j c∞, (13)

O2: (10)

squirmer decay rapidly with distance from the squirmer.

Non-Equilibrium Particle Dynamics

oxygen:

in the solution:

54

Finally, we assume that

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<sup>s</sup> 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<sup>þ</sup> s � �. Here, the plus sign emphasizes that c x<sup>þ</sup> 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=U0.) Accordingly, the local pressure P x� <sup>s</sup> ; <sup>η</sup> � � can be calculated from <sup>Φ</sup> and <sup>c</sup> <sup>x</sup><sup>þ</sup> s � �, where η is a coordinate defined at xs 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<sup>s</sup> 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:

$$\mathbf{v}\_{\mathfrak{s},a}(\mathbf{x}\_{\mathfrak{s}}) = -b(\mathbf{x}\_{\mathfrak{s}})\nabla\_{\parallel}\mathfrak{c}.\tag{14}$$

Here, the surface gradient operator is defined as ∇<sup>k</sup> � ð Þ� I � n^n^ ∇. The material-dependent parameter bð Þ x<sup>s</sup> 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 σ; σ<sup>0</sup> ð Þ and velocity fields u; u<sup>0</sup> ð Þ of two solutions to the Stokes equation within the same domain V:

$$\int\_{\mathcal{S}} \mathbf{u} \cdot \boldsymbol{\sigma}' \cdot \hat{\mathbf{n}} \cdot d\mathbf{S} = \int\_{\mathcal{S}} \mathbf{u}' \cdot \boldsymbol{\sigma} \cdot \hat{\mathbf{n}} \cdot d\mathbf{S},\tag{15}$$

where S is the boundary of V. For the N active particles in unbounded solution, <sup>S</sup> <sup>¼</sup> <sup>∪</sup><sup>N</sup> <sup>α</sup>¼<sup>1</sup>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 Fext,<sup>α</sup> ¼ 0 and τext,<sup>α</sup> ¼ 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 6N "primed" problems, indexed by j ¼ 1, 2, …, 6N. For problem j ¼ α, particle α is exposed to an external force with magnitude F<sup>0</sup> ext in the x^ direction. Each particle has unknown translational and rotational velocities U0ð Þ<sup>j</sup> <sup>α</sup> and <sup>Ω</sup>0ð Þ<sup>j</sup> <sup>α</sup> , and the velocity field u0ð Þ<sup>j</sup> is subject to no-slip boundary conditions on each particle:

$$\mathbf{u}'^{(j)}(\mathbf{x}\_{\mathfrak{s}}) = \mathbf{U}\_{a}'^{(j)} + \mathfrak{Q}\_{a}'^{(j)} \times (\mathbf{x}\_{\mathfrak{t}} - \mathbf{x}\_{a}), \ \mathbf{x}\_{\mathfrak{t}} \in \mathcal{S}\_{a\mathfrak{b}} \tag{16}$$

but the integral is simply the force F<sup>α</sup> on particle α. Since the particles are free of

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

where b is a vector containing the 6N integrals associated with the right hand side of Eq. 21, <sup>V</sup> is vector of all 6<sup>N</sup> velocity components ð Þ <sup>U</sup>α; <sup>Ω</sup><sup>α</sup> <sup>T</sup>, and <sup>R</sup> is the grand resistance matrix for N spheres at positions xα. Note that the arbitrary

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

> ∂ui ∂xj þ ∂uj ∂xi

δij þ 2μeij <sup>0</sup> � �eij

> eii þ 2μe 0 ijeij

eij <sup>¼</sup> <sup>1</sup> 2

ijeij:

ijeij ¼ �P<sup>0</sup>

¼ �P<sup>0</sup>

¼ 2μe 0 ijeij,

where we have used ∇ � u ¼ 0 to eliminate eii, and we assume the Einstein convention for summation over repeated indices. Similarly, we can obtain

σ0

and that, in index notation, the stress tensor is

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 6N velocities for any set of slip velocities vs,α. This, for instance, facilitates studying how various choices for bð Þ x<sup>s</sup> or Bn,<sup>α</sup> affect particle motion. Additionally, the "primed" problem for the resistance matrix R and stresses σ0ð Þ<sup>j</sup> 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 6N particle velocities in the presence of the force-free and torque-free constraints.

α ð Sα <sup>α</sup>ð Þj can be rearranged as

ð Þ� x<sup>s</sup> � x<sup>α</sup> σ � n^ dS, (20)

<sup>v</sup><sup>s</sup> � <sup>σ</sup>0ð Þ<sup>j</sup> � <sup>n</sup>^ dS , j <sup>¼</sup> <sup>1</sup>, …, <sup>6</sup><sup>N</sup> : (21)

R � V ¼ b, (22)

� �, (23)

(25)

σij ¼ �Pδij þ 2μeij (24)

external forces, F<sup>α</sup> ¼ 0. Likewise, the term involving Ω<sup>0</sup>

The Boundary Element Method for Fluctuating Active Colloids

DOI: http://dx.doi.org/10.5772/intechopen.86738

Ω0 <sup>α</sup>ð Þ� j ð Sα

but the integral is the torque τ<sup>α</sup> ¼ 0 on particle α.

α � � ¼ � <sup>∑</sup>

These 6N equations can be written in matrix form:

<sup>α</sup> <sup>þ</sup> <sup>Ω</sup><sup>α</sup> � <sup>τ</sup>0ð Þ<sup>j</sup>

∑ α

magnitudes F<sup>0</sup>

<sup>U</sup><sup>α</sup> � <sup>F</sup>0ð Þ<sup>j</sup>

ext and τ<sup>0</sup>

2.3.1 Proof of Lorentz reciprocal theorem

strain tensor eij is defined as

We consider the quantity σ<sup>0</sup>

ijeij, so that

σije<sup>0</sup>

57

ij ¼ 2μe<sup>0</sup>

Additionally, the flow field vanishes far away from the particles, i.e., <sup>u</sup>0ð Þ<sup>j</sup> ð Þ¼ <sup>j</sup>xj ! <sup>∞</sup> 0. The unprimed problem and primed problem <sup>α</sup> 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 τ<sup>0</sup> ext in the x^, ^y, and ^z directions, respectively, with the same boundary conditions. Each "primed" problem j can be solved for 6N unknown velocity components U0ð Þ<sup>j</sup> <sup>α</sup> and <sup>Ω</sup>0ð Þ<sup>j</sup> <sup>α</sup> .

For problem <sup>j</sup>, we substitute Eqs. 3 and 16 into Eq. 15 for <sup>u</sup> and <sup>u</sup><sup>0</sup> <sup>¼</sup> <sup>u</sup>0ð Þ<sup>j</sup> to obtain:

$$\sum\_{a} \int\_{S\_a} \left[ \mathbf{U}\_a + \mathbf{\Omega}\_a \times (\mathbf{x}\_t - \mathbf{x}\_a) + \mathbf{v}\_{\mathbf{s},a}(\mathbf{x}\_t) \right] \cdot \boldsymbol{\sigma}^{(j)} \cdot \hat{\mathbf{n}} \cdot d\mathbf{S} = \tag{17}$$

$$\sum\_{a} \int\_{\mathcal{S}\_{a}} \left[ \mathbf{U}\_{a}^{\prime}(j) + \mathbf{M}\_{a}^{\prime}(j) \times (\mathbf{x}\_{t} - \mathbf{x}\_{a}) \right] \cdot \boldsymbol{\sigma} \cdot \hat{\mathbf{n}} \,^{\prime} \, d\mathbf{S}. \tag{18}$$

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

$$\mathbf{U}'\_a(j) \cdot \int\_{S\_a} \boldsymbol{\sigma} \cdot \hat{\mathbf{n}} \ \mathrm{d}S,\tag{19}$$

### 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<sup>0</sup> ext in the x^ direction. Similarly, in primed problems α þ 1 and α þ 2, particle α is exposed to a force with magnitude F<sup>0</sup> 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 τ<sup>0</sup> ext in the x^, ^y, and ^z directions, respectively. Note that the primed and the 6N unprimed problems all have the same geometry, i.e., the same configuration of N spheres.

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 6N "primed" problems, indexed by j ¼ 1, 2, …, 6N. For problem j ¼ α,

the velocity field u0ð Þ<sup>j</sup> is subject to no-slip boundary conditions on each particle:

<sup>α</sup> <sup>þ</sup> <sup>Ω</sup>0ð Þ<sup>j</sup>

Additionally, the flow field vanishes far away from the particles, i.e.,

<sup>u</sup>0ð Þ<sup>j</sup> ð Þ¼ <sup>j</sup>xj ! <sup>∞</sup> 0. The unprimed problem and primed problem <sup>α</sup> 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

directions, respectively, with the same boundary conditions. Each "primed" prob-

For problem <sup>j</sup>, we substitute Eqs. 3 and 16 into Eq. 15 for <sup>u</sup> and <sup>u</sup><sup>0</sup> <sup>¼</sup> <sup>u</sup>0ð Þ<sup>j</sup> to

<sup>α</sup>ð Þ�j ð Þ x<sup>s</sup> � x<sup>α</sup>

It can be shown that the right hand side of this equation vanishes. Consider the

Illustration of the "unprimed" problem for the velocities of N active particles, and the "primed" problem α for

direction and the z^ direction, respectively. Moreover, in primed problems α þ 3, α þ 4, and α þ 5, particle α is

ext in the x^, ^y, and ^z directions, respectively. Note that the primed and the

the velocities of N inert particles when particle α is exposed to a force with magnitude F<sup>0</sup>

Similarly, in primed problems α þ 1 and α þ 2, particle α is exposed to a force with magnitude F<sup>0</sup>

6N unprimed problems all have the same geometry, i.e., the same configuration of N spheres.

<sup>U</sup><sup>α</sup> <sup>þ</sup> <sup>Ω</sup><sup>α</sup> � ð Þþ <sup>x</sup><sup>s</sup> � <sup>x</sup><sup>α</sup> <sup>v</sup>s,αð Þ <sup>x</sup><sup>s</sup> <sup>½</sup> � � <sup>σ</sup>0ð Þ<sup>j</sup> � <sup>n</sup>^ dS <sup>¼</sup> (17)

� � � <sup>σ</sup> � <sup>n</sup>^ dS: (18)

ext in the x^ direction.

ext in the x^, ^y, and ^z

<sup>α</sup> .

ext in the x^ direction.

ext in the ^y

<sup>α</sup> and <sup>Ω</sup>0ð Þ<sup>j</sup>

<sup>α</sup>ð Þj is a constant and can be

σ � n^ dS, (19)

<sup>α</sup> � ð Þ x<sup>s</sup> � x<sup>α</sup> , x<sup>s</sup> ∈ Sα, (16)

<sup>α</sup> and <sup>Ω</sup>0ð Þ<sup>j</sup>

<sup>α</sup> , and

particle α is exposed to an external force with magnitude F<sup>0</sup>

ð Þ¼ xs <sup>U</sup>0ð Þ<sup>j</sup>

j ¼ α þ 5, particle α is exposed to a torque with magnitude τ<sup>0</sup>

lem j can be solved for 6N unknown velocity components U0ð Þ<sup>j</sup>

<sup>α</sup>ð Þþj Ω<sup>0</sup>

<sup>α</sup>ð Þj . For each integral over Sα, U<sup>0</sup>

U0 <sup>α</sup>ð Þ� j ð Sα

u0ð Þ<sup>j</sup>

Non-Equilibrium Particle Dynamics

∑ α ð Sα

> ∑ α ð Sα U0

obtain:

Figure 2.

56

exposed to torques with magnitude τ<sup>0</sup>

term involving U<sup>0</sup>

moved out of the integral,

Each particle has unknown translational and rotational velocities U0ð Þ<sup>j</sup>

but the integral is simply the force F<sup>α</sup> on particle α. Since the particles are free of external forces, F<sup>α</sup> ¼ 0. Likewise, the term involving Ω<sup>0</sup> <sup>α</sup>ð Þj can be rearranged as

$$\mathbf{M}\_a'(j) \cdot \int\_{S\_a} (\mathbf{x}\_t - \mathbf{x}\_a) \times \boldsymbol{\sigma} \cdot \hat{\mathbf{n}} \ \mathrm{d}S,\tag{20}$$

but the integral is the torque τ<sup>α</sup> ¼ 0 on particle α.

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

$$\sum\_{a} \left( \mathbf{U}\_{a} \cdot \mathbf{F}\_{a}^{(j)} + \mathbf{\Omega}\_{a} \cdot \mathbf{\sigma}\_{a}^{(j)} \right) = - \sum\_{a} \int\_{S\_{a}} \mathbf{v}\_{i} \cdot \mathbf{\sigma}^{(j)} \cdot \hat{\mathbf{n}} \ \text{d}\mathbf{S} \ , j = \mathbf{1}, \dots, \text{6N} \ . \tag{21}$$

These 6N equations can be written in matrix form:

$$
\mathcal{R} \cdot \mathbf{V} = \mathbf{b},
\tag{22}
$$

where b is a vector containing the 6N integrals associated with the right hand side of Eq. 21, <sup>V</sup> is vector of all 6<sup>N</sup> velocity components ð Þ <sup>U</sup>α; <sup>Ω</sup><sup>α</sup> <sup>T</sup>, and <sup>R</sup> is the grand resistance matrix for N spheres at positions xα. Note that the arbitrary magnitudes F<sup>0</sup> ext and τ<sup>0</sup> 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 6N velocities for any set of slip velocities vs,α. This, for instance, facilitates studying how various choices for bð Þ x<sup>s</sup> or Bn,<sup>α</sup> affect particle motion. Additionally, the "primed" problem for the resistance matrix R and stresses σ0ð Þ<sup>j</sup> 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 6N particle velocities in the presence of the force-free and torque-free constraints.
