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 eij is defined as

$$\mathfrak{e}\_{ij} = \frac{1}{2} \left( \frac{\partial u\_i}{\partial \mathfrak{x}\_j} + \frac{\partial u\_j}{\partial \mathfrak{x}\_i} \right), \tag{23}$$

and that, in index notation, the stress tensor is

$$
\sigma\_{\vec{\eta}} = -P\delta\_{\vec{\eta}} + 2\mu e\_{\vec{\eta}} \tag{24}
$$

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

$$
\sigma'\_{\vec{ij}} e\_{\vec{ij}} = \left(-P' \delta\_{\vec{ij}} + 2\mu e\_{\vec{ij}}{}^{\prime}\right) e\_{\vec{ij}}
$$

$$
= -P' e\_{\vec{ii}} + 2\mu e'\_{\vec{ij}} e\_{\vec{ij}} \tag{25}
$$

$$
= 2\mu e'\_{\vec{ij}} e\_{\vec{ij}},
$$

where we have used ∇ � u ¼ 0 to eliminate eii, and we assume the Einstein convention for summation over repeated indices. Similarly, we can obtain σije<sup>0</sup> ij ¼ 2μe<sup>0</sup> ijeij, so that

$$
\sigma'\_{i\dot{\jmath}} e\_{i\dot{\jmath}} = \sigma\_{i\dot{\jmath}} e'\_{i\dot{\jmath}}.\tag{26}
$$

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 Nelm elements,

In order to obtain the BIE for the Laplace equation, we first consider the diver-

S

S

S

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 ¼ ϕ∇ψ,

<sup>∇</sup> � <sup>A</sup>dV ¼ �<sup>ð</sup>

ψ þ ∇ϕ � ∇ψ, and we obtain Green's first identity:

<sup>ψ</sup> <sup>þ</sup> <sup>∇</sup><sup>ϕ</sup> � <sup>∇</sup><sup>ψ</sup> � �dV ¼ �<sup>ð</sup>

<sup>ϕ</sup> <sup>þ</sup> <sup>∇</sup><sup>ψ</sup> � <sup>∇</sup><sup>ϕ</sup> � �dV ¼ �<sup>ð</sup>

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

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 Vp, and the interior of particle α by

<sup>α</sup>¼<sup>1</sup> <sup>S</sup>α. The observation point <sup>x</sup><sup>0</sup> can occur in <sup>V</sup>, in <sup>V</sup>, or on <sup>S</sup>; we show <sup>x</sup><sup>0</sup> in <sup>V</sup>.

<sup>α</sup>¼<sup>1</sup> <sup>V</sup>p,α. The solid and fluid domains are separated by the particle surfaces <sup>S</sup>α, with

elm � �, and the required computation

A � n^ dS, (32)

ϕ∇ψ � n^ dS: (33)

ψ∇ϕ � n^ dS: (34)

the required computer memory scales as O N<sup>2</sup>

The Boundary Element Method for Fluctuating Active Colloids

ð V

where ϕð Þ x and ψð Þ x are scalar fields, then the divergence term

We can also write Green's first identity for A ¼ ψ∇ϕ:

elm � �.

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

ð V ϕ∇<sup>2</sup>

ð V ψ∇<sup>2</sup>

time scales as O N<sup>3</sup>

gence theorem:

<sup>∇</sup> � <sup>A</sup> <sup>¼</sup> <sup>ϕ</sup>∇<sup>2</sup>

Figure 3.

<sup>S</sup> <sup>¼</sup> <sup>∪</sup><sup>N</sup>

59

<sup>V</sup>p,α, where <sup>V</sup><sup>p</sup> <sup>¼</sup> <sup>∪</sup><sup>N</sup>

We can also manipulate σ<sup>0</sup> ijeij as follows:

$$
\sigma\_{ij}' e\_{ij} = \frac{1}{2} \sigma\_{ij}' \frac{\partial u\_i}{\partial \mathbf{x}\_j} + \frac{1}{2} \sigma\_{ij}' \frac{\partial u\_j}{\partial \mathbf{x}\_i} \,. \tag{27}
$$

Swapping the two indices in the last term,

$$
\sigma\_{ij}' \varepsilon\_{ij} = \frac{1}{2} \sigma\_{ij}' \frac{\partial u\_i}{\partial \mathbf{x}\_j} + \frac{1}{2} \sigma\_{ji}' \frac{\partial u\_i}{\partial \mathbf{x}\_j} \,. \tag{28}
$$

But σ<sup>0</sup> ij ¼ σ<sup>0</sup> ji, giving

$$\begin{split} \sigma'\_{\vec{\eta}} \mathbf{e}\_{\vec{\eta}} &= \sigma'\_{\vec{\eta}} \frac{\partial u\_i}{\partial \mathbf{x}\_j} \\ &= \frac{\partial}{\partial \mathbf{x}\_j} \left( \sigma'\_{\vec{\eta}} u\_i \right) - \left( \frac{\partial \sigma'\_{\vec{\eta}}}{\partial \mathbf{x}\_j} \right), \end{split} \tag{29}$$

so that

