3.1.1. Standard multipole expansion for Laplace problem

The exact solution of the Laplace problems can in fact be recovered from (3.1), when the function $g(\boldsymbol {r},t)$ is taken as a (possibly infinite) set of singularities centred on each particle (Saffman Reference Saffman1973),

(3.2)\begin{equation} g(\boldsymbol{r},t) = \sum_{n=1}^N [q^M_n \delta(\boldsymbol{r}_n) {-} \boldsymbol{q}^D_n \boldsymbol{\cdot}\boldsymbol{\nabla}\delta(\boldsymbol{r}_n) + \cdots], \end{equation}

where $\delta (\boldsymbol {r}_n)$ is the Dirac delta distribution, and ($q_n^M$, $\boldsymbol {q}_n^D,\ldots$) are the intensities of the singularities associated with particle $n$, and are constant tensors of increasing order. Note that $\boldsymbol {\nabla }$ denotes here the gradient with respect to the observation position $\boldsymbol {r}$ and $\boldsymbol {r}_n=\boldsymbol {r}-\boldsymbol {Y}_n$. This equation can be solved explicitly for the concentration field $c$ as a multipole expansion for each particle in terms of source monopoles, dipoles, etc$\ldots$

(3.3)\begin{equation} c(\boldsymbol{r},t) = \sum_{n=1}^N [q^M_n G^M(\boldsymbol{r}_n) + \boldsymbol{q}^D_n \cdot \boldsymbol{G}^D(\boldsymbol{r}_n) + \cdots], \end{equation}

where $G^M$ and $\boldsymbol {G}^D$ are the monopole and dipole Green's functions and satisfy

(3.4a,b)\begin{equation} \nabla^2 G^M={-}\delta(\boldsymbol{r}_n),\quad \nabla^2 \boldsymbol{G}^D= \boldsymbol{\nabla}\delta(\boldsymbol{r}_n), \end{equation}

together with appropriate decay or boundary conditions on the domain's outer boundary. For unbounded domains with decaying conditions in the far field, the singular monopole and dipole Green's functions are simply

(3.5a,b)\begin{equation} G^M(\boldsymbol{r}_n) = \frac{1}{4{\rm \pi} r_n}\quad \textrm{and}\quad \boldsymbol{G}^D(\boldsymbol{r}_n)={-}\boldsymbol{\nabla} G^M =\frac{\boldsymbol{r}_n}{4{\rm \pi} r_n^3}\cdot. \end{equation}

The concentration distributions associated with these singular Green's functions are displayed in figure 3. Higher-order derivatives of $G^M(\boldsymbol {r})$, (3.5a,b), are also solutions of Laplace's equation, leading to singularities of increasing order (quadrupole, octupole,…).

Figure 3.

Singular (dotted lines, (3.5a,b)) and regularized (solid lines, (3.10)–(3.11)) concentration distributions along the axial polar direction associated with the Greens’ functions for the Laplace equation for (a) monopole terms and (b) dipole terms. The line $r/a=1$ represents the particle surface.

3.1.2. Truncated regularized multipole expansion

The previous approach, based on an infinite set of singular sources, is known as the standard multipole expansion of the Laplace problem. Although satisfying from a theoretical point of view, since it is able to recover an accurate representation of the analytical solution outside the particles for a large enough number of singular multipoles, it is not well suited for a versatile numerical implementation because of (i) the singular behaviour of the forcing terms in the modified Laplace equation, (3.1), and (ii) the a priori infinite set of singularities required for each particle.

To avoid the latter issue, the infinite expansion is truncated here after the first two terms, thus retaining the monopole and dipole contributions only. Physically, this amounts to retaining the two leading physical effects of the particle on the concentration field, i.e. a net emission with a front–back asymmetric distribution. In order to overcome the former problem, the standard FCM replaces the singular Dirac distributions $\delta (\boldsymbol {r})$ by regular Gaussian spreading functions $\varDelta (\boldsymbol {r})$:

(3.6)\begin{equation} \triangle(\boldsymbol{r}) = (2{\rm \pi}\sigma^2)^{{-}3/2} \exp\left(-\frac{r^2}{2\sigma^2}\right), \end{equation}

where $\sigma$ denotes the finite-size support of this envelop and acts as a smoothing parameter of the method, thus eliminating the singular behaviour of the delta distribution $\delta (\boldsymbol {r})$ near the origin, thereby allowing for a more accurate numerical treatment. The original singular distribution is recovered when $\sigma \ll r$, i.e. the solution of the regularized problem is an accurate representation of the true solution away from the particle. This approach using regular distributions allows for a more versatile and robust numerical solution of the physical equations than their singular counterparts (Maxey & Patel Reference Maxey and Patel2001; Lomholt & Maxey Reference Lomholt and Maxey2003).

Combining these two approximations, we therefore consider a truncated regularized expansion including only the monopole and the dipole terms as

(3.7)\begin{equation} g(\boldsymbol{r},t) = \sum_{n=1}^N [q^M_n \varDelta^M({r}_n) {-} \boldsymbol{q}^D_n \boldsymbol{\cdot}\boldsymbol{\nabla} \varDelta^D({r}_n)], \end{equation}

