Fluctuating hydrodynamics of an autophoretic particle near a permeable interface

Abstract

We study the autophoretic motion of a spherical active particle interacting chemically and hydrodynamically with its fluctuating environment in the limit of rapid diffusion and slow viscous flow. Then, the chemical and hydrodynamic fields can be expressed in terms of integrals. The resulting boundary-domain integral equations provide a direct way of obtaining the traction on the particle, requiring the solution of linear integral equations. An exact solution for the chemical and hydrodynamic problems is obtained for a particle in an unbounded domain. For motion near boundaries, we provide corrections to the unbounded solutions in terms of chemical and hydrodynamic Green's functions, preserving the dissipative nature of autophoresis in a viscous fluid for all physical configurations. Using this, we give the fully stochastic update equations for the Brownian trajectory of an autophoretic particle in a complex environment. First, we analyse the Brownian dynamics of particles capable of complex motion in the bulk. We then introduce a chemically permeable planar surface of two immiscible liquids in the vicinity of the particle and provide explicit solutions to the chemo-hydrodynamics of this system. Finally, we study the case of an isotropically phoretic particle hovering above an interface as a function of interfacial solute permeability and viscosity contrast.

1. Introduction

Autophoretic motion comprises the propulsion of particles due to self-generated gradients (Anderson 1989; Paxton et al. 2006; Ebbens & Howse 2010; Moran & Posner 2017), typically on an energy scale comparable to that of thermal fluctuations (Batchelor 1976; Graham 2018). This self-propulsion mechanism allows systems of phoretic particles to mimic the locomotion of microorganisms (Brennen & Winet 1977; Goldstein 2015), making them useful in the study of the fundamental principles of motility and collective behaviour (Palacci et al. 2013; Illien, Golestanian & Sen 2017; Shaebani et al. 2020; Zöttl & Stark 2023; Kumar et al. 2024). In particular, the study of interactions of autophoretic particles with nearby boundaries is relevant in micro-fluidics, biophysics and surface science (Kreuter et al. 2013; Ibrahim & Liverpool 2015; Uspal et al. 2015; Shen, Würger & Lintuvuori 2018; Thutupalli et al. 2018; Singh, Adhikari & Cates 2019).

Our goal is thus to formulate an effective description in which Brownian motion and autophoresis of active particles can be studied when suspended in a complex environment. Models for self-diffusiophoresis typically assume that chemical gradients generated by the particle induce an osmotic pressure, which is balanced by viscous stresses driving an effective slip flow confined to a thin layer at the surface of the particle (Anderson, Lowell & Prieve 1982). This sets the surrounding fluid in motion, with fluid stresses reacting back on the particle and setting it in motion. To compute the particle dynamics usually requires solving for the concentration field and the fluid flow in the bulk, and subsequently obtaining the stresses on the particle by matching all relevant boundary conditions (Golestanian, Liverpool & Ajdari 2005, 2007). Instead, by using a boundary-domain integral approach, we directly obtain the concentration distribution and the resulting traction (force per unit area) on the surface of the particle, obviating the need for solving the governing equations in the bulk. Compared with more conventional kinematic approaches (Lighthill 1952; Pak & Lauga 2014), it is then straightforward to incorporate thermal fluctuations in the surrounding fluid as Brownian stresses on the particle. The latter have been studied extensively for suspensions of colloidal particles (Einstein 1905; Zwanzig 1964; Chow 1973; Hinch 1975; Ermak & McCammon 1978; Ladd 1994; Cichocki et al. 2000; Keaveny 2014; Delmotte & Keaveny 2015; Singh & Adhikari 2017; Bao et al. 2018; Mozaffari et al. 2018; Elfring & Brady 2022; Turk, Singh & Adhikari 2022; Westwood, Delmotte & Keaveny 2022), highlighting that any acceptable approximation of the colloidal diffusion matrix in Brownian dynamics modelling must remain positive–definite for all physical configurations (Wajnryb, Szymczak & Cichocki 2004; Wajnryb et al. 2013). Based on a Galerkin–Jacobi iterative method, the analytical expressions we provide naturally satisfy this condition.

The fields generated by and the resulting stresses on autophoretic particles are well known in an unbounded fluid (Golestanian et al. 2007; Ebbens & Howse 2010; Illien et al. 2017; Lisicki, Reigh & Lauga 2018), or when confined either by no-slip walls that are impermeable to the solutes (Crowdy 2013; Ibrahim & Liverpool 2015; Uspal et al. 2015; Ibrahim & Liverpool 2016; Mozaffari et al. 2016; Daddi-Moussa-Ider et al. 2018; Kanso & Michelin 2019; Singh et al. 2019), or by chemically patterned boundaries (Uspal et al. 2019). In this paper we formulate a general framework for finding the full chemo-hydrodynamics of a particle in an arbitrary complex environment in terms of chemical and hydrodynamic Green's functions (the fields generated by Dirac delta function sources; Ladyzhenskaia 1969). Using this, we provide analytically the dynamics of a phoretic particle in the proximity of a chemically permeable liquid–liquid interface separating the suspending domain from a second, immiscible liquid phase. Assuming a large capillary number, we restrict our considerations to a planar interface. This is particularly relevant for studies on particle aggregation near fluid–fluid interfaces and free surfaces (Chen et al. 2015; Hokmabad et al. 2022), with a permeable interface being a plausible model of biofilms and hydrogels (Wichterle & Lím 1960; Berke et al. 2008).

The rest of the paper is organised as follows. In § 2, we review the chemo-hydrodynamic problem of autophoresis in a fluctuating environment and its formal solution via the boundary-domain integral representation of Laplace and Stokes equations. In § 3, we then use a Galerkin discretisation to project the formal solution onto a basis of tensor spherical harmonics (TSH), finding an exact and an approximate solution to the full chemo-hydrodynamic problem far away from and near boundaries, respectively. We provide the stochastic update equations for thermally agitated autophoresis in complex environments. In § 4, we apply these equations to the study of three representative examples. First, we consider systematically patterned particle surfaces, which we confirm can lead to complex phoretic motion even in the absence of boundaries (Lisicki et al. 2018). We study the effect of thermal fluctuations on the resulting particle motion in a bulk fluid. In the vicinity of the particle we then introduce the presence of a plane surface of two immiscible liquids that is permeable to the solutes. For this system, we obtain explicit forms of the relevant chemical and hydrodynamic connectors. We demonstrate our analytical results by numerically investigating the chemo-hydrodynamic effects the interface has on the dynamics of a nearby autophoretic particle. This includes an analysis of the hovering state of a phoretic particle above an interface as a function of particle activity, and interfacial properties. We conclude with a brief discussion of our results and potential future applications thereof in § 5.

2. Chemo-hydrodynamics

We consider a spherical autophoretic particle of radius $b$, suspended in an incompressible fluid ($\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {v}=0$, where $\boldsymbol {v}$ is the flow field) of viscosity $\eta$ at low Reynolds number. Thermal fluctuations of the fluid at equilibrium are modelled by a zero-mean Gaussian random field $\boldsymbol {\xi }$, the thermal force acting on the particle, whose variance is given by a fluctuation–dissipation relation (Zwanzig 1964; Fox & Uhlenbeck 1970; Hauge & Martin-Löf 1973; Bedeaux & Mazur 1974; Roux 1992). In table 1, we summarise the differential laws governing the chemo-hydrodynamics of this system. We denote fields defined on the surface of spherical particles as functions of the radius vector $\boldsymbol {b}$ of the sphere, where $\hat {\boldsymbol {b}}=\boldsymbol {b}/b$ is the unit outward normal to the surface, pointing into the fluid and with $b=|\boldsymbol{b}|$. We assume a negligibly small Péclet number, thus ignoring distortions induced by the flow on the solute concentration (Michelin, Lauga & Bartolo 2013; Morozov & Michelin 2019). Additionally, we assume that solute diffusion takes place on much shorter time scales than Brownian motion of the autophoretic particle, which in turn takes place on much shorter time scales than its rigid-body motion. The chemical problem is then represented by the Laplace equation for the concentration field $c$, for ideal solutions equivalent to a divergence-free chemical flux $\boldsymbol {j}$ in (Ia), where $D$ is the solute diffusivity in the fluid. In (Ic) the normal component of the flux at the surface of the particle $j^{\mathcal {A}}(\boldsymbol {b})$ is specified.

Table 1. Governing differential laws. This table summarises the chemo-hydrodynamic coupling at the surface of an autophoretic particle, an example of two three-dimensional partial differential equations, namely the Laplace equation for the concentration field (Ia) and the Stokes equation for the fluid flow and pressure in a fluid with thermal force density $\boldsymbol {\xi }$ (Id), coupling on a two-dimensional surface only (Melcher & Taylor 1969). The chemo-hydrodynamic coupling (Ig) leads to the specified active flux $j^{\mathcal {A}}$ driving a slip flow $\boldsymbol {v}^{\mathcal {A}}$ in a thin layer at the surface of the particle with specified phoretic mobility $\mu _{c}$, finally driving the fluid surrounding the particle and causing self-propulsion. A passive particle is a rigid sphere of radius $b$ with the boundary condition: $\boldsymbol {v}^{\mathcal {D}}(\boldsymbol {b})=\boldsymbol {V}+\boldsymbol {\varOmega }\times \boldsymbol {b}$, where $\boldsymbol {V}$ is the velocity and $\boldsymbol {\varOmega }$ is the angular velocity of the particle. An active particle is modelled as a sphere with boundary condition (If), which comprises both slip $\boldsymbol {v}^{\mathcal {A}}(\boldsymbol {b})$ and rigid-body motion $\boldsymbol {v}^{\mathcal {D}}(\boldsymbol {b})$ (Anderson 1989; Ebbens & Howse 2010).

Surface gradients of the generated concentration field induce a mass transport of solute, thus driving a fluid flow confined to a thin layer at the surface of the particle. This is modelled by a slip $\boldsymbol {v}^{\mathcal {A}}$ in the chemo-hydrodynamic coupling in (Ig). Here, $\mu _{c}$ is the particle-specific phoretic mobility, which incorporates the solute–colloid interactions. We assume that the solute is uncharged (neutral diffusiophoresis) (Prieve et al. 1984; Velegol et al. 2016; Yang, Rallabandi & Stone 2019). The slip is incorporated in the velocity boundary condition in (If), alongside rigid-body motion $\boldsymbol {v}^{\mathcal {D}}$ of the particle. Finally, the particle sets the surrounding fluid in motion (via the slip or rigid-body motion due to external forces and torques), hydrodynamically interacting with its surroundings via the Stokes equation (Id). Therein, we have defined the Cauchy stress tensor $\boldsymbol {\sigma }$, containing contributions from the isotropic fluid pressure $p$ and from spatial variations in the flow field. Here $\boldsymbol {I}$ is the identity tensor.

In table 2, we summarise the boundary-domain integral equations (BIEs) corresponding to the Laplace and Stokes equations, and their formal solution in terms of integral linear operators. The BIE (IIa) for the concentration at the surface of the particle is given in terms of a background concentration field $c^{\infty }(\boldsymbol {r})$, the single-layer operator $\mathcal {H}[j^{\mathcal {A}}]$ and the double-layer operator $\mathcal {L}[c]$. This naming convention of the integral operators is by analogy with potential theory (Jackson 1962; Kim & Karrila 2005). The integral kernels contain the concentration Green's function $H$ and its gradient $\boldsymbol {L}$. Due to linearity of the Laplace equation, we can find the solution in (IId) for the concentration, containing the operator $\zeta$ for the linear response to a background chemical field and the so-called elastance operator $\mathcal {E}$. The naming convention of the latter originates from Maxwell, who in his study of the capacitance of a system of spherical conductors coined the term elastance for the isotropic part of the tensor $\boldsymbol {\mathcal {E}}$ (Maxwell 1881).