$$
\frac{
\partial
}{
\partial \mathbf{x}\_j
} \left( \sigma\_{ij}' u\_i \right) - \left( \frac{
\partial \sigma\_{ij}'
}{
\partial \mathbf{x}\_j
} \right) u\_i = \frac{
\partial
}{
\partial \mathbf{x}\_j
} \left( \sigma\_{ij} u\_i' \right) - \left( \frac{
\partial \sigma\_{ij}
}{
\partial \mathbf{x}\_j
} \right) u\_i'.
\tag{30}
$$

If there are no point forces applied to the fluid in determination of u<sup>0</sup> and u, then ∇ � σ ¼ 0 and ∇ � σ<sup>0</sup> ¼ 0, and we obtain

$$\nabla \cdot (\mathbf{u} \cdot \boldsymbol{\sigma}') = \nabla \cdot (\mathbf{u}' \cdot \boldsymbol{\sigma}).\tag{31}$$

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

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

σ0

ijeij as follows:

σ0 ijeij <sup>¼</sup> <sup>1</sup> 2 σ0 ij ∂ui ∂xj þ 1 2 σ0 ij ∂uj ∂xi

σ0 ijeij <sup>¼</sup> <sup>1</sup> 2 σ0 ij ∂ui ∂xj þ 1 2 σ0 ji ∂ui ∂xj

σ0 ijeij ¼ σ<sup>0</sup> ij ∂ui ∂xj

> ¼ ∂ ∂xj σ0 ij ui � �

� <sup>∂</sup>σ<sup>0</sup> ij ∂xj !

2.4 Boundary integral formulation of the Laplace equation

ui <sup>¼</sup> <sup>∂</sup> ∂xj

If there are no point forces applied to the fluid in determination of u<sup>0</sup> and u, then

Integrating both sides over the volume V and applying the divergence theorem,

Even with the aid of the Lorentz reciprocal theorem, it is necessary to solve the

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

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

Stokes and (for self-phoretic particles) Laplace equations in a fluid domain

Swapping the two indices in the last term,

We can also manipulate σ<sup>0</sup>

Non-Equilibrium Particle Dynamics

But σ<sup>0</sup>

so that

we obtain Eq. 15.

the particles.

58

ij ¼ σ<sup>0</sup>

ji, giving

∂ ∂xj σ0 ijui � �

∇ � σ ¼ 0 and ∇ � σ<sup>0</sup> ¼ 0, and we obtain

ijeij ¼ σije

0

� <sup>∂</sup>σ<sup>0</sup> ij ∂xj !

> σiju<sup>0</sup> i � � � <sup>∂</sup>σij

,

∂xj � � u0 i

∇ � u � σ<sup>0</sup> ð Þ¼ ∇ � u<sup>0</sup> ð Þ � σ : (31)

ij: (26)

: (27)

: (28)

(29)

: (30)

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 Nelm elements, the required computer memory scales as O N<sup>2</sup> elm � �, and the required computation time scales as O N<sup>3</sup> elm � �.

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

$$\int\_{\mathcal{V}} \nabla \cdot \mathbf{A} \, dV = -\int\_{\mathcal{S}} \mathbf{A} \cdot \hat{\mathbf{n}} \, \, d\mathbf{S},\tag{32}$$

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 <sup>∇</sup> � <sup>A</sup> <sup>¼</sup> <sup>ϕ</sup>∇<sup>2</sup> ψ þ ∇ϕ � ∇ψ, and we obtain Green's first identity:

$$\int\_{\mathcal{V}} \left( \phi \nabla^2 \boldsymbol{\psi} + \nabla \phi \cdot \nabla \boldsymbol{\psi} \right) dV = - \int\_{\mathcal{S}} \phi \nabla \boldsymbol{\psi} \cdot \hat{\mathbf{n}} \, \, d\mathbf{S}.\tag{33}$$

We can also write Green's first identity for A ¼ ψ∇ϕ:

$$\int\_{\mathcal{V}} \left( \boldsymbol{\psi} \nabla^2 \phi + \nabla \boldsymbol{\psi} \cdot \nabla \phi \right) dV = - \int\_{\mathcal{S}} \boldsymbol{\psi} \nabla \phi \cdot \hat{\mathbf{n}} \cdot d\mathbf{S}.\tag{34}$$

Subtracting Eq. 34 from Eq. 33, we obtain Green's second 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 Vp, and the interior of particle α by <sup>V</sup>p,α, where <sup>V</sup><sup>p</sup> <sup>¼</sup> <sup>∪</sup><sup>N</sup> <sup>α</sup>¼<sup>1</sup> <sup>V</sup>p,α. The solid and fluid domains are separated by the particle surfaces <sup>S</sup>α, with <sup>S</sup> <sup>¼</sup> <sup>∪</sup><sup>N</sup> <sup>α</sup>¼<sup>1</sup> <sup>S</sup>α. The observation point <sup>x</sup><sup>0</sup> can occur in <sup>V</sup>, in <sup>V</sup>, or on <sup>S</sup>; we show <sup>x</sup><sup>0</sup> in <sup>V</sup>.