with the Gaussian spreading operators $\varDelta ^M$ and $\varDelta ^D$ defined as

(3.8a,b)\begin{equation} \varDelta^M(r) = (2{\rm \pi}\sigma_M^2)^{{-}3/2} \exp\left(-\frac{r^2}{2 \sigma_M^2}\right),\quad \varDelta^D(r) = (2{\rm \pi}\sigma_D^2)^{{-}3/2} \exp\left(-\frac{r^2}{2 \sigma_D^2}\right), \end{equation}

where $M$ and $D$ once again denote monopole and dipole, and $\sigma _M$ and $\sigma _D$ are the finite support of each regularized distribution and are free numerical parameters of the method that need to be calibrated. Note that, in general, these do not need to be identical (Lomholt & Maxey Reference Lomholt and Maxey2003).

The corresponding truncated regularized solution for $c$ is then finally obtained as

(3.9)\begin{equation} c(\boldsymbol{r},t) = \sum_{n=1}^N [q^M_n G^M(\boldsymbol{r}_n) + \boldsymbol{q}^D_n \cdot \boldsymbol{G}^D(\boldsymbol{r}_n)], \end{equation}

with the regularized monopole and dipole Green's functions

(3.10)$$\begin{gather} G^M(\boldsymbol{r}) = \frac{1}{4{\rm \pi} r} \mathrm{erf}\left( \frac{r}{{\sigma_M}\sqrt{2}} \right), \end{gather}$$

(3.11)$$\begin{gather}\boldsymbol{G}^D(\boldsymbol{r})= \frac{\boldsymbol{r}}{4{\rm \pi} r^3} \left[\mathrm{erf}\left(\frac{r}{\sigma_D \sqrt{2}}\right) -\sqrt{\frac{2}{\rm \pi}} \left(\frac{r}{\sigma_D}\right) \exp\left(-\frac{r^2}{2\sigma_D^2}\right)\right]. \end{gather}$$

These clearly match the behaviour of their singular counterpart, (3.5a,b), when $r$ is greater than a few $\sigma _M$ or $\sigma _D$, respectively, while still maintaining finite values within the particle (figure 3), e.g. $\boldsymbol {G}^D(\boldsymbol {r}=\boldsymbol {0})=\boldsymbol {0}$.

3.1.3. Finding the intensity of the singularities

Up to this point, no information was implemented regarding the surface boundary conditions on $c$ in (2.1). We now present how to determine the intensities of the monopole and dipole distributions associated with each particle, $q^M_n$ and $\boldsymbol {q}^D_n$, so as to satisfy a weak form of the Neuman boundary condition, (2.1), i.e. its first two moments over the particle's surface. Using the multipole expansion of the fundamental integral representation of the concentration (see Appendix A), the monopole and dipole intensities of particle $n$, $q^M_n$ and $\boldsymbol {q}^{D}_n$, are obtained as (Yan & Brady Reference Yan and Brady2016)

(3.12a,b)\begin{equation} q^M_n = \int_{S_n} \alpha_n \mathrm{d}S,\quad \boldsymbol{q}^{D}_n = a \int_{S_n} \alpha_n \boldsymbol{n}\,\mathrm{d}S + 4{\rm \pi} a^2 \langle c \boldsymbol{n} \rangle_n, \end{equation}

where the second term in $\boldsymbol {q}_n^D$ is proportional to the concentration polarity at the surface of particle $n$, i.e. its first moment $\langle c \boldsymbol {n} \rangle _n$, and is defined using the surface average operator $\langle \cdot \rangle _n$ over the surface of particle $n$. Note that the activity distribution at the particle's surface is known, and thus (2.9) explicitly provides the monopole intensity and the first term in the dipole intensity. The second contribution to the latter requires, however, knowledge of the solution on the particle's surface – which is not explicitly represented in the present FCM approach. This term therefore requires to be solved for as part of the general problem. In the previous equation, it should be noted that the dimensionless particle radius is $a=1$, but will be kept in the equations to emphasize the relative scaling of the numerical spreading envelopes (e.g. $\sigma _M$ and $\sigma _D$) with respect to the particle size.

Here, we use an iterative approach to solve this linear joint problem for the dipole intensity and concentration field, solving alternatively (3.7) and (3.12a,b) until convergence is reached, as defined by the following criterion between two successive iterations:

(3.13)\begin{equation} \left\| \frac{ \langle c\boldsymbol{n}\rangle^{k+1}- \langle c\boldsymbol{n}\rangle^{k}}{\langle c\boldsymbol{n}\rangle^{k+1}}\right\|_{\infty} < \epsilon, \end{equation}

where $\langle c\boldsymbol {n}\rangle ^{k}$ is the vector collecting the polarities of the $N$ particles at iteration $k$. For the results presented in this work, we set the tolerance to $\epsilon =10^{-10}$ in our calculations.

3.1.4. Regularized moments of the concentration distribution

Finding the dipole intensity, $\boldsymbol {q}^{D}_n$, requires computing the polarity $\langle c \boldsymbol {n} \rangle _n$, which is in principle defined at the particle's surface. To follow the spirit of FCM, and allow for efficient numerical treatment, this surface projection is replaced by a weighted projection over the entire volume $V_F$