Table 2. Governing integral laws. This table summarises the formal solutions to the boundary-domain integral equations corresponding to the Laplace equation for the concentration $c$ (IIa) and the Stokes equation for the traction (force per unit area) $\boldsymbol {f}$ (IIf) on the surface of the particle. Here, $\boldsymbol {r},\tilde {\boldsymbol {r}}\in S$, where $\boldsymbol {r}=\boldsymbol {R}+\boldsymbol {b}$ and $\tilde {\boldsymbol {r}}=\boldsymbol {R}+\boldsymbol {b}'$ are the field and source points at the surface $S$ of the particle centred at $\boldsymbol {R}$, respectively, and $\int \mathrm {d}S'$ implies an integration over $\tilde {\boldsymbol {r}}$. In (IId) and (IIj) we give the solutions for the concentration and the traction in terms of integral linear operators. We have used the fact that rigid-body motion $\boldsymbol {v}^{\mathcal {D}}$ lies in the eigenspectrum of the double-layer operator (IIh) with an eigenvalue $-1/2$ (Kim 2015). In the formal solutions we have introduced operators representing the linear response to a background concentration field $\zeta$, the so-called elastance $\mathcal {E}$, the rigid-body friction $\boldsymbol {\gamma }$ and the friction due to surface slip $\hat {\boldsymbol {\gamma }}$. Inserting the operator solution for the concentration in (IIe) into (Ig) in table 1 for the chemo-hydrodynamic coupling at the surface of an autophoretic particle, we find (IIl) for the surface slip $\boldsymbol {v}^{\mathcal {A}}(\boldsymbol {b})$.

The corresponding BIE of fluctuating Stokes flow (IIf) is a sum of the single-layer operator $\boldsymbol {\mathcal {G}}[\boldsymbol {f}]$ acting on the surface traction (force per unit area) on the particle, given by $\boldsymbol {f}=\boldsymbol {\sigma }\boldsymbol {\cdot }\hat {\boldsymbol {b}}$, the double-layer operator $\boldsymbol {\mathcal {K}}[\boldsymbol {v}]$ (Lorentz 1896; Odqvist 1930; Ladyzhenskaia 1969; Youngren & Acrivos 1975; Zick & Homsy 1982; Pozrikidis 1992; Muldowney & Higdon 1995; Cheng & Cheng 2005; Leal 2007; Singh, Ghose & Adhikari 2015) and the Brownian velocity field $\boldsymbol {u}[\boldsymbol {\xi }]$ (Singh & Adhikari 2017). The integral kernels contain the Green's function $\boldsymbol {G}$ of the Stokes equation and the stress tensor $\boldsymbol {K}$ associated with it. Linearity of the Stokes equation allows us to formally solve the BIE, introducing the friction operators $\boldsymbol {\gamma }$ and $\hat {\boldsymbol {\gamma }}$ due rigid-body motion and slip, respectively. They can be distinguished by a non-trivial contribution of the double-layer integral to the latter (Turk et al. 2022). Finally, the solutions to the chemical and hydrodynamic problems are coupled via the boundary condition (IIl).

In the following, an autophoretic particle is fully specified by its surface flux $j^{\mathcal {A}}$ and phoretic mobility $\mu _{c}$, as indicated in table 1. Our aim is to find its dynamics, governed by Newton's laws

(1.1a,b)\begin{equation} m\dot{\boldsymbol{V}}=\boldsymbol{F}^{H}+\boldsymbol{F}^{P}+\hat{\boldsymbol{F}},\quad I\dot{\boldsymbol{\varOmega}}=\boldsymbol{T}^{H}+\boldsymbol{T}^{P}+\hat{\boldsymbol{T}}.\end{equation}

Here, $m$ and $I$ are the particle mass and moment of inertia, respectively, and a dotted variable implies a time derivative. Body forces and torques are denoted by $\boldsymbol {F}^{P}$ and $\boldsymbol {T}^{P}$, and the hydrodynamic and fluctuating contributions are defined in terms of the traction on the particle

(1.2a–d)\begin{equation} \boldsymbol{F}^{H}=\int\boldsymbol{f}^{H}\,{\rm d}S,\quad\boldsymbol{T}^{H}= \int\boldsymbol{b}\times\boldsymbol{f}^{H}\,{\rm d}S,\quad \hat{\boldsymbol{F}}=\int\hat{\boldsymbol{f}}\,{\rm d}S,\quad \hat{\boldsymbol{T}}=\int\boldsymbol{b}\times\hat{\boldsymbol{f}}\,{\rm d}S,\end{equation}

where the total surface traction on the particle is the sum $\boldsymbol {f}=\boldsymbol {f}^{H}+\hat {\boldsymbol {f}}$. We define the hydrodynamic traction due to rigid-body and active interactions as $\boldsymbol {f}^{H}$ and the Brownian traction due to thermal fluctuations in the fluid as $\hat {\boldsymbol {f}}$ such that

(1.3a,b)\begin{equation} \boldsymbol{f}^{H} ={-}\boldsymbol{\gamma}\left[\boldsymbol{v}^{\mathcal{D}}\right]-\hat{\boldsymbol{\gamma}} \left[\boldsymbol{v}^{\mathcal{A}}\right],\quad\hat{\boldsymbol{f}}=\boldsymbol{\gamma}\left[\boldsymbol{u}\left[\boldsymbol{\xi}\right]\right]. \end{equation}

It is known that the latter are zero-mean random variables with variances fixed by a fluctuation–dissipation relation (Zwanzig 1964; Chow 1973). By linearity of the governing equations, the hydrodynamic and Brownian contributions can be solved for independently and the fluid degrees of freedom can be eliminated exactly, yielding the Brownian dynamics of the active particle.

3. Solution in an irreducible basis

In this section, we write the formal solutions to Laplace and Stokes equations in table 2, (IIe) and (IIk), in an irreducible basis, thus transforming the integral operator equations into linear systems, for which we give explicit solutions. We choose a basis of TSH, defined by

(3.1)\begin{equation} Y_{\alpha_{1}\dots\alpha_{l}}^{(l)}(\hat{\boldsymbol{b}})=(2l-1)!!\varDelta_{\alpha_{1}\dots\alpha_{l}, \beta_{1}\dots\beta_{l}}^{(l)}\hat{b}_{\beta_{1}}\dots\hat{b}_{\beta_{l}}=({-}1)^{l}\,b^{l+1}\, \boldsymbol{\nabla}_{\alpha_{1}}\dots\boldsymbol{\nabla}_{\alpha_{l}}\frac{1}{b}, \end{equation}

where $\boldsymbol {\varDelta }{}^{(l)}$ is a rank-$2l$ tensor, which projects a tensor of rank-$l$ onto its symmetric and traceless part (Hess 2015).

3.1. Chemical problem

To project (IIe) for the concentration at the surface of the particle onto a linear system, we expand the boundary fields

(3.2a,b)\begin{equation} c(\boldsymbol{b})=\sum_{q=0}^{\infty}w_{q}\boldsymbol{C}^{(q)}\odot\boldsymbol{Y}^{(q)}(\hat{\boldsymbol{b}}),\quad j^{\mathcal{A}}(\boldsymbol{b})=\sum_{q=0}^{\infty}\tilde{w}_{q}\boldsymbol{J}^{(q)} \odot\boldsymbol{Y}^{(q)}(\hat{\boldsymbol{b}}).\end{equation}

The product denoted by $\odot$ implies a maximal contraction of Cartesian indices (a $\boldsymbol {q}$-fold contraction between a tensor of rank-$q$ and another one of higher rank) and we have defined

(3.3a,b)\begin{equation} w_{q}=\frac{1}{q!\left(2q-1\right)!!},\quad\tilde{w}_{q}=\frac{2q+1}{4{\rm \pi} b^{2}}. \end{equation}

The expansion coefficients $\boldsymbol {C}^{(q)}$ and $\boldsymbol {J}^{(q)}$ are symmetric and traceless tensors of rank-$q$. The background concentration field $c^{\infty }(\boldsymbol {b})$ at the surface of the particle is expanded in an analogous manner to $c(\boldsymbol {b})$, with coefficients denoted by $\boldsymbol {C}^{\infty (q)}$. Linearity of the Laplace equation implies that the general solution in a basis of TSH can be written as

(3.4)\begin{equation} \boldsymbol{C}^{(q)}=\boldsymbol{\zeta}^{(q,q')}\odot\boldsymbol{C}^{\infty(q')}+\boldsymbol{\mathcal{E}}^{(q,q')}\odot\boldsymbol{J}^{(q')},\end{equation}

corresponding to (IIe), where the task now is to find the connecting tensors $\boldsymbol {\zeta }^{(q,q')}$ and $\boldsymbol {\mathcal {E}}^{(q,q')}$. In Appendix A, starting from the BIE for the surface concentration and using a Galerkin–Jacobi iterative method, we outline how to find approximate solutions, in leading powers of distance between the particle and surrounding boundaries, for these tensors in terms of a given Green's function $H$ of Laplace equation (Singh et al. 2019).

Any Green's function $H$ of Laplace equation can be written as the sum

(3.5)\begin{equation} H(\boldsymbol{R},\tilde{\boldsymbol{R}})=H^{o}(\boldsymbol{r})+H^{*}(\boldsymbol{R},\tilde{\boldsymbol{R}}),\end{equation}

with $\boldsymbol {r}=\boldsymbol {R}-\tilde {\boldsymbol {R}}$, where $\boldsymbol {R}$ and $\tilde {\boldsymbol {R}}$ are the field and the source point, respectively. Here, $H^{o}(\boldsymbol {r})=1/4{\rm \pi} Dr$ is the fundamental solution of Laplace equation in an unbounded domain. On the other hand, $H^{*}$ is an extra contribution needed to satisfy additional boundary conditions in the system. For the unbounded case, where $H=H^{o}(\boldsymbol {r})$, the single-layer and double-layer operators in (IIb) and (IIc) have singular integral kernels. However, due to translational invariance they can be evaluated using Fourier techniques, see Appendix A.1. We find that both integral operators diagonalise simultaneously in a basis of TSH, yielding

(3.6a,b)\begin{equation} \boldsymbol{\zeta}^{(q,q')}=\zeta_{q}\,\boldsymbol{I}^{(q,q')},\quad \boldsymbol{\mathcal{E}}^{(q,q')}=\mathcal{E}_{q}\,\boldsymbol{I}^{(q,q')},\end{equation}

where

(3.7a,b)\begin{equation} \zeta_{q}=\frac{2q+1}{q+1},\quad\mathcal{E}_{q}=\frac{b}{D}\frac{\tilde{w}_{q}}{w_{q}}\frac{1}{q+1},\end{equation}