$$\int\_{V} \left( \phi \nabla^{2} \boldsymbol{\varphi} - \boldsymbol{\varphi} \nabla^{2} \phi \right) dV = - \int\_{S} \left( \phi \nabla \boldsymbol{\varphi} - \boldsymbol{\varphi} \nabla \phi \right) \cdot \hat{\mathbf{n}} \cdot d\mathbf{S}.\tag{35}$$

Now, we let <sup>ϕ</sup> <sup>¼</sup> <sup>c</sup>ð Þ <sup>x</sup> , with <sup>∇</sup><sup>2</sup><sup>c</sup> <sup>¼</sup> 0. Furthermore, we let <sup>ψ</sup> <sup>¼</sup> <sup>G</sup>ð Þ <sup>x</sup>; <sup>x</sup><sup>0</sup> , where the Green's function Gð Þ x; x<sup>0</sup> satisfies Poisson's equation:

$$
\nabla^2 G(\mathbf{x}, \mathbf{x}\_0) + \delta(\mathbf{x} - \mathbf{x}\_0) = \mathbf{0}.\tag{36}
$$

a priori. Eq. 41 has an interesting physical interpretation: cð Þ x<sup>0</sup> 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

We still have two other options for where to place x0: inside the boundary S (i.e., outside the solution domain V) or somewhere on the boundary S. If we place

Placing x<sup>0</sup> 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

This is a boundary integral equation (BIE) because the left hand side is the concentration cð Þ x<sup>0</sup> 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

a point x<sup>0</sup> on the surface, we obtain the boundary integral equation:

∇G x; x<sup>j</sup>

The Nelm equations can be written in matrix form:

<sup>∇</sup>Gð Þ� <sup>x</sup>; <sup>x</sup><sup>0</sup> <sup>n</sup>^ dS � � � ð Þ <sup>∇</sup><sup>c</sup> � <sup>n</sup>^ <sup>i</sup>

� � � <sup>n</sup>^ dS � � � ð Þ <sup>∇</sup><sup>c</sup> � <sup>n</sup>^ <sup>i</sup>

Aij � <sup>1</sup> 2

Aij � ð Si

Choosing x<sup>0</sup> as the midpoint x<sup>j</sup> of element j, we can write Nelm equations:

<sup>δ</sup>ij � �cj <sup>¼</sup> Bijð Þ <sup>∇</sup><sup>c</sup> � <sup>n</sup>^ <sup>j</sup>

∇G x; x<sup>j</sup>

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

½cð Þ x ∇Gð Þ� x; x<sup>0</sup> Gð Þ x; x<sup>0</sup> ∇cð Þ x � � n^ dS ¼ 0: (42)

½cð Þ x ∇Gð Þ� x; x<sup>0</sup> Gð Þ x; x<sup>0</sup> ∇cð Þ x � � n^ dS: (43)

ð Si

ð Si

G x; x<sup>j</sup>

� � � <sup>n</sup>^ dS (47)

, (46)

<sup>G</sup>ð Þ <sup>x</sup>; <sup>x</sup><sup>0</sup> dS � � � � : (44)

� � dS � � � � : (45)

x<sup>0</sup> inside S, then the integral on the left hand side of Eq. 37 vanishes:

and orientation n^ on the boundaries.

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

ð S

1 2

cð Þ¼ x<sup>0</sup>

ð S

The Boundary Element Method for Fluctuating Active Colloids

we integrate over V:

obtained with Eq. 43.

cð Þ¼ x<sup>0</sup> ∑

cj ¼ ∑ Nelm i¼1 ci ð Si

Nelm i¼1 ci ð Si

1 2

> 1 2

where

and

61

We obtain:

$$\int\_{\mathcal{V}} c(\mathbf{x}) \, \nabla^2 G(\mathbf{x}, \mathbf{x}\_0) \, dV = -\int\_{\mathcal{S}} \left[ c(\mathbf{x}) \nabla G(\mathbf{x}, \mathbf{x}\_0) - G(\mathbf{x}, \mathbf{x}\_0) \nabla c(\mathbf{x}) \right] \cdot \hat{\mathbf{n}} \, \, d\mathbf{S}.\tag{37}$$

We have not yet specified the location of the pole x0. If x<sup>0</sup> 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<sup>0</sup> at a point x<sup>0</sup> ∈V:

$$\varepsilon(\mathbf{x}\_0) = \int\_{\mathcal{S}} \left[ \varepsilon(\mathbf{x}) \nabla G(\mathbf{x}, \mathbf{x}\_0) - G(\mathbf{x}, \mathbf{x}\_0) \nabla c(\mathbf{x}) \right] \cdot \hat{\mathbf{n}} \, dS. \tag{38}$$

Using the divergence theorem, can show that ∇<sup>2</sup> <sup>1</sup> ∣x�x0∣ � � ¼ �4πδð Þ <sup>x</sup> � <sup>x</sup><sup>0</sup> , so that

$$G(\mathbf{x}, \mathbf{x}\_0) = \frac{1}{4\pi|\mathbf{x} - \mathbf{x}\_0|}. \tag{39}$$