(3.14)\begin{equation} \langle c \boldsymbol{n} \rangle_n =\frac{1}{4{\rm \pi} a^2} \int_{S_n} c \boldsymbol{n} \mathrm{d}S \quad\longrightarrow\quad \{ c \boldsymbol{n} \}_n = \int_{V_F} c \boldsymbol{n}_n \varDelta^P(\boldsymbol{r}_n)\,\mathrm{d}V, \end{equation}

with $\boldsymbol {n}_n$ now defined as $\boldsymbol {n}_n={\boldsymbol {r}_n}/{r_n}$, and the regular averaging kernel $\varDelta ^P$ for the polarity as

(3.15)\begin{equation} \varDelta^P(\boldsymbol{r}) = \frac{r}{8{\rm \pi}\sigma_P^4} \exp\left(-\frac{r^2}{2 \sigma_P^2}\right). \end{equation}

Beyond its importance for determining the dipole intensity associated with a given particle, we will later show that the polarity of the concentration at the surface of particle $n$ is directly related to its self-induced phoretic velocity, (2.7a,b), and that, similarly, the self-induced hydrodynamic stresslet signature of the particle is in general associated with the first two moments of the surface concentration. Similarly to the polarity, the second surface moment, $\langle c(\boldsymbol {n}\boldsymbol {n}-{\boldsymbol{\mathsf{I}}}/3) \rangle _n$ will be replaced in our implementation by a weighted volume projection $\{c(\boldsymbol {n}\boldsymbol {n}-{\boldsymbol{\mathsf{I}}}/3)\}_n$

(3.16)\begin{align} \langle c(\boldsymbol{n}\boldsymbol{n}-{\boldsymbol{\mathsf{I}}}/3) \rangle_n &= \frac{1}{4{\rm \pi} a^2}\int_{S_n}c\left(\boldsymbol{n}\boldsymbol{n}-\frac{{\boldsymbol{\mathsf{I}}}}{3}\right)\mathrm{d}S \rightarrow \{ c(\boldsymbol{n}\boldsymbol{n}-{\boldsymbol{\mathsf{I}}}/3) \}_n \nonumber\\ &= \int_{V_F} c \left(\boldsymbol{n}_n\boldsymbol{n}_n-\frac{{\boldsymbol{\mathsf{I}}}}{3}\right) \varDelta^S(\boldsymbol{r}_n)\,\mathrm{d}V, \end{align}

where the projection kernel for the second moment of concentration, $\varDelta ^S$, is defined as

(3.17)\begin{equation} \varDelta^S(\boldsymbol{r})= \frac{r^2}{3(2{\rm \pi})^{{3}/{2}} \sigma_S^5} \exp\left(-\frac{r^2}{2 \sigma_S^2}\right). \end{equation}

The envelopes $\sigma _P$ and $\sigma _S$ are free parameters in the method that need to be calibrated. In our reactive FCM formulation, we use modified forms of the Gaussian operator $\varDelta$ as projection operators, (3.15) and (3.17), in order to ensure a fast numerical convergence of the integration for the first and second moment calculations, (3.14) and (3.16) respectively. The integrals over the entire volume $V_F$ of these averaging functions is still equal to one, and their weight is shifted from the particle centre toward the particle surface (figure 4), which is both numerically more accurate and more intuitive physically as these operators are used to obtain the properties of the particle on its surface.

Figure 4.

Averaging envelopes for the first and second moments of concentration, $\varDelta ^P$ (solid, (3.15)) and $\varDelta ^S$ (dashed, (3.17)) respectively. The numerical values for $\sigma _P$ and $\sigma _S$ are set from (3.21a,b) and (3.24).

3.1.5. Calibrating the spreading/averaging envelopes

Our method relies on four numerical parameters ($\sigma _M$, $\sigma _D$, $\sigma _P$, $\sigma _S$) that we choose to calibrate so as to ensure that several key results in reference configurations are obtained exactly. In particular, to properly account for the phoretic drift induced by the other particles, we ensure that the polarity $\langle c\boldsymbol {n}\rangle$ of an isolated particle placed in an externally imposed uniform gradient of concentration can be exactly recovered using the regular representation and averaging operators. A similar approach is then followed for the particle's second moment of concentration $\langle c(\boldsymbol {n}\boldsymbol {n}-{\boldsymbol{\mathsf{I}}}/3)\rangle$ in a quadratic externally imposed field.

Isolated passive particle in an external linear field – We first consider a single particle placed at the origin in an externally imposed linear concentration field so that for $r\gg a$, $c\approx c_E$ with

(3.18)\begin{equation} c_E = \boldsymbol{L}_E \cdot \boldsymbol{r}, \end{equation}

where $\boldsymbol {L}_E$ is the externally imposed uniform gradient. For a passive particle (i.e. $\alpha =0$), satisfying the boundary condition, (2.1), at the surface of the particle imposes that the exact concentration distribution around the particle is $c=c_E+c_I^o$, with ${c_I^o(\boldsymbol {r})=a^3\boldsymbol {L}_E\cdot \boldsymbol {r}/(2r^3)}$ a singular dipole-induced field. The polarity of the external and induced parts, $c_E$ and $c_I$, can be obtained analytically as