and $\boldsymbol {I}^{(q,q')}$ is a tensor with elements $\delta _{qq'}$, where the latter denotes a Kronecker delta. The expression for the elastance $\mathcal {E}_{q}$ is confirmed by previous results obtained by first solving the Laplace equation in the fluid volume and subsequently matching the boundary condition (Ic) for the surface flux (Jackson 1962; Golestanian et al. 2007). If the system contains additional boundaries, we find corrections to these diagonal expressions in terms of derivatives of $H^{*}$. To leading order, this yields

(3.8)\begin{equation} \boldsymbol{\zeta}^{(q,q')}=\zeta_{q}\left(\boldsymbol{I}^{(q,q')}+4{\rm \pi} bDw_{q'}\frac{q'}{q'+1}b^{q+q'}\boldsymbol{\nabla}^{(q)} \tilde{\boldsymbol{\nabla}}^{(q')}H^{*}(\boldsymbol{R},\tilde{\boldsymbol{R}})\right),\end{equation}

where $\boldsymbol {\nabla }_{\alpha _{1}\dots \alpha _{q}}^{(q)}=\boldsymbol {\nabla }_{\alpha _{1}}\dots \boldsymbol {\nabla }_{\alpha _{q}}$, and where we have introduced the short-hand notation $\boldsymbol {\nabla }_{\boldsymbol {R}}=\boldsymbol {\nabla }$ for derivatives with respect to the field point and $\boldsymbol {\nabla }_{\tilde {\boldsymbol {R}}}=\tilde {\boldsymbol {\nabla }}$ for the source point. Similarly, we find for the elastance

(3.9)\begin{equation} \boldsymbol{\mathcal{E}}^{(q,q')}=\mathcal{E}_{q}\left(\boldsymbol{I}^{(q,q')}+4{\rm \pi} bDw_{q}\frac{2q'+1}{q'+1}b^{q+q'}\boldsymbol{\nabla}^{(q)}\tilde{\boldsymbol{\nabla}}^{(q')}H^{*} (\boldsymbol{R},\tilde{\boldsymbol{R}})\right).\end{equation}

In these expressions, the point of evaluation, $\boldsymbol {R}=\tilde {\boldsymbol {R}}$, for the one-body problem, is left implicit for brevity.

3.2. Hydrodynamic problem and Brownian motion

Using the linearity of Stokes flow we solve for the hydrodynamic traction $\boldsymbol {f}^{H}$ in a basis of TSH. Upon eliminating the hydrodynamic problem, Newton's equations (1.1a,b) will reveal the Brownian motion of an active particle. First, to find the linear system corresponding to (IIk), we expand the slip and the hydrodynamic traction in a basis of TSH

(3.10a,b)\begin{equation} \boldsymbol{v}^{\mathcal{A}}(\boldsymbol{b})=\sum_{l=1}^{\infty}w_{l-1}\boldsymbol{V}^{(l)}\odot\boldsymbol{Y}^{(l-1)}(\hat{\boldsymbol{b}}),\quad \boldsymbol{f}^{H}(\boldsymbol{b})=\sum_{l=1}^{\infty}\tilde{w}_{l-1}\boldsymbol{F}^{(l)}\odot\boldsymbol{Y}^{(l-1)}(\hat{\boldsymbol{b}}).\end{equation}

The coefficients $\boldsymbol {V}^{(l)}$ and $\boldsymbol {F}^{(l)}$ are rank-$l$ tensors, symmetric and traceless in their last $l-1$ indices. They can be decomposed into irreducible representations, denoted by $\boldsymbol {V}^{(l\sigma )}$ (or $\boldsymbol {F}^{(l\sigma )}$ for the traction moments), where $\boldsymbol {V}^{(ls)}$ (symmetric and traceless), $\boldsymbol {V}^{(la)}$ (anti-symmetric) and $\boldsymbol {V}^{(lt)}$ (trace) are irreducible tensors of ranks $l$, $l-1$ and $l-2$, respectively (Singh et al. 2015). For slip restricted by mass conservation only, obeying $\int \boldsymbol {v}^{\mathcal {A}}\boldsymbol {\cdot }\hat {\boldsymbol {b}}\,{\rm d}S=0$, these irreducible components of $\boldsymbol {V}^{(l)}$ (and $\boldsymbol {F}^{(l)}$) are independent of each other. In terms of the common definitions for the velocity and angular velocity of an active particle in an unbounded domain (Anderson & Prieve 1991; Stone & Samuel 1996; Ghose & Adhikari 2014)

(3.11a,b)\begin{equation} \boldsymbol{V}^{\mathcal{A}}={-}\frac{1}{4{\rm \pi} b^{2}}\int\boldsymbol{v}^{\mathcal{A}}(\boldsymbol{b})\,{\rm d}S,\quad \boldsymbol{\varOmega}^{\mathcal{A}}={-}\frac{3}{8{\rm \pi} b^{3}}\int\hat{\boldsymbol{b}}\times\boldsymbol{v}^{\mathcal{A}} (\boldsymbol{b})\,{\rm d}S,\end{equation}

we have $\boldsymbol {V}^{(1s)}=-\boldsymbol {V}^{\mathcal {A}}$ and $\boldsymbol {V}^{(2a)}/2b=-\boldsymbol {\varOmega }^{\mathcal {A}}.$ Similarly, we have, for the hydrodynamic force and torque defined in (1.1a,b), $\boldsymbol {F}^{(1s)}=\boldsymbol {F}^{H}$ and $\boldsymbol {F}^{(2a)}=({1}/{b})\boldsymbol {T}^{H}$.

Linearity of the Stokes equation then allows us to write down the deterministic part of (IIk) in a basis of TSH

(3.12)\begin{equation} \boldsymbol{F}^{(l\sigma)} ={-}\boldsymbol{\gamma}^{(l\sigma,1s)}\boldsymbol{\cdot}\boldsymbol{V}-\boldsymbol{\gamma}^{(l\sigma,2a)} \boldsymbol{\cdot}\boldsymbol{\varOmega}-\hat{\boldsymbol{\gamma}}^{(l\sigma,l'\sigma')} \odot\boldsymbol{V}^{(l'\sigma')},\end{equation}

where $\boldsymbol {\gamma }^{(l\sigma,l'\sigma ')}$ and $\hat {\boldsymbol {\gamma }}^{(l\sigma,l'\sigma ')}$ are generalised friction tensors for rigid-body motion and slip, respectively. For the modes corresponding to rigid-body motion, it is known that $\boldsymbol {\gamma }^{(l\sigma,1s)}=\hat {\boldsymbol {\gamma }}^{(l\sigma,1s)}$ and $\boldsymbol {\gamma }^{(l\sigma,2a)}=\hat {\boldsymbol {\gamma }}^{(l\sigma,2a)}$ (Singh & Adhikari 2018; Turk et al. 2022). Therefore, we can write for the hydrodynamic force and torque

(3.13) \begin{multline} \begin{pmatrix}\boldsymbol{F}^{H}\\ \boldsymbol{T}^{H} \end{pmatrix}=-\boldsymbol{\varGamma}\boldsymbol{\cdot}\begin{pmatrix}\boldsymbol{V}-\boldsymbol{V}^{\mathcal{A}}\\ \boldsymbol{\varOmega}-\boldsymbol{\varOmega}^{\mathcal{A}} \end{pmatrix}-\sum_{l\sigma=2s}\hat{\boldsymbol{\varGamma}}^{(l\sigma)}\odot\boldsymbol{V}^{(l\sigma)},\\ \text{with}\quad\boldsymbol{\varGamma}=\begin{pmatrix}\boldsymbol{\gamma}^{TT} & \boldsymbol{\gamma}^{TR}\\ \boldsymbol{\gamma}^{RT} & \boldsymbol{\gamma}^{RR} \end{pmatrix},\quad\hat{\boldsymbol{\varGamma}}^{(l\sigma)}=\begin{pmatrix}\hat{\boldsymbol{\gamma}}^{(T,l\sigma)}\\ \hat{\boldsymbol{\gamma}}^{(R,l\sigma)} \end{pmatrix}\end{multline}

where the superscripts $T$ and $R$ imply $l\sigma =1s,2a$, respectively, to confirm with existing literature (Ladd 1988). The matrix $\boldsymbol {\varGamma }$ contains the friction on the particle due to rigid-body motion, and $\hat {\boldsymbol {\varGamma }}^{(l\sigma )}$ contains the friction due to higher modes of slip. This concludes the solution of the hydrodynamic problem without fluctuations.

In a thermally fluctuating fluid, the Brownian forces and torques obey the fluctuation–dissipation relations (Einstein 1905; Zwanzig 1964; Chow 1973; Singh & Adhikari 2017)

(3.14a,b)\begin{equation} \left\langle \begin{pmatrix}\hat{\boldsymbol{F}}(t)\\ \hat{\boldsymbol{T}}(t) \end{pmatrix}\right\rangle =\boldsymbol{0},\quad\left\langle \begin{pmatrix}\hat{\boldsymbol{F}}(t)\\ \hat{\boldsymbol{T}}(t) \end{pmatrix}\begin{pmatrix}\hat{\boldsymbol{F}}(t')\\ \hat{\boldsymbol{T}}(t') \end{pmatrix}^{\text{tr}}\right\rangle =2k_{B}T\,\boldsymbol{\varGamma}\,\delta(t-t'),\end{equation}

where angled brackets denote ensemble averages, $k_{B}$ is the Boltzmann constant and $T$ is the temperature, while the transpose is defined as $(A_{\alpha \beta })^{\textrm {tr}}=A_{\beta \alpha }$. Inserting the above equations for the deterministic and stochastic forces and torques into Newton's equations (1.1a,b) yields the Langevin equation

(3.15)\begin{align} \begin{pmatrix}m\dot{\boldsymbol{V}}\\ I\dot{\boldsymbol{\varOmega}} \end{pmatrix}=\begin{pmatrix}\boldsymbol{F}^{P}\\ \boldsymbol{T}^{P} \end{pmatrix}-\boldsymbol{\varGamma}\boldsymbol{\cdot}\begin{pmatrix}\boldsymbol{V}-\boldsymbol{V}^{\mathcal{A}}\\ \boldsymbol{\varOmega}-\boldsymbol{\varOmega}^{\mathcal{A}} \end{pmatrix}-\sum_{l\sigma=2s}\hat{\boldsymbol{\varGamma}}^{(l\sigma)}\odot\boldsymbol{V}^{(l\sigma)}+\sqrt{2k_{B}T\,\boldsymbol{\varGamma}} \boldsymbol{\cdot}\begin{pmatrix}\boldsymbol{\xi}^{T}\\ \boldsymbol{\xi}^{R} \end{pmatrix}.\end{align}

The parameters $\boldsymbol {\xi }^{\alpha }$ are random variables with zero mean and unit variance. In the inertial equation (3.15) the noise is not multiplicative since $\boldsymbol {\varGamma }$ is configuration dependent, but not velocity dependent. With the particle centre of mass $\boldsymbol {R}$ and its unit orientation vector $\boldsymbol {e}$ (its orientation is governed by the rotational dynamics $\dot {\boldsymbol {\varTheta }}=\boldsymbol {\varOmega }$, where $\boldsymbol {\varTheta }$ is an arbitrary set of angles), we can find its Brownian trajectory by integrating

(3.16a,b)\begin{equation} \dot{\boldsymbol{R}}=\boldsymbol{V},\quad\dot{\boldsymbol{e}}=\boldsymbol{\varOmega}\times\boldsymbol{e},\end{equation}

over time. In colloidal systems the inertia of both the particles and the fluid are typically negligible. This corresponds to the Smoluchowski limit of (3.15). Adiabatic elimination of the momentum variables in phase space then directly leads to the following update equations in Itô form (Ermak & McCammon 1978; Gardiner 1984; Wajnryb et al. 2004; Volpe & Wehr 2016):

(3.17a) \begin{align} &\boldsymbol{R}(t+\Delta t) =\boldsymbol{R}(t)+\Delta\hat{\boldsymbol{R}}\nonumber\\ &\quad+\Big\{\boldsymbol{V}^{\mathcal{A}}+\boldsymbol{\mu}^{TT}\boldsymbol{\cdot}\boldsymbol{F}^{P}+\boldsymbol{\mu}^{TR}\boldsymbol{\cdot}\boldsymbol{T}^{P}+\sum_{l\sigma=2s}\boldsymbol{\rm \pi}^{(T,l\sigma)}\odot\boldsymbol{V}^{(l\sigma)}+k_{B}T\,\tilde{\boldsymbol{\nabla}}\boldsymbol{\cdot}\boldsymbol{\mu}^{TT}\Big\}\Delta t, \end{align}

(3.17b) \begin{align} &\boldsymbol{e}(t+\Delta t) =\boldsymbol{e}(t)+\Delta\hat{\boldsymbol{e}}\nonumber\\ &+\Big\{\boldsymbol{\varOmega}^{\mathcal{A}}+\boldsymbol{\mu}^{RT}\boldsymbol{\cdot}\boldsymbol{F}^{P}+\boldsymbol{\mu}^{RR}\boldsymbol{\cdot}\boldsymbol{T}^{P}+\sum_{l\sigma=2s}\boldsymbol{\rm \pi}^{(R,l\sigma)}\odot\boldsymbol{V}^{(l\sigma)}+k_{B}T\,\tilde{\boldsymbol{\nabla}}\boldsymbol{\cdot}\boldsymbol{\mu}^{RT}\Big\}\Delta t\times\boldsymbol{e}(t), \end{align}

with $\Delta \hat {\boldsymbol {e}}=\Delta \hat {\boldsymbol {\varTheta }}(t)\times \boldsymbol {e}(t)+\tfrac {1}{2} \Delta \hat {\boldsymbol {\varTheta }}(t)\boldsymbol {\cdot }[\boldsymbol {e}(t) \Delta \hat {\boldsymbol {\varTheta }}(t)-\Delta \hat {\boldsymbol {\varTheta }}(t)\boldsymbol {e}(t)]$, while

(3.18a,b)\begin{equation} \left\langle \begin{pmatrix}\Delta\hat{\boldsymbol{R}}\\ \Delta\hat{\boldsymbol{\varTheta}} \end{pmatrix}\right\rangle =\boldsymbol{0},\quad\left\langle \begin{pmatrix}\Delta\hat{\boldsymbol{R}}\\ \Delta\hat{\boldsymbol{\varTheta}} \end{pmatrix}\begin{pmatrix}\Delta\hat{\boldsymbol{R}}\\ \Delta\hat{\boldsymbol{\varTheta}} \end{pmatrix}^{\text{tr}}\right\rangle =2k_{B}T\,\mathbb{M}\,\Delta t.\end{equation}

It is clear that the grand mobility matrix $\mathbb {M}$ and the grand propulsion tensor $\boldsymbol {\varPi }^{(l\sigma )}$ satisfy

(3.19a,b)\begin{equation} \mathbb{M}=\begin{pmatrix}\boldsymbol{\mu}^{TT} & \boldsymbol{\mu}^{TR}\\ \boldsymbol{\mu}^{RT} & \boldsymbol{\mu}^{RR} \end{pmatrix}=\boldsymbol{\varGamma}^{{-}1},\quad\boldsymbol{\varPi}^{(l\sigma)}=\begin{pmatrix} \boldsymbol{\rm \pi}^{(T,l\sigma)}\\ \boldsymbol{\rm \pi}^{(R,l\sigma)} \end{pmatrix}={-}\mathbb{M}\boldsymbol{\cdot}\hat{\boldsymbol{\varGamma}}^{(l\sigma)}. \end{equation}

Onsager–Casimir symmetry implies symmetry of the mobility matrix, and we can identify the so-called propulsion tensors as $\boldsymbol {{\rm \pi} }^{(\alpha,l\sigma )}=-\boldsymbol {\mu }^{\alpha T}\boldsymbol {\cdot }\hat {\boldsymbol {\gamma }}^{(T,l\sigma )}-\boldsymbol {\mu }^{\alpha R}\boldsymbol {\cdot }\hat {\boldsymbol {\gamma }}^{(R,l\sigma )},$ with $\alpha \in \{T,R\}$ (Singh & Adhikari 2018). The convective terms in the update equations constitute the thermal drift, which arises from a simple forward Euler integration scheme of the Langevin equations. The occurring derivative $\tilde {\boldsymbol {\nabla }}$ is the standard spatial gradient (with respect to the source point). If the mobilities depend on the particle orientation, additional orientational convective terms must be included. For the spherical particles considered here, however, these terms do not contribute. The quadratic term in $\Delta \boldsymbol {\varTheta }$ in $\Delta \hat {\boldsymbol {e}}$ is needed to preserve the condition $|\boldsymbol {e}|=1$, as discussed in Makino & Doi (2004) and De Corato et al. (2015).

As the Stokes equation defines a dissipative system, any acceptable approximation of $\mathbb {M}$ must remain positive–definite for all physical configurations, e.g. when a simulated particle does not overlap with nearby boundaries (Cichocki et al. 2000). In Appendix B, starting from the BIE of Stokes flow and using a Galerkin–Jacobi iterative method, we outline how to find such solutions, in principle to arbitrary accuracy in the distance between the particle and surrounding boundaries, for the mobility and propulsion tensors in terms of the Green's function $\boldsymbol {G}$ of Stokes flow. For this, we write the Green's function as the sum (Smoluchowski 1911)

(3.20)\begin{equation} \boldsymbol{G}(\boldsymbol{R},\tilde{\boldsymbol{R}})=\boldsymbol{G}^{o}(\boldsymbol{r})+\boldsymbol{G}^{*}(\boldsymbol{R},\tilde{\boldsymbol{R}}),\end{equation}

where $\boldsymbol {r}=\boldsymbol {R}-\tilde {\boldsymbol {R}}$, and $\boldsymbol {G}^{o}(\boldsymbol {r})=(\boldsymbol {I}+\hat {\boldsymbol {r}}\hat {\boldsymbol {r}})/8{\rm \pi} \eta r$ is the Oseen tensor for unbounded Stokes flow (Oseen 1927; Pozrikidis 1992). The term $\boldsymbol {G}^{*}$ is the correction necessary to satisfy additional boundary conditions in the system. In the unbounded domain, where ${\boldsymbol {G}=\boldsymbol {G}^{o}(\boldsymbol {r})}$, the mobility matrix $\mathbb {M}$ diagonalises and the propulsion tensors vanish identically,

(3.21a–d)\begin{equation} \boldsymbol{\mu}^{TT}=\mu_{T}\,\boldsymbol{I},\quad\boldsymbol{\mu}^{R}=\mu_{R}\,\boldsymbol{I}, \quad\boldsymbol{\mu}^{TR}=\boldsymbol{\mu}^{RT}=\boldsymbol{0},\quad \boldsymbol{\rm \pi}^{(\alpha,l\sigma)}=\boldsymbol{0}. \end{equation}

Here, $\mu _{T}=(6{\rm \pi} \eta b)^{-1}$ and $\mu _{R}=(8{\rm \pi} \eta b^{3})^{-1}$ are the well-known mobility coefficients for translation and rotation of a sphere of radius $b$ in an unbounded fluid of viscosity $\eta$ (Stokes 1850). For a system containing additional boundaries, we obtain corrections to the above expressions in terms of derivatives of $\boldsymbol {G}^{*}$. As shown in the Appendix, to leading order in the Jacobi iteration the mobilities are

(3.22a–c)\begin{align} \boldsymbol{\mu}^{TT}=\mu_{T}\left(\boldsymbol{I}+6{\rm \pi}\eta b\,\mathcal{F}^{0}\tilde{\mathcal{F}}^{0}\boldsymbol{G}^{*}\right),\quad \boldsymbol{\mu}^{TR}=\tfrac{1}{2}\mathcal{F}^{0}\,\tilde{\boldsymbol{\nabla}}\times\boldsymbol{G}^{*}, \quad\boldsymbol{\mu}^{RR}=\mu_{R}\,\boldsymbol{I}+\tfrac{1}{2}\boldsymbol{\nabla}\times\boldsymbol{\mu}^{TR}, \end{align}

where we have defined the differential operators $\mathcal{F}^{l}=\big(1+b^{2}/(4l+6)\nabla^{2}\big)$ and $\tilde{\mathcal{F}}^{l}=\big(1+b^{2}/(4l+6)\tilde{\nabla}^{2}\big)$.

Governed by the particle's activity, we choose to retain the leading symmetric and polar modes of the slip. As demonstrated in the next section, this requires the following propulsion tensors:

(3.23a–c) \begin{align} \boldsymbol{\rm \pi}^{(T,2s)}=\tfrac{10{\rm \pi}\eta b^{2}}{3}\mathcal{F}^{0}\tilde{\mathcal{F}}^{1}\left[\tilde{\boldsymbol{\nabla}}\boldsymbol{G}^{*}\right.&\left.+(\tilde{\boldsymbol{\nabla}}\boldsymbol{G}^{*})^{\text{tr}}\right],\qquad\boldsymbol{\rm \pi}^{(T,3t)}=-\tfrac{2{\rm \pi}\eta b^{3}}{5}\mathcal{F}^{0}\tilde{\nabla}^{2}\boldsymbol{G}^{*},\nonumber\\[1em] &\boldsymbol{\rm \pi}^{(T,4t)}=-\tfrac{2{\rm \pi}\eta b^{4}}{63}\mathcal{F}^{0}\tilde{\boldsymbol{\nabla}}\tilde{\nabla}^{2}\boldsymbol{G}^{*},\end{align}

given to leading order in the Jacobi iteration. The structure of the problem implies that $\boldsymbol {{\rm \pi} }^{(R,l\sigma )}=\tfrac {1}{2}(\boldsymbol {\nabla }\times \boldsymbol {{\rm \pi} }^{(T,l\sigma )}).$ To the given order these have been first obtained in Singh & Adhikari (2018).

3.3. Chemo-hydrodynamic coupling and resulting propulsion

We now consider the boundary condition (IIl), coupling the hydrodynamic to the chemical problem. We observe that the differential operator $\boldsymbol {\chi }$ defined in (Ig) implies tangential slip such that $\hat {\boldsymbol {b}}\boldsymbol {\cdot }\boldsymbol {v}^{\mathcal {A}}=0$, i.e. chemical gradients at the surface of the particle can only drive tangential slip flows. Satisfying this condition, we write the tangential modes in the expansion of the slip in (3.10a,b) with a subscript $s$ as $\boldsymbol {V}_{s}^{(l\sigma )}$. In order to obey the tangential slip condition, the symmetric and trace modes of the slip expansion coefficients have to satisfy

(3.24)\begin{equation} \boldsymbol{V}_{s}^{([l+2]t)}={-}l(2l+3)\boldsymbol{V}_{s}^{(ls)}.\end{equation}

This means that, whenever a $\boldsymbol {V}_{s}^{(ls)}$ mode is generated, a $\boldsymbol {V}_{s}^{([l+2]t)}$ mode of strength given by (3.24) will be generated too. For the anti-symmetric modes $\boldsymbol {V}_{s}^{(la)}$ there is no such condition as they produce tangential slip flow by definition (Singh et al. 2015).

Finally, to express the boundary condition (Ig) in a basis of TSH, we expand the phoretic mobility as

(3.25)\begin{equation} \mu_{c}(\boldsymbol{b})=\sum_{q=0}^{\infty}\tilde{w}_{q}\boldsymbol{M}^{(q)}\odot\boldsymbol{Y}^{(q)}(\hat{\boldsymbol{b}}).\end{equation}

The coefficients $\boldsymbol {M}^{(q)}$ are symmetric and traceless tensors of rank-$q$. This yields the linear system corresponding to (IIl)

(3.26)\begin{equation} \boldsymbol{V}_{s}^{(l\sigma)} =\boldsymbol{\chi}^{(l\sigma,q)}\odot\boldsymbol{C}^{(q)}.\end{equation}

The coupling tensor $\boldsymbol {\chi }^{(l\sigma,q)}$ is given in Appendix A.3, and satisfies $\boldsymbol {\chi }^{(l\sigma,0)}=\boldsymbol {0}$, i.e. a uniform surface concentration does not induce slip.

In principle, any form of tangential slip can be generated by the chemo-hydrodynamic coupling in (3.26). Here, we only consider the leading polar ($\boldsymbol {V}_{s}^{(3t)}$), chiral ($\boldsymbol {V}_{s}^{(2a)}$) and symmetric ($\boldsymbol {V}_{s}^{(2s)}$) modes. Using (3.24), we can identify

(3.27a–c)\begin{equation} \boldsymbol{V}^{\mathcal{A}}={-}\boldsymbol{V}_{s}^{(1s)}=\tfrac{1}{5}\boldsymbol{V}_{s}^{(3t)},\quad 2b\boldsymbol{\varOmega}^{\mathcal{A}}={-}\boldsymbol{V}_{s}^{(2a)},\quad \boldsymbol{V}_{s}^{(2s)}={-}\tfrac{1}{14}\boldsymbol{V}_{s}^{(4t)}.\end{equation}

In the following, we therefore parametrise polar, chiral and symmetric slip by $\boldsymbol {V}^{\mathcal {A}}$, $\boldsymbol {\varOmega }^{\mathcal {A}}$ and $\boldsymbol {V}_{s}^{(2s)}$, respectively. With this, the propulsion terms in the update equations (3.17) are

(3.28)\begin{equation} \sum_{l\sigma=2s}\boldsymbol{\varPi}^{(l\sigma)}\odot\boldsymbol{V}^{(l\sigma)}=5\begin{pmatrix} \boldsymbol{\rm \pi}^{(T,3t)}\\ \boldsymbol{\rm \pi}^{(R,3t)} \end{pmatrix}\boldsymbol{\cdot}\boldsymbol{V}^{\mathcal{A}}+\begin{pmatrix}\boldsymbol{\rm \pi}^{(T,2s)}-14\boldsymbol{\rm \pi}^{(T,4t)}\\ \boldsymbol{\rm \pi}^{(R,2s)}-14\boldsymbol{\rm \pi}^{(R,4t)} \end{pmatrix}\boldsymbol{\colon}\boldsymbol{V}_{s}^{(2s)},\end{equation}

for an autophoretic particle, and (3.26) yields

(3.29a)$$\begin{gather} \boldsymbol{V}^{\mathcal{A}} ={-}\frac{1}{4{\rm \pi} b^{3}}\sum_{q=1}^{\infty}\left[\frac{q+1}{2q+1} \boldsymbol{M}^{(q-1)}-q(q+1)\boldsymbol{M}^{(q+1)}\right]\odot\boldsymbol{C}^{(q)}, \end{gather}$$

(3.29b)$$\begin{gather}\boldsymbol{\varOmega}^{\mathcal{A}} ={-}\frac{3}{8{\rm \pi} b^{4}}\sum_{q=1}^{\infty}q \boldsymbol{M}^{(q)}\times'\boldsymbol{C}^{(q)}, \end{gather}$$

(3.29c)$$\begin{gather}\boldsymbol{V}_{s}^{(2s)} =\frac{3}{4{\rm \pi} b^{3}}\sum_{q=1}^{\infty}\left[\frac{q+1}{4q^{2}-1} \boldsymbol{M}^{(q-2)}+\frac{3q}{2q+3}\boldsymbol{M}_{{sym}}^{(q)}-q(q+1)(q+2) \boldsymbol{M}^{(q+2)}\right]\odot\boldsymbol{C}^{(q)}. \end{gather}$$

For brevity, we have left the solution for $\boldsymbol {C}^{(q)}$ of (3.4) implicit. Here, we have defined a cross-product for irreducible tensors as $(\boldsymbol {M}^{(q)}\times '\boldsymbol {C}^{(q)}){}_{\alpha }=\epsilon _{\alpha \beta \gamma }M_{\beta ({\scriptscriptstyle Q-1})}^{(q)}C_{\gamma ({\scriptscriptstyle Q-1})}^{(q)}$ and a symmetric and traceless product contracting $(q-1)$ indices, $(\boldsymbol {M}_{{sym}}^{(q)}\odot \boldsymbol {C}^{(q)})_{\alpha \beta }=\varDelta _{\alpha \beta,\alpha '\beta '}^{(2)}M_{\alpha '({\scriptscriptstyle Q-1})}^{(q)}C_{\beta '({\scriptscriptstyle Q-1})}^{(q)}$, where we have used the short-hand notation $Q=\gamma _{1}\gamma _{2}\dots \gamma _{q}$ for Cartesian indices (Damour & Iyer 1991).

4. Applications

In this section, we demonstrate the methodology introduced thus far with the help of three examples. First, we discuss model design to achieve certain types of motion and the effect of thermal fluctuations in the bulk fluid. With the help of the appropriate Green's functions, we then provide the chemical and hydrodynamic connectors necessary to describe the dynamics of a phoretic particle near a plane, chemically permeable surface of two immiscible liquids. In a representative example, we investigate some of the chemo-hydrodynamic effects this interface has on the motion of a self-rotating autophoretic particle. Finally, we discuss the hovering state of an isotropic chemical source particle above an interface as a function of particle activity, and the chemo-hydrodynamic properties of the interface.

We use the following notation for the uniaxial parameterisation of the $q$th modes of the phoretic mobility and surface flux:

(4.1a,b)\begin{equation} \boldsymbol{M}^{(q)}=M_{q}\boldsymbol{Y}^{(q)}(\boldsymbol{p}_{q}),\quad \boldsymbol{J}^{(q)}=J_{q}\boldsymbol{Y}^{(q)}(\boldsymbol{e}_{q}).\end{equation}

Here, $M_{q}$ and $J_{q}$ are constants representing the strength of the $q$th mode, while $\boldsymbol {p}_{q}$ and $\boldsymbol {e}_{q}$ are unit vectors.

4.1. Programmed Brownian motion in the bulk

In the bulk, far away from boundaries, we can simplify (3.17). First, we non-dimensionalise the equations by rescaling velocities by the speed of a particle with constant phoretic mobility, namely $4{\rm \pi} b^{2}\,V=\mu _{c}J_{1}/D$. Angular velocities are rescaled by $V/b$. Renaming rescaled variables such that (3.16a,b) reads the same, we obtain

(4.2a,b)\begin{equation} \partial_{t}\boldsymbol{R}=\boldsymbol{V}^{\mathcal{A}}+\sqrt{2D_{T}}\,\boldsymbol{\xi}_{T}, \quad\partial_{t}\boldsymbol{e}=\left(\boldsymbol{\varOmega}^{\mathcal{A}}+\sqrt{2D_{R}}\, \boldsymbol{\xi}_{R}\right)\times\boldsymbol{e},\end{equation}

where for a spherical body in an unbounded fluid the translational diffusivity $D_{T}=\mathcal {B}/6$ and the rotational diffusivity $D_{R}=\mathcal {B}/8$ are isotropic and defined in terms of the Brown number

(4.3)\begin{equation} \mathcal{B}=\frac{k_{B}T}{{\rm \pi}\eta b^{2}V}, \end{equation}

the ratio of Brownian to hydrodynamic forces. Analogously, a particle Péclet number can be defined by ${Pe}=1/\mathcal {B}$ (Mozaffari et al. 2018). For a model including modes up to second order in both the phoretic mobility and the flux expansion, the velocity and angular velocity read

(4.4a)$$\begin{gather} \boldsymbol{V}^{\mathcal{A}} ={-}\boldsymbol{e}_{1}+3m_{2}\left[3\boldsymbol{p}_{2}\left(\boldsymbol{p}_{2} \boldsymbol{\cdot}\boldsymbol{e}_{1}\right)-\boldsymbol{e}_{1}\right]-6m_{1}j_{2} \left[3\boldsymbol{e}_{2}\left(\boldsymbol{e}_{2}\boldsymbol{\cdot}\boldsymbol{p}_{1}\right)- \boldsymbol{p}_{1}\right], \end{gather}$$

(4.4b)$$\begin{gather}\boldsymbol{\varOmega}^{\mathcal{A}} ={-}\tfrac{9}{4}m_{1} \left(\boldsymbol{p}_{1}\times\boldsymbol{e}_{1}\right)-270m_{2}j_{2}\left(\boldsymbol{p}_{2} \boldsymbol{\cdot}\boldsymbol{e}_{2}\right)\left(\boldsymbol{p}_{2}\times\boldsymbol{e}_{2}\right), \end{gather}$$

where $m_{i}=M_{i}/M_{0}$ and $j_{2}=J_{2}/J_{1}$. As will be convenient later, we define the angles $\boldsymbol {e}_{1}\boldsymbol {\cdot }\boldsymbol {p}_{i}=\cos \alpha _{i}$, $\boldsymbol {e}_{1}\boldsymbol {\cdot }\boldsymbol {e}_{2}=\cos \beta$ and $\boldsymbol {e}_{2}\boldsymbol {\cdot }\boldsymbol {p}_{2}=\cos \gamma$. Without loss of generality, we choose $\boldsymbol {e}_{1}$ as the orientation of the particle. This constitutes a minimal model capable of modelling the five distinct types of motion (Lisicki et al. 2018): (i) pure translation, (ii) pure rotation (spinning), (iii) parallel rotation and translation, (iv) perpendicular rotation and translation (circular swimming) and (v) helical motion. In the following we briefly discuss particle designs for each type of motion and analyse how thermal noise affects the dynamics by computing the mean-squared displacement (MSD) $\left\langle \Delta\boldsymbol{}{r}^{2}(t)\right\rangle = \left\langle [\boldsymbol{r}(t)-\boldsymbol{r}(0)]^2\right\rangle$ of selected examples.

Pure translation is the simplest kind of motion and is achieved by choosing $m_{1}=m_{2}=j_{2}=0$. The update equations are those of an active Brownian particle (ABP) with swimming direction $-\boldsymbol {e}_{1}$

(4.5a)$$\begin{gather} \boldsymbol{R}(t+\Delta t) =\boldsymbol{R}(t)-\boldsymbol{e}_{1}\Delta t+\sqrt{2D_{T}}\Delta\boldsymbol{W}_{T}, \end{gather}$$

(4.5b)$$\begin{gather}\boldsymbol{e}_{1}(t+\Delta t) =\boldsymbol{e}_{1}(t)+\sqrt{2D_{R}}\Delta\boldsymbol{W}_{R}\times\boldsymbol{e}_{1}(t), \end{gather}$$

where $\Delta \boldsymbol {W}_{T}$ and $\Delta \boldsymbol {W}_{R}$ are increments of mutually independent Wiener processes (Gardiner 1985). In figure 1(a), we show the simulated (markers) and theoretical (dashed lines) MSDs for such a particle at various temperatures. The MSD for an ABP is known exactly and is given by (Fodor & Marchetti 2018)

(4.6)\begin{equation} \left\langle \Delta\boldsymbol{r}^{2}(t)\right\rangle _{{tra}}=6\left(D_{T}+D_{A}\right)t+2(V\tau)^{2}\left({\rm e}^{{-}t/\tau}-1\right), \end{equation}

where $D_{A}=\mu _{T}V^{2}\tau /3$ is the active diffusion coefficient and $\tau ^{-1}=2D_{R}$ is the persistence time due to rotational noise. The persistence time indicates a transition from a ballistic to a diffusive dynamics, clearly visible in the figure. In the limit of zero temperature, i.e. $\mathcal {B}=0$, the MSD reduces to

(4.7)\begin{equation} \left\langle \Delta\boldsymbol{r}^{2}(t)\right\rangle _{{tra}}^{\mathcal{B}=0}=(Vt)^{2}, \end{equation}

as indicated in the figure.

Figure 1. Mean-squared displacement of the programmed Brownian motion of an autophoretic particle. In panels (a–c) we compare the non-dimensionalised MSDs of translational, circular and helical swimming computed for Brownian simulations with their theoretical predictions (see the main text for the latter) at various temperatures characterised by a particle Péclet number ${Pe}$. We define ${Pe}=\infty$ (no noise), ${Pe}=100$ (moderate noise) and ${Pe}=10$ (strong noise) following Mozaffari et al. (2018). For all three types of motion a diffusive regime ${\sim }t$ can be identified above a certain persistence time broadly determined by the amount of rotational diffusivity. The insets show the respective trajectories over an arbitrary time $T$ in the limit of zero temperature.

A spinning particle (pure rotation) can be modelled by choosing $\alpha _{2}={\rm \pi} /2$, $m_{2}=-1/3$ and $j_{2}=0$ while $\sin \alpha _{1}\neq 0$. This is captured by the update equations

(4.8a)$$\begin{gather} \boldsymbol{R}(t+\Delta t) =\boldsymbol{R}(t)+\sqrt{2D_{T}}\Delta\boldsymbol{W}_{T}, \end{gather}$$

(4.8b)$$\begin{gather}\boldsymbol{e}_{1}(t+\Delta t) =\boldsymbol{e}_{1}(t)+\tfrac{9}{4}m_{1}\left(\boldsymbol{p}_{1}(t)-\boldsymbol{e}_{1}(t) \cos\alpha_{1}\right)\Delta t+\sqrt{2D_{R}}\Delta\boldsymbol{W}_{R}\times\boldsymbol{e}_{1}(t). \end{gather}$$

Additional translation parallel to rotation, on the other hand, occurs for the parameter values $\alpha _{1}=0$, $\alpha _{2}=\beta ={\rm \pi} /2$ and $\sin (2\gamma )\neq 0$, with $\boldsymbol {e}_{1}$ as the translation and rotation axis of the update equations

(4.9a)$$\begin{gather} \boldsymbol{R}(t+\Delta t) =\boldsymbol{R}(t)-\boldsymbol{e}_{1}\left[1+3m_{2}-6m_{1}j_{2}\right]\Delta t+ \sqrt{2D_{T}}\Delta\boldsymbol{W}_{T}, \end{gather}$$

(4.9b)$$\begin{gather}\boldsymbol{e}_{1}(t+\Delta t) =\boldsymbol{e}_{1}(t)+\sqrt{2D_{R}}\Delta\boldsymbol{W}_{R}\times \boldsymbol{e}_{1}(t). \end{gather}$$

Circular swimming (perpendicular rotation and translation) is obtained by choosing $m_{2}=j_{2}=0$ and $\sin \alpha _{1}\neq 0$. For such a self-rotating circle swimmer one can compute the MSD exactly if the Brownian motion is confined to the plane perpendicular to $\boldsymbol {\varOmega }^{\mathcal {A}}$ (van Teeffelen & Löwen 2008)

(4.10)\begin{align} &\left\langle \Delta\boldsymbol{r}^{2}(t)\right\rangle _{{circ}}\nonumber\\ &\quad =2\lambda^{2}\left\{ \varOmega^{2}-D_{R}^{2}+D_{R}(D_{R}^{2}+\varOmega^{2})t+{\rm e}^{{-}D_{R}t}\left[(D_{R}^{2}-\varOmega^{2})\cos\varOmega t-2D_{R}\varOmega\sin\varOmega t\right]\right\}\nonumber\\ &\qquad +4D_{T}t, \end{align}

where $\lambda =V/(D_{R}^{2}+\varOmega ^{2})$. Here, $\varOmega$ represents the circular frequency. In the limit of zero temperature, this reduces to

(4.11)\begin{equation} \left\langle \Delta\boldsymbol{r}^{2}(t)\right\rangle _{{circ}}^{\mathcal{B}=0}=2\left(\frac{V}{\varOmega}\right)^{2}(1-\cos\varOmega t), \end{equation}

where $V/\varOmega$ is the radius of the circular motion. Since in this case we can restrict our attention to the $x$–$z$ plane, we can define the planar polar angle $\vartheta$ such that $\boldsymbol {e}_{1}=\cos \vartheta \,\hat {\boldsymbol {x}}+\sin \vartheta \,\hat {\boldsymbol {z}}.$ The update equations then take the simplified form

(4.12a)$$\begin{gather} x(t+\Delta t) =x(t)-\cos\vartheta\Delta t+\sqrt{2D_{T}}\Delta W_{x}, \end{gather}$$

(4.12b)$$\begin{gather}z(t+\Delta t) =z(t)-\sin\vartheta\Delta t+\sqrt{2D_{T}}\Delta W_{z}, \end{gather}$$

(4.12c)$$\begin{gather}\vartheta(t+\Delta t) =\vartheta(t)+\tfrac{9}{4}m_{1}\sin\alpha_{1}\Delta t+\sqrt{2D_{R}}\Delta W_{\vartheta}, \end{gather}$$

where $\Delta W_{x}$, $\Delta W_{z}$ and $\Delta W_{\vartheta }$ are increments of mutually independent Wiener processes. In figure 1(b) we compare the MSD obtained from simulations with the theoretical expressions.

Helical motion occurs for all other parameter values, representing a general non-axisymmetric phoretic particle. A simple example is given by choosing $j_{2}=0$ and $\cos \beta \neq 0$. The corresponding update equations are

(4.13a)$$\begin{gather} \boldsymbol{R}(t+\Delta t) =\boldsymbol{R}(t)+\left(9m_{2}\boldsymbol{p}_{2}\cos\beta-\boldsymbol{e}_{1} \left(1+3m_{2}\right)\right)\Delta t+\sqrt{2D_{T}}\Delta\boldsymbol{W}_{T}, \end{gather}$$

(4.13b)$$\begin{gather}\boldsymbol{e}_{1}(t+\Delta t) =\boldsymbol{e}_{1}(t)+\tfrac{9}{4}m_{1}\left(\boldsymbol{p}_{1}(t)- \boldsymbol{e}_{1}(t)\cos\alpha_{1}\right)\Delta t+\sqrt{2D_{R}}\Delta\boldsymbol{W}_{R}\times\boldsymbol{e}_{1}(t). \end{gather}$$

At zero temperature, the pitch angle $\psi$ of the resulting helix is given by the simple expression

(4.14)\begin{equation} \frac{\boldsymbol{V}^{\mathcal{A}}}{\left|\boldsymbol{V}^{\mathcal{A}}\right|}\boldsymbol{\cdot} \frac{\boldsymbol{\varOmega}^{\mathcal{A}}}{\left|\boldsymbol{\varOmega}^{\mathcal{A}}\right|}={-}\frac{9m_{2}\sin(2\beta)}{2\sqrt{1+m_{2}(6-\cos^{2}\beta)+3m_{2}^{2}(3-26\cos^{2}\beta)}}=\cos\psi. \end{equation}

In figure 1(c) we approximate the corresponding MSD by a superposition of translational and circular Brownian motion discussed in the previous paragraphs such that

(4.15)\begin{equation} \left\langle \Delta\boldsymbol{r}^{2}(t)\right\rangle _{{hel}}=\left\langle \Delta\boldsymbol{r}^{2}(t)\right\rangle _{{tra}}^{V_{{\perp}}}+\left\langle \Delta\boldsymbol{r}^{2}(t)\right\rangle _{{circ}}^{V_{{\parallel}}}-4D_{T}t, \end{equation}

where the last term is introduced to avoid accounting for translational noise twice. The superscripts of the MSD terms indicate which component of the velocity enters the respective terms. Here, $V_{\perp }$ is the component of the velocity perpendicular to the plane of the circular motion, and $V_{\parallel }$ the component within that plane. Naturally, in the limit of zero temperature this reduces to the exact deterministic MSD for a helix

(4.16)\begin{equation} \left\langle \Delta\boldsymbol{r}^{2}(t)\right\rangle_{{hel}}^{{det}}=(V_{{\perp}}t)^{2}+2\left(\frac{V_{{\parallel}}}{\varOmega} \right)^{2}(1-\cos\varOmega t). \end{equation}

There is good agreement between this approximation and the MSD computed from simulated Brownian trajectories. Compared with the MSD of an ABP in figure 1(a), a kink indicating the period $2{\rm \pi} /\varOmega$ of the circular part of the motion is clearly visible.

4.2. Autophoresis near a permeable interface

We now introduce a plane surface of two immiscible liquids of viscosity ratio $\lambda ^{f}=\eta _{2}/\eta _{1}$ and solute diffusivity ratio $\lambda ^{c}=D_{2}/D_{1}$ in the vicinity of the particle, see figure 2. The interface is characterised by the Green's functions in table 3. Here, $H$ satisfies the boundary condition of continuous normal flux $\hat {\boldsymbol {z}}\boldsymbol {\cdot }\boldsymbol {j}^{(1)}=\hat {\boldsymbol {z}} \boldsymbol {\cdot }\boldsymbol {j}^{(2)}$ across the interface, and $\boldsymbol {G}$ arises from the boundary conditions of continuous tangential flow $v_{\rho }^{(1)}=v_{\rho }^{(2)},$ vanishing normal flow $v_{z}^{(1)}=v_{z}^{(2)}=0$ and continuous tangential stress $\eta _{1}\sigma _{\rho z}^{(1)}=\eta _{2}\sigma _{\rho z}^{(2)}$ across the interface, where the index $\rho =x,y$ lies in the plane of the interface (Jones, Felderhof & Deutch 1975; Aderogba & Blake 1978). The superscripts label whether the quantity of interest is above or below the interface, where $(1)$ refers to the positive half-space $z>0$.

Figure 2. Half-coated phoretic particle near the surface of two immiscible liquids. Schematic representation of a half-coated phoretic particle (Janus particle) of radius $b$ and with orientation $\boldsymbol {e}_{1}$ at a height $h$ above an interface of viscosity ratio $\lambda ^{f}=\eta _{2}/\eta _{1}$ and diffusivity ratio $\lambda ^{c}=D_{2}/D_{1}$. The latter effectively measures how permeable the interface is to the solutes, where $\lambda ^{c}=0$ implies an impermeable and $\lambda ^{c}=1$ a perfectly permeable interface. Due to the cylindrical symmetry of the system we can restrict our attention to the $x$–$z$ plane, where the symbols $\parallel$ and $\perp$ imply motion parallel and perpendicular to the interface, respectively.

Table 3. Green's functions for a plane interface. The Green's functions for the concentration (Laplace equation, left) and velocity (Stokes equation, right) fields near a plane, chemically permeable fluid–fluid interface of solute diffusivity ratio $\lambda ^{c}=D_{2}/D_{1}$ and viscosity ratio $\lambda ^{f}=\eta _{2}/\eta _{1}$ at $z=0$ are given. The particle is in the positive half-space $z>0$, where the solute diffusivity is $D_{1}$ and the fluid viscosity is $\eta _{1}$. We use $\boldsymbol {r}=\boldsymbol {R}-\tilde {\boldsymbol {R}}$ and $\boldsymbol {r}^{*}=\boldsymbol {R}-\tilde {\boldsymbol {R}}^{*}$, with the field point $\boldsymbol {R}=(x,y,z)^{\textrm {tr}}$ , the source point $\tilde {\boldsymbol {R}}=(\tilde {x},\tilde {y},\tilde {z})^{\textrm {tr}}$ and the relation $\tilde {\boldsymbol {R}}^{*}=\boldsymbol {\mathcal {M}}\boldsymbol {\cdot }\tilde {\boldsymbol {R}}$ between the physical position vector $\tilde {\boldsymbol {R}}$ and the position vector of the image singularities $\tilde {\boldsymbol {R}}^{*}$. The height of the centre of the particle above the interface is $\tilde {z}=h$. With this, we have $\tilde {\boldsymbol {R}}-\tilde {\boldsymbol {R}}^{*}=(0,0,2h)^{\textrm {tr}}.$ We define $\boldsymbol {\nabla }^{*}=\boldsymbol {\nabla }_{\boldsymbol {r}^{*}}$, $\mathcal {M}_{\beta \gamma }^{f}=(({1-\lambda ^{f}})/({1+\lambda ^{f}}))\delta _{\beta \rho }\delta _{\rho \gamma }-\delta _{\beta z}\delta _{z\gamma }$, the mirroring operator $\mathcal {M}_{\beta \gamma }=\delta _{\beta \rho }\delta _{\rho \gamma }-\delta _{\beta z}\delta _{z\gamma }$ and the index $\rho =x,y$ in the plane of the interface.

To discuss the chemo-hydrodynamic effect a plane interface has on autophoresis and Brownian motion, we choose a simple non-axisymmetric particle model. We truncate the expansions of the phoretic mobility and surface flux each at linear order and choose $J_{0}/3=J_{1}=J$, so that the particle has one inert pole ($j^{\mathcal {A}}=0$) and one active pole ($j^{\mathcal {A}}>0$ ($j^{\mathcal {A}}<0$) for $J>0$ ($J<0$), corresponding to a source (sink) of chemical reactants). This is a first-order approximation to a Janus swimmer so that for $J>0$ we have an inert-side-forward swimmer. We define $\boldsymbol {e}_{1}\boldsymbol {\cdot }\boldsymbol {p}_{1}=\cos \alpha$, where as before, $\boldsymbol {e}_{1}$ shall serve as the orientation of the particle. For $\sin \alpha \neq 0$ the particle is capable of phoretic self-rotation, see the discussion on circular swimming in the previous section.

We choose to truncate the generated concentration field at the surface of the particle at second order with the coefficients determined by (3.4). Since slip is proportional to gradients in the surface concentration we can ignore the terms $\boldsymbol {\zeta }^{(0,q)}$, $\boldsymbol {\zeta }^{(q,0)}$ and $\boldsymbol {\mathcal {E}}^{(0,q)}$. Cylindrical symmetry of the system can then be used to write the remaining non-zero chemical tensors in terms of scalar coefficients, which are given in table 4. For simplicity, we assume the absence of any background concentration field.

Table 4. Chemical coefficients. Scalar coefficients for the linear response to a background concentration field and for the elastance tensors with $\varLambda ^{c}=(1-\lambda ^{c})/(1+\lambda ^{c})$, where $\lambda ^{c}=D_{2}/D_{1}$, and the height above the interface $\hat {h}=h/b$. We have used the exact unbounded coefficients in (3.6a,b) such that $\zeta _{1}=3/2$, $\mathcal {E}_{1}=3/8{\rm \pi} bD_{1}$ and $\mathcal {E}_{2}=5/2{\rm \pi} bD_{1}$. Cylindrical symmetry of the system implies, for the coupling to a linear background concentration field, $\boldsymbol {\zeta }^{(1,1)}=(\boldsymbol {I}-\hat {\boldsymbol {z}}\hat {\boldsymbol {z}})\,\zeta _{\parallel }^{(1,1)}+\hat {\boldsymbol {z}}\hat {\boldsymbol {z}}\, \zeta _{\perp }^{(1,1)}$. The relevant elastance tensors are $\boldsymbol {\mathcal {E}}^{(1,0)}=\hat {\boldsymbol {z}}\,\mathcal {E}^{(1,0)},$ $\boldsymbol {\mathcal {E}}^{(1,1)}=(\boldsymbol {I}-\hat {\boldsymbol {z}}\hat {\boldsymbol {z}})\, \mathcal {E}_{\parallel }^{(1,1)}+\hat {\boldsymbol {z}}\hat {\boldsymbol {z}}\,\mathcal {E}_{\perp }^{(1,1)},$ $\boldsymbol {\mathcal {E}}^{(2,0)}=-(3\hat {\boldsymbol {z}}\hat {\boldsymbol {z}}-\boldsymbol {I})\,\mathcal {E}^{(2,0)}$ and $\mathcal {E}_{\alpha \beta \gamma }^{(2,1)}=\mathcal {E}^{(2,1)}[(\delta _{\gamma \alpha }-\delta _{\gamma z}\delta _{\alpha z})\delta _{\beta z}+(\delta _{\gamma \beta }-\delta _{\gamma z}\delta _{\beta z})\delta _{\alpha z}+(3\delta _{\beta z}\delta _{\alpha z}-\delta _{\alpha \beta })\delta _{\gamma z}]$.

The induced slip sets the surrounding fluid in motion. The fluid then reacts back on the particle and causes rigid-body motion, governed by the equations of motion in (3.17), mediated by mobility and propulsion tensors. Again, the cylindrical symmetry of the system allows us to write the mobility and propulsion tensors in terms of scalar coefficients, summarised in table 5.

Table 5. Hydrodynamic coefficients. Scalar coefficients for the mobility matrices and the relevant propulsion tensors with $\lambda ^{f}=\eta _{2}/\eta _{1}$ and the height above the interface $\hat {h}=h/b$. Cylindrical symmetry of the system allows us to write, for the translational mobilities, $\boldsymbol {\mu }^{TT}\!=(\boldsymbol {I}-\hat {\boldsymbol {z}}\hat {\boldsymbol {z}})\mu _{\parallel }^{TT}+ \hat {\boldsymbol {z}}\hat {\boldsymbol {z}}\mu _{\perp }^{TT}$ and $\boldsymbol {\mu }^{TR}\!=(\boldsymbol {\mu }^{RT})^{\textrm {tr}}\!=\mu ^{TR}\boldsymbol {\epsilon }\boldsymbol {\cdot }\hat {\boldsymbol {z}}$. The mobility $\boldsymbol {\mu }^{RR}$ has the same structure as $\boldsymbol {\mu }^{TT}$ with the corresponding coefficients $\mu _{\parallel }^{RR}$ and $\mu _{\perp }^{RR}$. The propulsion tensor for the leading symmetric slip mode is ${{\rm \pi} _{\alpha \beta \gamma }^{(T,2s)}}={\rm \pi} _{1}^{(T,2s)}[(\delta _{\gamma \alpha }-\delta _{\gamma z}\delta _{\alpha z})\delta _{\beta z}+(\delta _{\beta \alpha }-\delta _{\beta z}\delta _{\alpha z})\delta _{\gamma z}]+{\rm \pi} _{2}^{(T,2s)}(\delta _{\gamma \beta }-3\delta _{\gamma z}\delta _{\beta z})\delta _{\alpha z}$. The propulsion tensor $\boldsymbol {{\rm \pi} }^{(T,4t)}$ is of the same structure as $\boldsymbol {{\rm \pi} }^{(T,2s)}$ with the corresponding coefficients ${\rm \pi} _{1}^{(T,4t)}$ and ${\rm \pi} _{2}^{(T,4t)}$. Since $\boldsymbol {{\rm \pi} }^{(T,3t)}$ has the same structure as $\boldsymbol {\mu }^{TT}$, we adopt an analogous notation with ${\rm \pi} _{\parallel }^{(T,3t)}$ and ${\rm \pi} _{\perp }^{(T,3t)}$. For the propulsion tensors contributing to the particle's rotational dynamics we obtain ${\rm \pi} _{\alpha \beta \gamma }^{(R,2s)}={\rm \pi} ^{(R,2s)}(\delta _{\beta z}\epsilon _{z\gamma \alpha }+\delta _{\gamma z}\epsilon _{z\beta \alpha })$ and $\boldsymbol {{\rm \pi} }^{(R,3t)}={\rm \pi} ^{(R,3t)}\boldsymbol {\epsilon }\boldsymbol {\cdot }\hat {\boldsymbol {z}}$. The propulsion tensor $\boldsymbol {{\rm \pi} }^{(R,4t)}$ is of the same structure as $\boldsymbol {{\rm \pi} }^{(R,2s)}$ with the corresponding coefficient ${\rm \pi} ^{(R,4t)}$.

The full set of mobility coefficients for a wall, a free surface or, indeed, a fluid–fluid interface have been obtained in the literature to a high degree of accuracy (Brenner 1961; Goldman, Cox & Brenner 1967; Felderhof 1976; Lee, Chadwick & Leal 1979; Lee & Leal 1980; Perkins & Jones 1990, 1992). In Appendix C we show that, for the plane boundary, the diffusion terms $\propto \sqrt {2k_{B}T\,\mathbb {M}}$ in (3.17) take a particularly simple analytic form and that, with the coefficients given in table 5, this diffusion matrix is inherently positive–definite for all physical configurations. The propulsion tensors are a unique feature of active particles and have not been obtained in this form in the literature before.

Assuming there is no external torque rotating the particle out of plane, i.e. $\boldsymbol {T}^{P}\boldsymbol {\cdot }\hat {\boldsymbol {z}}\,{=}\,0$, we can once again restrict our attention to the $x$–$z$ plane for which we define the planar polar angle $\vartheta$ such that $\boldsymbol {e}_{1}=\cos \vartheta \,\hat {\boldsymbol {x}}+\sin \vartheta \,\hat {\boldsymbol {z}}$. Autophoretic particles in typical experiments are neither force nor torque free due to mismatches between particle and solvent densities and between gravitational and geometric centres (Drescher et al. 2010; Ebbens & Howse 2010; Palacci et al. 2010, 2013; Buttinoni et al. 2013). Since the resulting forces and torques become dominant, at long distances, over active contributions, it is crucial to include their effects in our analysis. In simulating (3.16a,b) we therefore assume a bottom-heavy Janus particle (the chemically active coating, blue in figures, is assumed to be slightly heavier than the inert side, white in figures). Therefore, we have to take into account gravity in negative $z$-direction and a gravitational torque given by

(4.17a,b)\begin{equation} \boldsymbol{F}^{P}={-}mg\hat{\boldsymbol{z}},\quad\boldsymbol{T}^{P}=\kappa(\hat{\boldsymbol{z}}\times \boldsymbol{e})=\kappa\cos\vartheta\,\hat{\boldsymbol{y}}, \end{equation}

with $m$ the buoyancy-corrected mass of the particle, $g$ the gravitational constant and $\kappa$ a constant parametrising bottom heaviness. Inserting this into the update equations (3.17) with $\boldsymbol {R}=(x,y,h)^{\textrm {tr}}$, where $h(t)$ is the height of the particle above the boundary, we can now simulate the time evolution of this bottom-heavy, non-axisymmetric phoretic particle. In figure 3 we show typical trajectories near a free surface (e.g. an air–water surface) and a fluid–fluid interface of $\lambda ^{f}=50$ (e.g. an oil–water interface). We probe the effect of the nearby boundary on the dynamics of the autophoretic particle by truncating the dynamical system in (3.16a,b) at various orders in $h^{-1}$ and comparing the results; see Appendix C for the truncated expressions.

Figure 3. Chemo-hydrodynamic effects of a boundary. We compare typical trajectories of a bottom-heavy non-axisymmetric active particle (see the main text for the chosen particle specifications) near two types of chemically impermeable ($\lambda ^{c}=0$) liquid–liquid interfaces. In panels (a,b) the interface is chosen to be a free air–water surface ($\lambda ^{f}=0$) and a water–oil interface ($\lambda ^{f}\approx 50$), respectively. The temperature in these panels is zero and the starting point (orientation) of the particle is indicated by a grey disk (line). In the upper panels we compare the trajectories of the particle obtained to various orders in the inverse distance to the wall $h^{-1}$, thus gradually including more interactions with the wall. Here, ‘simulation’ refers to the full dynamical system at the accuracy obtained in this paper. The lower panels show the orientational evolution of the particle in the various approximations in a polar plot, parametrised by the angle $\vartheta$ and the distance to the starting point. There are two main features to be observed when comparing (a,b). First, the particle tends to stay further away from the water–oil interface when compared with the air–water surface. This is expected as the particle's mobility perpendicular to the interface will be reduced with increasing viscosity ratio, see table 5. Second, the approximations are better for higher viscosity ratios. Again, this can be understood intuitively when considering the example that fluid flows produced near a wall decay faster in the far field when compared with fluid flows produced near a free surface (Aderogba & Blake 1978). The slower decay of the latter means that, to achieve the same quality of approximation as with interfaces of higher viscosity ratios, higher orders in $h^{-1}$ are necessary. Panel (c) shows the same system as in panel (b) but with a finite particle Péclet number of ${Pe}=100$, representing experimentally relevant noise levels. It is clear that the translational as well as the rotational diffusion is affected by the presence of the interface. Thermal diffusion also induces a net repulsive effect between the particle and the interface.

At order $h^{-1}$ only hydrodynamic interactions with the boundary due to the gravitational force manifest themselves. It is at the next order, $h^{-2}$, that the gravitational torque and active effects become apparent. The latter comprise hydrodynamic interactions from symmetric propulsion via $\boldsymbol {{\rm \pi} }^{(T,2s)}$ and a purely monopolar chemical interaction with the interface. At this order in the approximation, fore–aft symmetry breaking of the particle is no longer necessary for self-propulsion near a boundary; see Appendix C. An isotropic particle with uniform phoretic mobility $\mu _{c}$ and surface flux $j^{\mathcal {A}}$ will get repelled (attracted) to the interface depending on whether it is a source or sink of chemical reactants and depending on the diffusivity ratio $\lambda ^{c}$ of the interface, see § 4.3 for a detailed discussion. Furthermore, at this order in the approximation the thermal advective term in (3.17) proportional to $\boldsymbol {\nabla }\boldsymbol {\cdot }\mathbb {M}$ starts to affect the dynamics. It is worth noting that, at order $h^{-3}$, our analytical results match those obtained in Ibrahim & Liverpool (2015), using a method of reflections, for a Janus particle of trivial phoretic mobility near an inert no-slip wall.

The system parameters in figure 3 are chosen as follows. The starting position of the particle is at a height $\hat {h}_{0}=h(t=0)/b=2$ and an angle $\vartheta _{0}=-3{\rm \pi} /4$ to the wall. For the surface flux of the particle we choose the dimensionless control parameter $j_{1}=J_{1}/J_{0}=1/3$, modelling a source particle. Its phoretic mobility distribution is specified by the dimensionless parameter $m_{1}=M_{1}/M_{0}=0.7$, implying a significant non-isotropy ($m_{1}=0$ specifies a trivial phoretic mobility). The angle between the axes of surface flux and phoretic mobility is chosen such that $\alpha ={\rm \pi} /2$. In figure 3(c) the Brown number is set to $\mathcal {B}\sim 10^{-2}$, roughly matching a set of experiments on Janus colloids (Jiang, Yoshinaga & Sano 2010; Palacci et al. 2013) with a bead size $b\sim 1\ \mathrm {\mu } \textrm {m}$ and speed $v_{s}\sim 10\ \mathrm {\mu } \textrm {m}\ \textrm {s}^{-1}$. The ratios of gravitational to active forces and torques are chosen such that $mg/F_{A}\sim 10^{-1}$ and $\kappa /T_{A}\sim 10^{-2}$, respectively. Finally, inertial effects decay on the time scale of momentum relaxation, typically $\tau _{T}=m/6{\rm \pi} \eta b$ and $\tau _{R}=I/8{\rm \pi} \eta b^{3}=m/20{\rm \pi} \eta b$ for translational and rotational effects, respectively. The time step $\Delta t$ in our simulation is chosen such that $\tau _{T}/\Delta t\sim 10^{-4}$ and $\tau _{R}/\Delta t\sim 10^{-4}$, ensuring that the Smoluchowski limit of the dynamics provides an appropriate description.

4.3. Hovering above a permeable interface

As mentioned in the previous section, if chemo-hydrodynamic particle–boundary interactions of order $h^{-2}$ and higher are considered, fore–aft symmetry breaking of the particle is no longer necessary for self-propulsion near a boundary. Indeed, self-propulsion of isotropic particles near a boundary has been observed in light-activated phoretic swimmers (Palacci et al. 2013). We therefore consider a particle that is an isotropic source of reactants ($\mu _{c}=\textrm {const.}>0, j^{\mathcal {A}}=\textrm {const.}>0$) and investigate how its dynamics is affected by the viscosity ratio $\lambda^f$ and the diffusivity ratio $\lambda^c$ of the nearby interface in the limit of zero temperature. The particle is assumed to be driven towards the interface due to gravity. With the rotational dynamics and the translational dynamics parallel to the plane vanishing by symmetry, the particle is expected to hover above the interface at a height that can be found by solving $\dot {h}=0$, where $\dot {h}$ is given in Appendix C. Rescaling heights by $b$, mobilities by $\mu _{T}$ and velocities by $\mu _{T}mg$ and renaming the thus non-dimensionalised variables such that they read the same, we obtain the hovering condition

(4.18)\begin{equation} 0={-}\mu_{{\perp}}^{TT}+\tfrac{1}{4}\varLambda^{c}\mathcal{A}_{G}\left[h^{{-}2}(1+5{\rm \pi}_{{\perp}}^{(T,3t)})-3h^{{-}3}({\rm \pi}_{2}^{(T,2s)}-14{\rm \pi}_{2}^{(T,4t)})\right].\end{equation}

Hovering is thus characterised by only one dimensionless number, $\mathcal {A}_{G}=\mu _{c}j^{\mathcal {A}}/D\mu _{T}mg$, defined as the ratio of the speed of a uniformly coated phoretic particle in a uniform concentration gradient $j^{\mathcal {A}}/D$, namely $\mu _{c}j^{\mathcal {A}}/D$, to the settling velocity under gravity $\mu _{T}mg$. This is a measure of activity.

In figure 4, we show how in our approximation the hovering height $h$ of the isotropic particle is affected by its activity, the diffusivity ratio and the viscosity ratio of the boundary. We limit our results to $h>h_{{min}}=1.3$. This is because, very close to the interface, other effects such as electric double-layer and Van der Waals interactions are expected to dominate (Wu & Bevan 2005; Verweij et al. 2020). As expected, a higher chemical diffusivity ratio of the interface, and thus decreased chemical repulsion from it, requires higher particle activities for hovering to remain possible; see figure 5 for an illustration of this effect using the method of images for the concentration field. It is also evident that lower levels of activity are sufficient for hovering above a free surface as compared with a solid wall. This is intuitive when considering that due to increased fluid internal friction flows near a wall decay faster than near a free surface (Aderogba & Blake 1978) and so it is easier for a particle near a free surface to drive flows that make it hover. Using the method of images for Stokes flows this is illustrated in figure 6.

Figure 4. Hovering above an interface. We show the hovering height $h$ (pseudo-colour map) for an active particle at zero temperature as a function of its activity $\mathcal {A}_{G}$ and the diffusivity ratio $\lambda ^{c}$ or the viscosity ratio $\lambda ^{f}$ of the interface. In panel (a) we consider the particle hovering above a wall ($\lambda ^{f}\rightarrow \infty$) that is permeable to the solutes. For a boundary between regions of equal solute diffusivity ($\lambda ^{c}=1$) there exists no chemical repulsion, leading to the particle inevitably crashing into the boundary. Therefore, we only consider values $\lambda ^{c} \, \leq 0.99$. The red line shows the particle's minimum activity to hover at a minimum height of $h_{{min}}=1.3$ as a function of the diffusivity ratio. The white region below the red line indicates physics that is not accessible in our simplified model and we assume that the particle crashes into the boundary due to gravity, i.e. we set $h=1$. The two dashed grey lines indicate the limiting values $\mathcal {A}_{G}^{1}\approx 500$ and $\mathcal {A}_{G}^{2}\approx 10^{5}$ that are required to hover above an impermeable and a highly permeable wall. In panel (b) we consider the particle hovering above an impermeable ($\lambda ^{c}=0$) interface of varying viscosity ratio. The two horizontal dashed grey lines indicate the limiting values $\mathcal {A}_{G}^{3}\approx 3$ and $\mathcal {A}_{G}^{1}$ that are required to hover above a free surface and a solid wall, respectively.

Figure 5. Solute concentration for a particle hovering above a permeable boundary. The chemical field generated by a hovering source particle (green) is shown above ($z>0$) and below ($z<0$) the permeable ($\lambda ^{c}=0.3$) boundary as a pseudo-colour map, where contours indicate regions of constant concentration. We emphasise that the net concentration (right panel) in the region containing the particle is a superposition of the source without the boundary present (a) and its image below the boundary (b). The net concentration below the permeable boundary is then generated by the appropriate boundary conditions. To imply a change in the chemical diffusivity $D_{2}=\lambda ^{c}D_{1}$ the colour map for the region $z<0$ is shown with slightly reduced opacity. Since source particles want to move down chemical gradients (anti-chemotaxis), it is clear that the image creates a repulsive concentration field in $z>0$, making it possible for the particle to hover above the boundary. It is worth noting that, for an impermeable boundary ($\lambda ^{c}=0$), the contour lines meet the boundary at a right angle and the corresponding vector field ($\boldsymbol {\nabla }c$) becomes purely tangential to this ‘no-flux’ boundary. In the limit of rapid solute diffusion the viscosity ratio of the interface has no influence on the chemical field.

Figure 6. Fluid flow for a particle hovering above a fluid–fluid interface. The fluid flow (laboratory frame) generated by a hovering source particle (green) is shown above ($z>0$) and below ($z<0$) a chemically permeable fluid–fluid interface ($\lambda ^{c}=0.3, \lambda ^{f}=10$). The direction of the fluid flow is indicated by black arrows, while its relative magnitude is implied by the overlaid pseudo-colour map. The net flow (c) in the region containing the particle is a superposition of the source without the boundary conditions on the fluid flow and stress (a) and its image below the boundary (b). Note that for the source the chemical boundary conditions must still be satisfied (see the net chemical field in figure 5), inducing a non-trivial slip on the otherwise isotropic particle. The net flow below the interface is then driven by the appropriate velocity and stress boundary conditions. To imply a change in the viscosity $\eta _{2}=\lambda ^{f}\eta _{1}$ the colour map for the region $z<0$ is shown with slightly reduced opacity. It is clear that the image produces an upwards flow in $z>0$ which balances gravity and makes the particle hover away from the interface. In the net flow this creates convection rolls between the particle and the interface which in turn drive convection rolls below the interface.

5. Discussion

In this paper, we have presented a simultaneous solution of the BIEs describing the chemical and the fluid flow around an autophoretic particle in a fluctuating environment. This has been achieved in a basis of TSH. Compared with the common squirmer model approach to active particles (Lighthill 1952; Blake 1971; Pak & Lauga 2014; Pedley, Brumley & Goldstein 2016), our boundary-domain integral method offers the distinct advantage of obtaining the traction on the particle directly in a complete orthonormal basis. This provides a naturally kinetic approach via Newton's equations in which thermal fluctuations manifest themselves as fluctuating stresses. The Brownian motion of an autophoretic particle is obtained in terms of coupled roto-translational stochastic update equations containing mobility and propulsion tensors. The latter are found to arise from chemical activity of the particle and the chemo-hydrodynamic coupling at the particle's surface, inducing a coupling of slip modes. We have obtained exact and leading-order solutions for both the chemical and the fluctuating hydrodynamic problems far away from and in the vicinity of boundaries, respectively. Studying the Brownian motion of particles in the bulk, some of the flexibility of our method in particle design has been demonstrated. In the case of autophoresis near a plane interface, characterised by its solute diffusivity and viscosity ratios, we have provided analytical expressions for the chemo-hydrodynamic coupling tensors. The given mobilities ensure a positive–definite diffusion matrix in stochastic simulations. Finally, we have studied the hovering state of an isotropic phoretic particle above an interface as a function of particle activity, and the diffusivity and viscosity ratios of the interface. In doing so, we have provided numerical results as well as physical insights into the repulsive chemo-hydrodynamic particle–interface interactions.

We have given the leading-order results for the chemical and hydrodynamic coupling tensors. In principle, these can be obtained to arbitrary accuracy, and the general iterative solutions are given in the Appendix. This non-trivial computation will be the topic of future work. While our results in § 3.2 are guaranteed to provide dissipative motion for physical configurations, in Brownian simulations, unphysical situations with a non-zero particle–boundary overlap may occur on occasion (the probability of which can be lowered by imposing a short-range repulsive potential between the particle and the boundary). In this case, one can either impose an ad hoc regularisation on the mobility (Wajnryb et al. 2013; Balboa Usabiaga, Delmotte & Donev 2017; Singh & Adhikari 2017) or use a bounce-back condition, effectively implementing a reflective boundary condition in the simulation (Volpe, Gigan & Volpe 2014).

It is helpful to compare our results with previous work on chemical and hydrodynamic interactions of an active particle in a fluctuating fluid. We have shown that simultaneous harmonic expansions of the surface fields provide a high degree of flexibility in particle design, comparable to previous models capable of motion in three dimensions (Lisicki et al. 2018). Additionally, our framework has been shown to provide a straightforward way of studying the Brownian dynamics of particles that, even in the limit of zero temperature, perform complex motion (van Teeffelen & Löwen 2008; Mozaffari et al. 2018; Bailey et al. 2024). To the best of our knowledge, this is the first work to obtain the elastance for an active particle near an interface separating two fluids of arbitrary diffusivity ratios. The corresponding Green's function, which is given in table 3, does not appear anywhere in the literature, although its derivation is straightforward given the correct boundary conditions. This paper is also the first to simultaneously study the chemo-hydrodynamics of an autophoretic particle near a planar surface of two immiscible fluids of arbitrary ratio of viscosities and diffusivities. While previous works have studied phoretic particles hovering above a chemically impermeable solid wall as a function of particle coverage (Uspal et al. 2015; Ibrahim & Liverpool 2016), we investigated the hovering state as a function of the properties of the interface, relevant for studies on particle aggregation near fluid–fluid or free surfaces (Chen et al. 2015; Hokmabad et al. 2022).

We briefly discuss the level of approximation of explicit results provided in this paper. For a passive particle, the mobility of sedimentation towards a plane interface is known exactly (Brenner 1961). In the absence of an exact solution for motion parallel to a boundary, careful examination of the case when the particle–interface gap distance is much smaller than the particle radius is necessary. In 1967, Goldman et al. (Goldman, Cox & Brenner 1967) used a lubrication approximation to derive an asymptotic solution for this case. However, matching the asymptotic solutions for the near- ($h/b\ll 1$) and far-field ($h/b\gg 1$) limits can be challenging in dynamic simulations (Brady & Bossis 1988; Ichiki 2002). It has been confirmed experimentally that, for parallel motion, the order to which the mobilities are given in this paper provides a good approximation even when the particle–interface gap distance is only a fraction of the particle radius (Choudhury et al. 2017). So while an approach using lubrication theory is appropriate for general motion very close to a plane (Villa et al. 2020, 2023), the given approximation in the mobilities arising from a series expansion can still be expected to be of interest to a wide range of experimental settings in which colloidal particles are studied near a plane boundary. Thus, our work also adds to the existing literature on the mobility (Brenner 1961; Goldman et al. 1967; Felderhof 1976; Lee et al. 1979; Lee & Leal 1980; Perkins & Jones 1990, 1992; Swan & Brady 2007; Michailidou et al. 2009; Daddi-Moussa-Ider et al. 2018) and diffusion (Ermak & McCammon 1978; Wajnryb et al. 2004; Rogers et al. 2012; Delong, Balboa Usabiaga & Donev 2015; Lisicki, Cichocki & Wajnryb 2016) of a sedimenting particle near a boundary. Explicit expressions for propulsion tensors and mobility matrices are given in table 5, while table 4 contains the corresponding chemical connectors near an interface. These will be helpful in Langevin simulations of autophoretic particles in various experimentally realisable settings and for studying fluctuating trajectories of an active particle including both chemical and hydrodynamic interactions.

Aside from its intrinsic theoretical significance, the single-body solution (exact away from boundaries and approximate in complex environments) holds potential value in numerically solving the BIE for multiple particles. This is due to the ability to initiate numerical iterations with the single-body solution. For problems falling under this category, discretised versions of the BIEs result in diagonally dominant linear systems. Notably, the one-body solution serves as the solution in cases where hydrodynamic interactions are disregarded. This implies that starting iterations from the one-body solution can lead to rapid convergence towards diagonally dominant numerical solutions (Singh & Adhikari 2018). In scenarios involving multiple interacting particles, utilising a basis of TSH for expanding surface fields offers distinct advantages over other bases like spherical or vector spherical harmonics, including reduced computational cost due to covariance under rotations (Greengard & Rokhlin 1987; Damour & Iyer 1991; Applequist 2002; Turk 2023). The condition for tangential slip flow in terms of TSH in (3.24) now connects in a straightforward way the formalism for general slip (restricted by mass conservation only) used in previous works (Ghose & Adhikari 2014; Singh et al. 2015; Singh & Adhikari 2018; Singh et al. 2019; Turk et al. 2022) to the present and other problems in which tangential slip is considered, e.g. active drops.

In future work we will analytically and numerically build upon the theoretical results contained in this paper. A detailed analysis of the dynamical system in (3.17) governing autophoresis near an interface might reveal features such as intricate, thermally limited bound states (Mozaffari et al. 2018; Bolitho, Singh & Adhikari 2020) potentially relevant to the study of biofilm formation in bacteria (Wilking et al. 2011; Persat et al. 2015). Removing the assumption of rapid diffusion gives rise to nonlinear advection–diffusion coupling, uncovering a range of potential applications such as the intricate dynamics of active droplets (Herminghaus et al. 2014; Michelin 2023). These and more provide exciting avenues for future research.