We recall from electrostatics that Gð Þ x; x<sup>0</sup> represents the electrostatic potential at x from a point charge of unit strength located at x0. Within the context of the steady diffusion equation <sup>∇</sup><sup>2</sup><sup>c</sup> <sup>¼</sup> 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<sup>0</sup> and with unit strength (i.e., unit number density flux per unit time). One can take derivatives of Gð Þ x; x<sup>0</sup> with respect to x<sup>0</sup> to obtain higher order multipole singularities. For instance, we can obtain the Green's function Gdpð Þ x; x<sup>0</sup> for a source/sink dipole located at x<sup>0</sup> as

$$\mathbf{G}\_{\rm dp}(\mathbf{x}, \mathbf{x}\_0) \equiv \nabla\_{\mathbf{x}\_0} G(\mathbf{x}, \mathbf{x}\_0) = \frac{(\mathbf{x} - \mathbf{x}\_0)}{4\pi |\mathbf{x} - \mathbf{x}\_0|^3}. \tag{40}$$

As Gdpð Þ x; x<sup>0</sup> 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<sup>0</sup> = Gð Þ x0; x , so we can rewrite Eq. 38 as

$$c(\mathbf{x}\_0) = \int\_{\mathcal{S}} \left[ c(\mathbf{x}) \nabla G(\mathbf{x}\_0, \mathbf{x}) - G(\mathbf{x}\_0, \mathbf{x}) \nabla c(\mathbf{x}) \right] \cdot \hat{\mathbf{n}} \ \mathrm{d}\mathbf{S}.\tag{41}$$

Interestingly, we have obtained an expression for cð Þ x<sup>0</sup> 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<sup>0</sup> , 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^

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

ð V ϕ∇<sup>2</sup>

Non-Equilibrium Particle Dynamics

We obtain:

dipole located at x<sup>0</sup> as

60

strength and orientation of the dipole.

Gð Þ x0; x , so we can rewrite Eq. 38 as

cð Þ¼ x<sup>0</sup>

ð S

ð V <sup>c</sup>ð Þ <sup>x</sup> <sup>∇</sup><sup>2</sup>

<sup>ψ</sup> � <sup>ψ</sup>∇<sup>2</sup> <sup>ϕ</sup> � �dV ¼ �

the Green's function Gð Þ x; x<sup>0</sup> satisfies Poisson's equation:

∇2

ð S

Gð Þ x; x<sup>0</sup> dV ¼ �

cð Þ¼ x<sup>0</sup>

ð S

Using the divergence theorem, can show that ∇<sup>2</sup> <sup>1</sup>

Gð Þ¼ x; x<sup>0</sup>

Gdpð Þ� x; x<sup>0</sup> ∇<sup>x</sup>0Gð Þ¼ x; x<sup>0</sup>

As Gdpð Þ x; x<sup>0</sup> 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

By inspection, the Green's function obeys the symmetry property Gð Þ x; x<sup>0</sup> =

Interestingly, we have obtained an expression for cð Þ x<sup>0</sup> 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<sup>0</sup> , 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^

ð S

Now, we let <sup>ϕ</sup> <sup>¼</sup> <sup>c</sup>ð Þ <sup>x</sup> , with <sup>∇</sup><sup>2</sup><sup>c</sup> <sup>¼</sup> 0. Furthermore, we let <sup>ψ</sup> <sup>¼</sup> <sup>G</sup>ð Þ <sup>x</sup>; <sup>x</sup><sup>0</sup> , where

We have not yet specified the location of the pole x0. If x<sup>0</sup> 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<sup>0</sup> at a point x<sup>0</sup> ∈V:

ðϕ∇ψ � ψ∇ϕÞ � n^ dS: (35)

Gð Þþ x; x<sup>0</sup> δð Þ¼ x � x<sup>0</sup> 0: (36)

½cð Þ x ∇Gð Þ� x; x<sup>0</sup> Gð Þ x; x<sup>0</sup> ∇cð Þ x � � n^ dS: (37)

½cð Þ x ∇Gð Þ� x; x<sup>0</sup> Gð Þ x; x<sup>0</sup> ∇cð Þ x � � n^ dS: (38)

∣x�x0∣ � �

ð Þ x � x<sup>0</sup> 4πj j x � x<sup>0</sup>

½cð Þ x ∇Gð Þ� x0; x Gð Þ x0; x ∇cð Þ x � � n^ dS: (41)

1 4π∣x � x0∣

We recall from electrostatics that Gð Þ x; x<sup>0</sup> represents the electrostatic potential at x from a point charge of unit strength located at x0. Within the context of the steady diffusion equation <sup>∇</sup><sup>2</sup><sup>c</sup> <sup>¼</sup> 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<sup>0</sup> and with unit strength (i.e., unit number density flux per unit time). One can take derivatives of Gð Þ x; x<sup>0</sup> with respect to x<sup>0</sup> to obtain higher order multipole singularities. For instance, we can obtain the Green's function Gdpð Þ x; x<sup>0</sup> for a source/sink

¼ �4πδð Þ x � x<sup>0</sup> , so that

<sup>3</sup> : (40)

: (39)