(3.19a,b)\begin{equation} \langle c_E \boldsymbol{n} \rangle = \frac{a}{3}\boldsymbol{L}_E,\quad \langle c_I^o \boldsymbol{n} \rangle = \frac{a}{6}\boldsymbol{L}_E. \end{equation}

Following the framework presented above, the regularized solution can be written $c=c_E+c_I^r$ with $c_i^r$ a regularized dipole, and the corresponding regularized-volume moments based on (3.14) are obtained using (3.11), as

(3.20a,b)\begin{equation} \{ c_E \boldsymbol{n} \} = \sqrt{\frac{\rm \pi}{8}} \sigma_P \boldsymbol{L}_E,\quad \{ c_I^r \boldsymbol{n} \} = \frac{a^3\sigma_P}{12(\sigma_D^2+\sigma_P^2)^{{3}/{2}}} \boldsymbol{L}_E. \end{equation}

Identification of the regularized result (3.20a,b) with the true solution (3.19a,b), determines $\sigma _P$ and $\sigma _D$ uniquely as

(3.21a,b)\begin{equation} \frac{\sigma_P}{a}=\frac{1}{3}\sqrt{\frac{8}{\rm \pi}} \approx 0.5319,\quad \frac{\sigma_D}{a} = \sqrt{\left(\frac{\sigma_P}{2a}\right)^{2/3} - \left(\frac{\sigma_P}{a}\right)^{2}} \approx 0.3614. \end{equation}

Isolated passive particle in an external quadratic field – Similarly, in an external quadratic field $c_E$ of the form

(3.22)\begin{equation} c_E(\boldsymbol{r})= \boldsymbol{r} \cdot {\boldsymbol{\mathsf{Q}}}_E \cdot \boldsymbol{r}, \end{equation}

with ${\boldsymbol{\mathsf{Q}}}_E$ a second-order symmetric and traceless tensor, the concentration distribution around a passive particle ($\alpha =0$) takes the form $c=c_E+c_I^o$ with $c_I^o(\boldsymbol {r})$ an induced singular quadrupole. The exact and regularized second moments of the external field $c_E$ at the particle surface are equal to

(3.23a,b)\begin{equation} \langle c_E (\boldsymbol{n}\boldsymbol{n}-{\boldsymbol{\mathsf{I}}}/3) \rangle = \frac{2a^2}{15}{\boldsymbol{\mathsf{Q}}}_E,\quad \{c_E (\boldsymbol{n}\boldsymbol{n}-{\boldsymbol{\mathsf{I}}}/3) \} = \frac{2\sigma_S^2}{3}{\boldsymbol{\mathsf{Q}}}_E. \end{equation}

Identifying both results determines the size of the averaging envelope for the second moment uniquely, as

(3.24)\begin{equation} \frac{\sigma_S}{a}=\sqrt{\frac{1}{5}} \approx 0.4472. \end{equation}

Note that we do not enforce here a constraint on the representation of the second moment of the induced field $c_I$, since the particles’ representations do not include a regularized quadrupole in our method.

The value $\sigma _M$ remains as a free parameter at this point and cannot be calibrated with a similar approach. In the following, in order to minimize the number of distinct numerical parameters and to minimize the departure of the regularized solution from its singular counterpart, we set its value equal to the smallest envelope size, namely $\sigma _M=\sigma _D$. These specific values of the parameters were used in figures 3 and 4.

3.2.1. FCM for passive suspensions

With hydrodynamic FCM, the effect of the particles on the fluid is accounted for through a forcing term $\boldsymbol {f}$ applied to the dimensionless Stokes equations

(3.25)\begin{equation} \boldsymbol{\nabla} p - \nabla^2 \boldsymbol{u}= \boldsymbol{f}(\boldsymbol{r},t) \quad \mathrm{in}\ V_F. \end{equation}

As for reactive FCM, this forcing arises from a truncated regularized multipolar expansion up to the dipole level

(3.26)\begin{equation} \boldsymbol{f}(\boldsymbol{r},t) = \sum_{n=1}^N [\boldsymbol{F}_n \varDelta (r_n) + {\boldsymbol{\mathsf{D}}}_n \boldsymbol{\cdot}\boldsymbol{\nabla} \varDelta^*(r_n)], \end{equation}

where the spreading envelopes are defined by

(3.27a,b)\begin{equation} \varDelta(r) = (2{\rm \pi}\sigma^2)^{{-}3/2} \exp\left(-\frac{r^2}{2 \sigma^2}\right),\quad \varDelta^*(r) = (2{\rm \pi}\sigma_*^2)^{{-}3/2} \exp\left(-\frac{r^2}{2 \sigma_*^2}\right). \end{equation}

Here, $\boldsymbol {F}_n$ and ${\boldsymbol{\mathsf{D}}}_n$ are the force monopole and dipole applied to particle $n$. The force dipole can be split into a symmetric part, the stresslet ${\boldsymbol{\mathsf{S}}}$, and an antisymmetric one related to the external torque $\boldsymbol {T}$

(3.28)\begin{equation} {\boldsymbol{\mathsf{D}}}_n = {\boldsymbol{\mathsf{S}}}_n + \tfrac{1}{2} \boldsymbol{\epsilon} \cdot \boldsymbol{T}_n, \end{equation}

with $\boldsymbol{\epsilon}$ the third-order permutation tensor. The corresponding regularized solution for the fluid velocity $\boldsymbol {u}$ is then obtained as

(3.29)\begin{equation} \boldsymbol{u} = \boldsymbol{u}(\boldsymbol{r}) = \sum_{n=1}^{N} [\boldsymbol{F}_n \cdot {\boldsymbol{\mathsf{J}}} (\boldsymbol{r}_n)+ {\boldsymbol{\mathsf{D}}}_n : {\boldsymbol{\mathsf{R}}}^*(\boldsymbol{r}_n)]. \end{equation}

For unbounded domains with vanishing perturbations in the far field (i.e. $\|\boldsymbol {u}\|\rightarrow 0$ when $r \rightarrow \infty$), the regularized Green's function ${\boldsymbol{\mathsf{J}}}(\boldsymbol {r})$ reads

(3.30)\begin{equation} {\boldsymbol{\mathsf{J}}}(\boldsymbol{r}) = \frac{1}{8 {\rm \pi}r} \left(A(r) {\boldsymbol{\mathsf{I}}} + B(r) \frac{\boldsymbol{r}\boldsymbol{r}}{r^2}\right), \end{equation}

with

(3.31)$$\begin{gather} A(r)= \left(1+\frac{\sigma^2}{r^2}\right) \mathrm{erf}\left(\frac{r}{\sigma \sqrt{2}}\right) - \frac{\sigma}{r} \sqrt{\frac{2}{\rm \pi}} \exp\left(-\frac{r^2}{2\sigma^2}\right), \end{gather}$$

(3.32)$$\begin{gather}B(r)= \left(1-\frac{3\sigma^2}{r^2}\right) \mathrm{erf}\left(\frac{r}{\sigma \sqrt{2}}\right) + \frac{3\sigma}{r}\sqrt{\frac{2}{\rm \pi}} \exp\left(-\frac{r^2}{2\sigma^2}\right), \end{gather}$$

and ${\boldsymbol{\mathsf{R}}}^*=\boldsymbol {\nabla } {\boldsymbol{\mathsf{J}}}^*$ is the FCM dipole Green's function evaluated with the parameter $\sigma _*$.

The particle's translational and angular velocities, $\boldsymbol {U}_n$ and $\boldsymbol {\varOmega }_n$, are obtained from a volume-weighted average of the local fluid velocity and vorticity

(3.33a,b)\begin{equation} \boldsymbol{U}_n = \int_{V_F} \boldsymbol{u} \varDelta(\boldsymbol{r}_n)\,\mathrm{d}V,\quad \boldsymbol{\varOmega}_n = \tfrac{1}{2}\int_{V_F} [\boldsymbol{\nabla} \times \boldsymbol{u}] \varDelta^*(\boldsymbol{r}_n)\,\mathrm{d}V. \end{equation}

The Gaussian parameters, $\sigma$ and $\sigma ^*$, are calibrated to recover the correct Stokes drag, ${\boldsymbol {F}=6{\rm \pi} a \mu \boldsymbol {U}}$, and viscous torque, $\boldsymbol {T}=8{\rm \pi} a^3 \mu \boldsymbol {\varOmega }$, of an isolated particle (Maxey & Patel Reference Maxey and Patel2001; Lomholt & Maxey Reference Lomholt and Maxey2003), leading to

(3.34a,b)\begin{equation} \frac{\sigma}{a} =\frac{1}{\sqrt{\rm \pi}} \approx 0.5641,\quad \frac{\sigma_*}{a}=\frac{1}{(6\sqrt{\rm \pi})^{1/3}} \approx 0.4547. \end{equation}

The rigidity of the particle is similarly weakly enforced by imposing that the volume-averaged strain rate ${\boldsymbol{\mathsf{E}}}_n$ over the envelope of particle $n$ vanishes

(3.35)\begin{equation} {\boldsymbol{\mathsf{E}}}_n = \tfrac{1}{2} \int_{V_F} [\boldsymbol{\nabla} \boldsymbol{u} + (\boldsymbol{\nabla} \boldsymbol{u})^\mathrm{T}] \varDelta^*(\boldsymbol{r}_n)\,\mathrm{d}V = \boldsymbol{0}, \end{equation}

which determines the stresslet ${\boldsymbol{\mathsf{S}}}_n$ induced by particle $n$. Note that, unlike forces and torques which are typically set by external or inter-particle potentials, the stresslets result from the constraint on the flow given by (3.35) and, consequently, need to be solved for as part of the general flow problem. The resulting linear system for the unknown stresslet coefficients is solved directly or iteratively, with the conjugate gradients method, depending on the number of particles considered (Lomholt & Maxey Reference Lomholt and Maxey2003; Yeo & Maxey Reference Yeo and Maxey2010). In the following, we consider pairs of particles (see § 4) and therefore use direct inversion.