a priori. Eq. 41 has an interesting physical interpretation: cð Þ x<sup>0</sup> 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 x0: inside the boundary S (i.e., outside the solution domain V) or somewhere on the boundary S. If we place x<sup>0</sup> inside S, then the integral on the left hand side of Eq. 37 vanishes:

$$\int\_{\mathcal{S}} \left[ \boldsymbol{c}(\mathbf{x}) \nabla G(\mathbf{x}, \mathbf{x}\_0) - G(\mathbf{x}, \mathbf{x}\_0) \nabla c(\mathbf{x}) \right] \cdot \hat{\mathbf{n}} \ \mathrm{d}\mathbf{S} = \mathbf{0}.\tag{42}$$

Placing x<sup>0</sup> 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:

$$\frac{1}{2}c(\mathbf{x}\_0) = \int\_{\mathcal{S}} \left[ c(\mathbf{x}) \nabla G(\mathbf{x}, \mathbf{x}\_0) - G(\mathbf{x}, \mathbf{x}\_0) \nabla c(\mathbf{x}) \right] \cdot \hat{\mathbf{n}} \, d\mathbf{S}.\tag{43}$$

This is a boundary integral equation (BIE) because the left hand side is the concentration cð Þ x<sup>0</sup> 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<sup>0</sup> on the surface, we obtain the boundary integral equation:

$$\frac{1}{2}c(\mathbf{x}\_0) = \sum\_{i=1}^{N\_{\text{clm}}} \left[ c\_i \left( \int\_{S\_i} \nabla G(\mathbf{x}, \mathbf{x}\_0) \cdot \hat{\mathbf{n}} \, \, d\mathbf{S} \right) - (\nabla c \cdot \hat{\mathbf{n}})\_i \left( \int\_{S\_i} G(\mathbf{x}, \mathbf{x}\_0) \, \, d\mathbf{S} \right) \right]. \tag{44}$$

Choosing x<sup>0</sup> as the midpoint x<sup>j</sup> of element j, we can write Nelm equations:

$$\frac{1}{2}c\_{j} = \sum\_{i=1}^{N\_{\text{clm}}} \left[ c\_{i} \left( \int\_{S\_{i}} \nabla G(\mathbf{x}, \mathbf{x}\_{j}) \cdot \hat{\mathbf{n}} \, \, d\mathbf{S} \right) - (\nabla c \cdot \hat{\mathbf{n}})\_{i} \left( \int\_{S\_{i}} G(\mathbf{x}, \mathbf{x}\_{j}) \, \, d\mathbf{S} \right) \right]. \tag{45}$$

The Nelm equations can be written in matrix form:

$$\left(A\_{\vec{\eta}} - \frac{1}{2}\delta\_{\vec{\eta}}\right)c\_{\vec{\jmath}} = B\_{\vec{\eta}}(\nabla c \cdot \hat{\mathbf{n}})\_{\vec{\jmath}}.\tag{46}$$

where

$$\mathcal{A}\_{\vec{\eta}} \equiv \int\_{\mathcal{S}\_i} \nabla G(\mathbf{x}, \mathbf{x}\_{\vec{\eta}}) \cdot \hat{\mathbf{n}} \, d\mathbf{S} \tag{47}$$

and

Non-Equilibrium Particle Dynamics

$$B\_{\vec{\eta}} \equiv \int\_{S\_i} G\left(\mathbf{x}, \mathbf{x}\_{\vec{\eta}}\right) \, d\mathcal{S}.\tag{48}$$

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:

> ΣijkFk � � <sup>¼</sup> <sup>∂</sup>

∂xj

<sup>G</sup>ikσij � � � <sup>∂</sup>

Gikð Þ x; x<sup>0</sup> σij � uiΣijkð Þ x; x<sup>0</sup>

∂ ∂xj

∂xj

� �nj dS, (58)

<sup>G</sup>ikð Þ <sup>x</sup>0; <sup>x</sup> <sup>σ</sup>ij <sup>þ</sup> uiΣijkð Þ <sup>x</sup>0; <sup>x</sup> � �nj dS: (59)

<sup>G</sup>ikð Þ <sup>x</sup>0; <sup>x</sup> <sup>σ</sup>ijð Þþ <sup>x</sup> <sup>Σ</sup>ijkð Þ <sup>x</sup>0; <sup>x</sup> uið Þ <sup>x</sup> � �nj dS: (60)

uiΣijk � � � �dV (57)

σijGikFk

� � (55)

Gikσij � � (56)

∂ ∂xj

uiΣijk � � <sup>þ</sup> uiFiδð Þ¼ <sup>x</sup> � <sup>x</sup><sup>0</sup> Fk

ð ∂ ∂xj

where the negative sign appears because of our convention that n^ points into V. We note that Gikð Þ¼ x; x<sup>0</sup> Gikð Þ x0; x and Σijkð Þ¼� x; x<sup>0</sup> Σijkð Þ x0; x . Additionally, the

If we choose to place x<sup>0</sup> in V, we obtain a boundary integral representation for

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

If we place x<sup>0</sup> outside V, i.e., inside the space V<sup>p</sup> enclosed by the particles, we