Note that the averaging envelopes used to recover the translational and rotational velocities, $\triangle _n$ and $\triangle ^*_n$, are exactly the same as the spreading operators in (3.26), all of them Gaussian functions. As a result, the spreading and averaging operators are adjoints to one another. Also note that only two envelope lengths are required for the hydrodynamic problem: $\sigma$ and $\sigma _*$. In contrast, the new reactive FCM extension presented in § 3.1 uses spreading and averaging operators that are not adjoint. To recover the first (3.14) and second (3.16) moments of concentration we have two non-Gaussian averaging envelopes ($\varDelta ^P$ and $\varDelta ^S$), that differ from the Gaussian spreading envelopes ($\varDelta ^M$ and $\varDelta ^D$) in (3.7). While having adjoint operators is crucial in hydrodynamic FCM to satisfy the fluctuation–dissipation balance, the lack of adjoint properties for the Laplace problem does not raise any issue in the deterministic setting.

3.2.2. Active hydrodynamic FCM

In recent years, FCM has been extended to handle suspensions of active particles, such as microswimmers. In addition to undergoing rigid body motion in the absence of applied forces or torques, active and self-propelled particles are also characterized by the flows they generate. These flows can be incorporated into FCM by adding an appropriate set of regularized multipoles to the Stokes equations. This problem was solved previously for the classical squirmer model (Delmotte et al. Reference Delmotte, Keaveny, Plouraboué and Climent2015), a spherical self-propelled particle that swims using prescribed distortions of its surface. In the most common case where radial distortions are ignored, the squirmer generates a tangential slip velocity on its surface, just like phoretic particles, which can be expanded into spherical harmonics modes (Blake Reference Blake1971; Pak & Lauga Reference Pak and Lauga2014). Consistently with the phoretic problem presented above, only the first two modes are included in the following.

The FCM force distribution produced by $N$ microswimmers self-propelling with a surface slip velocity is given by

(3.36)\begin{equation} \boldsymbol{f}(\boldsymbol{r},t) = \sum_{n=1}^N [{\boldsymbol{\mathsf{S}}}_n \boldsymbol{\cdot} \boldsymbol{\nabla} \varDelta^*(r_n) + {\boldsymbol{\mathsf{S}}}^a_n \boldsymbol{\cdot}\boldsymbol{\nabla} \varDelta(r_n) + \boldsymbol{H}^a_n \nabla^2 \varDelta^*(r_n)], \end{equation}

where ${\boldsymbol{\mathsf{S}}}^a_n$ is the active stresslet and $\boldsymbol {H}^a_n$ is the active potential dipole associated with the swimming disturbances of swimmer $n$. The latter is defined as

(3.37)\begin{equation} \boldsymbol{H}^a_n ={-}2 {\rm \pi}a^3 \boldsymbol{U}^a_n, \end{equation}

where $\boldsymbol {U}^a_n$ is the swimming velocity arising from the slip velocity on the swimmer surface $\boldsymbol {u}^s$ (2.7a,b). Note that the rigidity stresslet ${\boldsymbol{\mathsf{S}}}_n$ is included in (3.36) to enforce the absence of deformation of the swimmers, (3.35). The resulting velocity field reads

(3.38)\begin{equation} \boldsymbol{u}(\boldsymbol{r},t) = \sum_{n=1}^{N} [{\boldsymbol{\mathsf{S}}}_n : {\boldsymbol{\mathsf{R}}}^*(\boldsymbol{r}_n) + {\boldsymbol{\mathsf{S}}}^a_n : {\boldsymbol{\mathsf{R}}}(\boldsymbol{r}_n) + \boldsymbol{H}^a_n \cdot {\boldsymbol{\mathsf{A}}}^*(\boldsymbol{r}_n) ], \end{equation}

where ${\boldsymbol{\mathsf{R}}}$ is the FCM dipole Green's function evaluated with the parameter $\sigma$ instead of $\sigma _*$. The second-order tensor ${\boldsymbol{\mathsf{A}}}^*$ is the FCM Green's function for the potential dipole

(3.39)\begin{equation} {\boldsymbol{\mathsf{A}}}^*(\boldsymbol{r}) = \frac{1}{4{\rm \pi} r^3} \left[{\boldsymbol{\mathsf{I}}} - \frac{3\boldsymbol{r}\boldsymbol{r}}{r^2}\right] \mathrm{erf}\left(\frac{r}{\sigma_* \sqrt{2}}\right) - \frac{1}{\mu} \left[\left({\boldsymbol{\mathsf{I}}} - \frac{\boldsymbol{r}\boldsymbol{r}}{r^2} \right) + \left({\boldsymbol{\mathsf{I}}} - \frac{3\boldsymbol{r}\boldsymbol{r}}{r^2} \right) \left(\frac{\sigma_*}{r}\right)^2 \right] \varDelta^*(r). \end{equation}

The particles’ velocity, angular velocity and mean strain rate are then computed as

(3.40)$$\begin{gather} \boldsymbol{U}_n = \boldsymbol{U}^a_n -\boldsymbol{W}_n + \int_{V_F} \boldsymbol{u}\varDelta(\boldsymbol{r}_n) \,\mathrm{d}V, \end{gather}$$

(3.41)$$\begin{gather}\boldsymbol{\varOmega}_n = \boldsymbol{\varOmega}^a_n + \tfrac{1}{2}\int_{V_F} [\boldsymbol{\nabla} \times \boldsymbol{u}] \varDelta^*(\boldsymbol{r}_n)\,\mathrm{d}V, \end{gather}$$

(3.42)$$\begin{gather}{\boldsymbol{\mathsf{E}}}_n ={-}{\boldsymbol{\mathsf{K}}}_n + \tfrac{1}{2} \int_{V_F} [\boldsymbol{\nabla} \boldsymbol{u} + (\boldsymbol{\nabla} \boldsymbol{u})^\mathrm{T}] \varDelta^*(\boldsymbol{r}_n)\,\mathrm{d}V = \boldsymbol{0}, \end{gather}$$

where the active swimming velocities $\boldsymbol {U}^a_n$ and rotation rates $\boldsymbol {\varOmega }^a_n$ correspond to the intrinsic velocities of particle $n$, if it was alone (i.e. in the absence of external flows or other particles), and $\boldsymbol {W}_n$ and ${\boldsymbol{\mathsf{K}}}_n$ are defined as

(3.43)$$\begin{gather} \boldsymbol{W}_n = \int_{V_F} [\boldsymbol{H}^a_n\cdot{\boldsymbol{\mathsf{A}}}^*(\boldsymbol{r}_n)] \varDelta(\boldsymbol{r}_n)\,\mathrm{d}V, \end{gather}$$

(3.44)$$\begin{gather}{\boldsymbol{\mathsf{K}}}_n = \tfrac{1}{2} \int_{V_F} [{\boldsymbol{\mathsf{S}}}^a_n:\boldsymbol{\nabla} {\boldsymbol{\mathsf{R}}}(\boldsymbol{r}_n) + ({\boldsymbol{\mathsf{S}}}_n^a:\boldsymbol{\nabla} {\boldsymbol{\mathsf{R}}}(\boldsymbol{r}_n))^\mathrm{T}] \varDelta^*(\boldsymbol{r}_n)\,\mathrm{d}V, \end{gather}$$

and are included to subtract away the spurious self-induced velocities and local rates of strain arising from the integration of the full velocity field $\boldsymbol {u}$, which already includes the contribution of $\boldsymbol {H}^a_n$ and ${\boldsymbol{\mathsf{S}}}_n^a$ (Delmotte et al. Reference Delmotte, Keaveny, Plouraboué and Climent2015).

3.3.1. DFCM: coupling reactive and hydrodynamic FCM

The active swimming speed $\boldsymbol {U}^a_n$ involved in the potential dipole $\boldsymbol {H}^a_n$, (3.37), is the phoretic response of particle $n$ to the chemical field, if it was hydrodynamically isolated (i.e. neglecting the presence of other particles in solving the swimming problem). It thus includes its self-induced velocity (i.e. the response to the concentration contrasts induced by its own activity) and the drift velocity induced by the activity of the other particles. The swimming problem for a hydrodynamically isolated particle in unbounded flows can be solved directly using the reciprocal theorem (Stone & Samuel Reference Stone and Samuel1996), and using the definition of the phoretic slip flow

(3.45)\begin{equation} \boldsymbol{U}^a_n ={-} \langle \boldsymbol{u}^s \rangle_n ={-} \langle M \nabla_{{\parallel}} c \rangle_n. \end{equation}

After substitution of the mobility distribution at the surface of particle $n$, (2.9), using a truncated multipolar expansion of the surface concentration on particle $n$ (up to its second-order moment) and integration by parts, the intrinsic swimming velocity is obtained in terms of the first two surface concentration moments (see Appendix B for more details)

(3.46)\begin{equation} \boldsymbol{U}^a_n ={-}\frac{2\bar{M}_n}{a} \langle c\boldsymbol{n} \rangle_n -\frac{15 M^*_n}{8a} [2 \langle c(\boldsymbol{n}\boldsymbol{n}-{\boldsymbol{\mathsf{I}}}/3) \rangle_n \cdot \boldsymbol{p}_n + \left( \langle c(\boldsymbol{n}\boldsymbol{n}-{\boldsymbol{\mathsf{I}}}/3) \rangle_n : \boldsymbol{p}_n\boldsymbol{p}_n \right) \boldsymbol{p}_n]. \end{equation}

Similarly, the active stresslet ${\boldsymbol{\mathsf{S}}}^a_n$, is defined as in (2.8),

(3.47)\begin{equation} {\boldsymbol{\mathsf{S}}}^a_n ={-}10{\rm \pi} a^2 \langle \boldsymbol{n} \boldsymbol{u}^s +\boldsymbol{u}^s \boldsymbol{n} \rangle_n ={-}10{\rm \pi} a^2 \langle M ( \boldsymbol{n} \nabla_{{\parallel}} c + (\nabla_{{\parallel}} c) \boldsymbol{n} ) \rangle_n, \end{equation}