Finally, if we place x<sup>0</sup> on the boundary S, we obtain a boundary integral

For rigid body motion, including the 6N "primed" problems for the Lorentz reciprocal theorem, the double layer can be eliminated from the boundary integral

<sup>G</sup>ikð Þ <sup>x</sup>0; <sup>x</sup> <sup>σ</sup>ijð Þþ <sup>x</sup> <sup>Σ</sup>ijkð Þ <sup>x</sup>0; <sup>x</sup> uið Þ <sup>x</sup> � �nj dS <sup>¼</sup> <sup>0</sup>: (61)

<sup>G</sup>ikð Þ <sup>x</sup>0; <sup>x</sup> <sup>σ</sup>ijð Þþ <sup>x</sup> <sup>Σ</sup>ijkð Þ <sup>x</sup>0; <sup>x</sup> uið Þ <sup>x</sup> � �nj dS: (62)

∂ ∂xj

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

Fk ∂ ∂xj

ð V

Fk ð V

> ð V

force dipoles [59].

obtain

equation:

63

ukð Þ x<sup>0</sup> :

uiΣijkFk � � � ui

The Boundary Element Method for Fluctuating Active Colloids

We integrate both sides over the domain V:

Now we apply the divergence theorem:

ukδð Þ x � x<sup>0</sup> dV ¼ �Fk

ukδð Þ <sup>x</sup> � <sup>x</sup><sup>0</sup> dV ¼ �<sup>ð</sup>

ukð Þ¼� x<sup>0</sup>

� ð S

ukð Þ¼� x<sup>0</sup>

ð S

1 2

choice of Fk was arbitrary. We can therefore write:

ð S

ukFkδð Þ x � x<sup>0</sup> dV ¼ Fk

ð S

S