and rewrites in terms of the moments of concentration (see Appendix B for more details)

(3.48)\begin{align} {\boldsymbol{\mathsf{S}}}^a_n ={-}60 {\rm \pi}a \bar{M}_n \langle c(\boldsymbol{n}\boldsymbol{n}-{\boldsymbol{\mathsf{I}}}/3) \rangle_n + \frac{15 {\rm \pi}a M^*_n}{2} [(\langle c\boldsymbol{n} \rangle_n \cdot \boldsymbol{p}_n) ({\boldsymbol{\mathsf{I}}}-\boldsymbol{p}_n\boldsymbol{p}_n) - \langle c\boldsymbol{n} \rangle_n \boldsymbol{p}_n - \boldsymbol{p}_n \langle c\boldsymbol{n} \rangle_n]. \end{align}

Finally, the active rotation $\boldsymbol {\varOmega }^a_n$, (2.7a,b), is obtained in terms of the moments of concentration and the mobility contrast (see Appendix B)

(3.49)\begin{equation} \boldsymbol{\varOmega}^a_n = -\frac{9M^*_n}{4a^2} \boldsymbol{p}_n \times \langle c\boldsymbol{n} \rangle_n. \end{equation}

For uniform mobility, the swimming velocity and stresslet are directly related to the first and second of surface concentrations, but non-uniform mobility introduces a coupling of the different concentration moments. Here, the surface concentration is expanded up to its second-order moment only.

In our regularized approach, the surface concentration moments appearing in the previous equations will conveniently be computed as weighted volume averages over the entire domain $V_F$ as detailed in (3.14) and (3.16).

Computing the second moment of concentration, however, requires an additional step: as detailed in § 3.1.5, the second moment of concentration in an external field arises from the second gradient of that external field, and includes both an externally induced component $\langle c_E (\boldsymbol {n}\boldsymbol {n}-{\boldsymbol{\mathsf{I}}}/3) \rangle _n$ (i.e. the moment of that externally imposed field) and a self-induced component which corresponds to the second moment of the induced field generated by the particle to ensure that the correct flux boundary condition is satisfied at the particles’ surface. For a chemically inert particle ($\alpha =0$), the self-induced contribution is obtained exactly as ${\langle c_I^o (\boldsymbol {n}\boldsymbol {n}-{\boldsymbol{\mathsf{I}}}/3) \rangle _n = \frac {2}{3}\langle c_E (\boldsymbol {n}\boldsymbol {n}-{\boldsymbol{\mathsf{I}}}/3) \rangle }$.

Our representation of the particles in the chemical problem is, however, truncated at the dipole level, (3.9), and as a result, the quadrupolar response of the particle to the external field cannot be accounted for directly. To correct for this shortcoming, we first compute the external second moment produced by the other particles on particle $n$ using (3.16) and (3.9), and multiply the resulting value by $5/3$ to account for the full second moment induced by the concentration field indirectly.

Finally, the particles are themselves active and may generate an intrinsic quadrupole. Its effect on the second surface concentration moment can be added explicitly in terms of the second activity moment, so that the total second moment on particle $n$ is finally evaluated as

(3.50)\begin{equation} \langle c (\boldsymbol{n}\boldsymbol{n}-{\boldsymbol{\mathsf{I}}}/3) \rangle_n = \frac{5}{3}\{ c_E (\boldsymbol{n}\boldsymbol{n}-{\boldsymbol{\mathsf{I}}}/3)\}_n + \frac{a}{D}\langle \alpha (\boldsymbol{n}\boldsymbol{n}-{\boldsymbol{\mathsf{I}}}/3) \rangle_n. \end{equation}

In summary, at a given time step, the particles’ velocities are obtained from their instantaneous position and orientation as follows. The first two surface concentration moments are first obtained using our new reactive FCM framework by solving the Poisson problem (3.7). These moments are then used to compute the phoretic intrinsic translation and rotation velocities (3.46) and (3.49), as well as the active stresslets and potential dipoles (3.48) and (3.37). The Stokes equations forced by the swimming singularities (3.36), and subject to the particle rigidity constraint (3.42), are finally solved to obtain the total particle velocities (3.40) and (3.41).

3.3.2. Numerical details

The volume integrals required to compute the concentration moments and the hydrodynamic quantities are performed with a Riemann sum on Cartesian grids centred at each particle position. To ensure a sufficient resolution, the grid size, $\varDelta x$, is chosen so that the smallest envelope size $\sigma _D$ satisfies ${\sigma _{D} = 1.5\varDelta x = 0.3614a}$, which corresponds to approximately 4 grid points per radius. Owing to the fast decay of the envelopes, the integration domain is truncated so that the widest envelope (that with the largest $\sigma$) essentially vanishes on the boundary of the domain, $\varDelta (r) < \gamma = 10^{-16}$, which, given the grid resolution, requires 39 integration points in each direction. Doing so, the numerical integrals yield spectral accuracy. Setting instead $\gamma = \epsilon = 10^{-10}$, where $\epsilon$ is the relative tolerance for the polarity in the iterative procedure (3.13), reduces that number to 31 integration points along each axis while maintaining spectral convergence.