Given either a specification of either cj (Dirichlet boundary conditions) or ð Þ ∇c � n^ <sup>j</sup> (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 Bii: the evaluation point x<sup>i</sup> 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:

$$
\nabla \cdot \boldsymbol{\sigma} = -\nabla P + \mu \nabla^2 \mathbf{u} = \mathbf{0}.\tag{49}
$$

We can define a Green's function Gð Þ x; x<sup>0</sup> as the solution u xð Þ� Gð Þ� x; x0 F to the Stokes equation with a body point force F located at x0:

$$-\nabla P + \mu \nabla^2 \mathbf{u} + \mathbf{F} \delta(\mathbf{x} - \mathbf{x}\_0) = \mathbf{0},\tag{50}$$

or

$$\nabla \cdot \boldsymbol{\sigma} = -\mathbf{F} \delta(\mathbf{x} - \mathbf{x}\_0). \tag{51}$$

It can be shown that the Green's function is

$$\mathcal{G}\_{\vec{\eta}}(\mathbf{x}, \mathbf{x}\_0) = \frac{1}{8\pi\mu r} \left(\delta\_{\vec{\eta}} + \frac{\tilde{\mathbf{x}}\_i \tilde{\mathbf{x}}\_j}{r^2}\right),\tag{52}$$

where r � ∣x � x0∣ and where x~<sup>j</sup> � xj � x0,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<sup>0</sup> F, where

$$\mathcal{P}\_j = \frac{\ddot{\mathfrak{x}}\_j}{4\pi r^3}.\tag{53}$$

The stress in the fluid is given by σ � Σð Þ� x; x<sup>0</sup> F, where

$$
\Sigma\_{ijk} = -\Im \frac{\tilde{\varkappa}\_i \tilde{\varkappa}\_j \tilde{\varkappa}\_k}{4\pi r^5}. \tag{54}
$$

Now we wish to apply Eq. 30. We specify the "primed" fields u<sup>0</sup> and σ<sup>0</sup> as the fields due to a point force F at x<sup>0</sup> in unbounded fluid. For the "unprimed" fields u Bij � ð Si

2.5 Boundary integral formulation of the Stokes equation

the Stokes equation with a body point force F located at x0:

It can be shown that the Green's function is

�∇<sup>P</sup> <sup>þ</sup> <sup>μ</sup>∇<sup>2</sup>

Gijð Þ¼ x; x<sup>0</sup>

The stress in the fluid is given by σ � Σð Þ� x; x<sup>0</sup> F, where

solved numerically with standard methods.

Non-Equilibrium Particle Dynamics

in Ref. [51].

[52, 54, 61].

or

62

Pð Þ� x Pð Þ� x; x<sup>0</sup> F, where

Stokes equation is:

G x; x<sup>j</sup>

Given either a specification of either cj (Dirichlet boundary conditions) or ð Þ ∇c � n^ <sup>j</sup> (Neumann boundary conditions), the algebraic system in Eq. 46 can be

A certain difficulty becomes apparent when we consider the element Bii: the evaluation point x<sup>i</sup> 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

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 ε

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

We can define a Green's function Gð Þ x; x<sup>0</sup> as the solution u xð Þ� Gð Þ� x; x0 F to

1 8πμr

<sup>P</sup><sup>j</sup> <sup>¼</sup> <sup>x</sup>~<sup>j</sup>

Σijk ¼ �3

where r � ∣x � x0∣ and where x~<sup>j</sup> � xj � x0,j. Eq. 52 is commonly called the Oseen tensor or "Stokeslet". The fluid pressure in response to the point force is given by

> x~i x~jx~<sup>k</sup>

Now we wish to apply Eq. 30. We specify the "primed" fields u<sup>0</sup> and σ<sup>0</sup> as the fields due to a point force F at x<sup>0</sup> in unbounded fluid. For the "unprimed" fields u

δij þ

x~ix~<sup>j</sup> r2 � �

<sup>∇</sup> � <sup>σ</sup> ¼ �∇<sup>P</sup> <sup>þ</sup> <sup>μ</sup>∇<sup>2</sup>

� � dS: (48)

u ¼ 0: (49)

, (52)

u þ Fδð Þ¼ x � x<sup>0</sup> 0, (50)

<sup>4</sup>πr<sup>3</sup> : (53)

<sup>4</sup>πr<sup>5</sup> : (54)

∇ � σ ¼ �Fδð Þ x � x<sup>0</sup> : (51)

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:

$$\frac{\partial}{\partial \mathbf{x}\_j} \left( u\_i \Sigma\_{ijk} F\_k \right) - u\_i \frac{\partial}{\partial \mathbf{x}\_j} \left( \Sigma\_{ijk} F\_k \right) = \frac{\partial}{\partial \mathbf{x}\_j} \left( \sigma\_{\vec{\eta}} \mathcal{G}\_{ik} F\_k \right) \tag{55}$$

$$F\_k \frac{\partial}{\partial \mathbf{x}\_j} \left( u\_i \Sigma\_{\vec{\eta}k} \right) + u\_i F\_i \delta(\mathbf{x} - \mathbf{x}\_0) = F\_k \frac{\partial}{\partial \mathbf{x}\_j} \left( \mathcal{G}\_{ik} \sigma\_{\vec{\eta}} \right) \tag{56}$$

We integrate both sides over the domain V:

$$\int\_{\mathcal{V}} \mu\_k F\_k \delta(\mathbf{x} - \mathbf{x}\_0) dV = F\_k \int \left[ \frac{\partial}{\partial \mathbf{x}\_j} \left( \mathcal{G}\_{ik} \sigma\_{ij} \right) - \frac{\partial}{\partial \mathbf{x}\_j} \left( \mu\_i \Sigma\_{ijk} \right) \right] dV \tag{57}$$

Now we apply the divergence theorem:

$$F\_k \int\_{\mathcal{V}} u\_k \delta(\mathbf{x} - \mathbf{x}\_0) dV = -F\_k \int\_{\mathcal{S}} \left[ \mathcal{G}\_{ik}(\mathbf{x}, \mathbf{x}\_0) \sigma\_{\vec{\eta}} - u\_i \Sigma\_{\vec{\eta}k}(\mathbf{x}, \mathbf{x}\_0) \right] \eta\_j dS,\tag{58}$$

where the negative sign appears because of our convention that n^ points into V. We note that Gikð Þ¼ x; x<sup>0</sup> Gikð Þ x0; x and Σijkð Þ¼� x; x<sup>0</sup> Σijkð Þ x0; x . Additionally, the choice of Fk was arbitrary. We can therefore write:

$$\int\_{\mathcal{V}} \boldsymbol{\mu}\_{k} \delta(\mathbf{x} - \mathbf{x}\_{0}) dV = -\int\_{\mathcal{S}} \left[ \mathcal{G}\_{ik}(\mathbf{x}\_{0}, \mathbf{x}) \, \boldsymbol{\sigma}\_{ij} + \boldsymbol{\mu}\_{i} \boldsymbol{\Sigma}\_{j\bar{k}}(\mathbf{x}\_{0}, \mathbf{x}) \right] \boldsymbol{\eta}\_{j} d\mathbf{S}.\tag{59}$$

If we choose to place x<sup>0</sup> in V, we obtain a boundary integral representation for ukð Þ x<sup>0</sup> :

$$u\_k(\mathbf{x}\_0) = -\int\_{\mathcal{S}} \left[ \mathcal{G}\_{ik}(\mathbf{x}\_0, \mathbf{x}) \, \sigma\_{\vec{\eta}}(\mathbf{x}) + \Sigma\_{\vec{\eta}k}(\mathbf{x}\_0, \mathbf{x}) u\_i(\mathbf{x}) \right] n\_{\vec{\eta}} d\mathcal{S}.\tag{60}$$

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<sup>0</sup> outside V, i.e., inside the space V<sup>p</sup> enclosed by the particles, we obtain

$$-\int\_{\mathcal{S}} \left[ \mathcal{G}\_{ik}(\mathbf{x}\_0, \mathbf{x}) \, \sigma\_{ij}(\mathbf{x}) + \Sigma\_{ijk}(\mathbf{x}\_0, \mathbf{x}) \, u\_i(\mathbf{x}) \right] n\_j dS = \mathbf{0}.\tag{61}$$

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

$$\frac{1}{2}u\_k(\mathbf{x}\_0) = -\int\_S \left[\mathcal{G}\_{ik}(\mathbf{x}\_0, \mathbf{x})\,\sigma\_{\vec{\eta}}(\mathbf{x}) + \Sigma\_{\vec{\eta}\vec{k}}(\mathbf{x}\_0, \mathbf{x})u\_i(\mathbf{x})\right] \eta\_{\vec{\eta}}d\mathbf{S}.\tag{62}$$

For rigid body motion, including the 6N "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 Vp. Within Vp, the flow field u is simply the flow field uRBM for rigid body motion, since it must obey no-slip on the surface S. Now we apply Eq. 59 for the field uRBM inside V<sup>p</sup> and x<sup>0</sup> ∈V, noting that we must use a normal n^<sup>0</sup> ¼ �n^ pointing into Vp:

$$-\int\_{\mathcal{S}} \left[ \mathcal{G}\_{ik}(\mathbf{x}\_0, \mathbf{x}) \, \sigma\_{\vec{\eta}}^{\text{RBM}}(\mathbf{x}) + \Sigma\_{\vec{\eta}\vec{k}}(\mathbf{x}\_0, \mathbf{x}) u\_i^{\text{RBM}}(\mathbf{x}) \right] n\_j' d\mathbf{S} = \mathbf{0}.\tag{63}$$

For rigid body motion, there is no shear stress and the pressure is uniform, i.e., σRBM ij ð Þ x n^<sup>0</sup> ¼ �p0n^<sup>0</sup> . The first term is simply the integral of Gikð Þ� x0; x n^<sup>0</sup> over the surface S for x<sup>0</sup> ∈V, which vanishes identically by incompressibility. This leaves:

$$-\int\_{S} \left[\Sigma\_{ijk}(\mathbf{x}\_{0}, \mathbf{x}) u\_{i}^{\text{RBM}}(\mathbf{x})\right] n\_{j}^{\prime} d\mathbf{S} = \mathbf{0}.\tag{64}$$

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

So far, we have considered the deterministic contributions to the 6N 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

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

confined geometry can often be constructed by the method of images.

The Boundary Element Method for Fluctuating Active Colloids

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

long-time behavior of an active colloid, it is absolutely essential.

Langevin equation for the generalized, 6N-component coordinate q:

dt <sup>¼</sup> <sup>V</sup> <sup>þ</sup> kBTð Þþ <sup>∇</sup> �<sup>M</sup> ffiffiffiffiffiffiffiffiffiffi

Wiener processes. Discretizing time in steps of Δt, one can write a generalized

<sup>Δ</sup><sup>q</sup> <sup>¼</sup> <sup>V</sup>Δ<sup>t</sup> <sup>þ</sup> kBTð Þ <sup>∇</sup> �<sup>M</sup> <sup>Δ</sup><sup>t</sup> <sup>þ</sup> ffiffiffiffiffiffiffiffiffiffiffi

stochastic drift term ð Þ ∇ �M is a consequence of having a configuration-dependent

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

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

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

displacement Δq as the following Euler-Maruyama equation [34, 64]:

where Δw is a stochastic variable with h i Δw ¼ 0 and Δwi Δwj

mobility tensor in the framework of the Itô interpretation.

of the chemical field cð Þ x in the fluid domain V.

where V is the deterministic contribution of activity to the generalized velocity ð Þ <sup>U</sup>α; <sup>Ω</sup><sup>α</sup> <sup>T</sup>, i.e., the solution to the problem outlined above; <sup>M</sup> is the grand mobility

; <sup>B</sup> satisfies <sup>B</sup> � <sup>B</sup><sup>T</sup> <sup>¼</sup> <sup>M</sup>; and <sup>W</sup> is a collection of independent

<sup>2</sup>kBT <sup>p</sup> <sup>B</sup> � <sup>W</sup>, (67)

<sup>2</sup>kBT <sup>p</sup> <sup>B</sup> � <sup>Δ</sup>w, (68)

� � <sup>¼</sup> <sup>δ</sup>ijΔt. The

dq

2.7 Thermal fluctuations

matrix <sup>M</sup> <sup>¼</sup> <sup>R</sup>�<sup>1</sup>

reference frames [67].

field on the boundary S [71].

65

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

$$\mu\_k(\mathbf{x}\_0) = -\int\_{\mathcal{S}} \left[ \mathcal{G}\_{ik}(\mathbf{x}\_0, \mathbf{x}) \sigma\_{\vec{\eta}}(\mathbf{x}) \right] \eta\_j d\mathbf{S}, \ \mathbf{x}\_0 \in \mathcal{V}. \tag{65}$$

In order to obtain a single-layer boundary integral equation for x<sup>0</sup> ∈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:

$$\mu\_k(\mathbf{x}\_0) = -\int\_{\mathcal{S}} \left[ \mathcal{G}\_{ik}(\mathbf{x}\_0, \mathbf{x}) \sigma\_{\hat{\eta}}(\mathbf{x}) \right] \eta\_{\hat{\eta}} d\mathbf{S}, \ \mathbf{x}\_0 \in \mathcal{S}. \tag{66}$$

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.
