Hostname: page-component-745bb68f8f-hvd4g Total loading time: 0 Render date: 2025-02-07T12:59:15.082Z Has data issue: false hasContentIssue false

Artificial chemotaxis of phoretic swimmers: instantaneous and long-time behaviour

Published online by Cambridge University Press:  12 October 2018

Maria Tătulea-Codrean
Affiliation:
Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK
Eric Lauga*
Affiliation:
Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK
*
Email address for correspondence: e.lauga@damtp.cam.ac.uk

Abstract

Phoretic swimmers are a class of artificial active particles that has received significant attention in recent years. By making use of self-generated gradients (e.g. in temperature, electric potential or some chemical product) phoretic swimmers are capable of self-propulsion without the complications of mobile body parts or a controlled external field. Focusing on diffusiophoresis, we quantify in this paper the mechanisms through which phoretic particles may achieve chemotaxis, both at the individual and the non-interacting population level. We first derive a fully analytical law for the instantaneous propulsion and orientation of a phoretic swimmer with general axisymmetric surface properties, in the limit of zero Péclet number and small Damköhler number. We then apply our results to the case of a Janus sphere, one of the most common designs of phoretic swimmers used in experimental studies. We next put forward a novel application of generalised Taylor dispersion theory in order to characterise the long-time behaviour of a population of non-interacting phoretic swimmers. We compare our theoretical results with numerical simulations for the mean drift and anisotropic diffusion of phoretic swimmers in chemical gradients. Our results will help inform the design of phoretic swimmers in future experimental applications.

Type
JFM Papers
Copyright
© 2018 Cambridge University Press 

1 Introduction

The academic community has recently taken significant interest in the better understanding of the locomotion of microorganisms (Lauga & Powers Reference Lauga and Powers2009; Lauga Reference Lauga2016) and the design of biomimetic devices (Nelson, Kaliakatsos & Abbott Reference Nelson, Kaliakatsos and Abbott2010). One popular type of synthetic swimmer is those able to self-propel through autophoresis – a type of design that escapes the technical challenges associated with manufacturing motile body parts on a small scale. This class of devices makes use of self-generated gradients in temperature (thermophoresis), electric potential (electrophoresis) or concentration of some chemical species (diffusiophoresis) in order to induce motion. While currently at the centre of active research, the physical ideas behind autophoresis are by no means novel and a comprehensive overview of the theory of phoretic transport can be found in the classical article by Anderson (Reference Anderson1989).

Recent technological advances have facilitated the manufacturing of diffusiophoretic swimmers in a variety of shapes and chemical properties, starting with pioneering experimental work done by Paxton et al. (Reference Paxton, Kistler, Olmeda, Sen, St. Angelo, Cao, Mallouk, Lammert and Crespi2004) and followed by a suite of detailed experiments (Howse et al. Reference Howse, Jones, Ryan, Gough, Vafabakhsh and Golestanian2007; Ebbens & Howse Reference Ebbens and Howse2011; Ebbens et al. Reference Ebbens, Tu, Howse and Golestanian2012, Reference Ebbens, Gregory, Dunderdale, Howse, Ibrahim, Liverpool and Golestanian2014). There has also been much interest recently in the study of thermophoretic particles (Jiang, Yoshinaga & Sano Reference Jiang, Yoshinaga and Sano2010; Golestanian Reference Golestanian2012; Bickel, Majee & Würger Reference Bickel, Majee and Würger2013) as well as in the understanding of Marangoni self-propelled droplets (Thutupalli, Seemann & Herminghaus Reference Thutupalli, Seemann and Herminghaus2011; Schmitt & Stark Reference Schmitt and Stark2013; Izri et al. Reference Izri, van der Linden, Michelin and Dauchot2014).

In the present paper, we focus on the mechanism of diffusiophoresis due to concentrations of non-ionic species. Many recent papers investigate ways of exploiting geometric asymmetry (Shklyaev, Brady & Córdova-Figueroa Reference Shklyaev, Brady and Córdova-Figueroa2014; Michelin & Lauga Reference Michelin and Lauga2015), asymmetry of surface properties (Golestanian, Liverpool & Ajdari Reference Golestanian, Liverpool and Ajdari2005, Reference Golestanian, Liverpool and Ajdari2007) or a combination of the two (Popescu et al. Reference Popescu, Dietrich, Tasinkevych and Ralston2010) in order to generate a local imbalance of concentration in an otherwise uniform medium and induce motion.

In this paper we quantify the mechanisms through which diffusiophoretic particles may achieve artificial chemotaxis, i.e. the directed motion along an external chemical gradient, both at the individual and the non-interacting population level. It is known that an asymmetric swimmer placed in a uniform background gradient will undergo active reorientation and experience a torque that seeks to align its axis of symmetry with the direction of the gradient (Bickel, Zecua & Würger Reference Bickel, Zecua and Würger2014; Saha, Golestanian & Ramaswamy Reference Saha, Golestanian and Ramaswamy2014). Although the mechanisms through which diffusiophoretic swimmers achieve chemotaxis are qualitatively different from the techniques used by living organisms (see Berg Reference Berg1975), the prospect of achieving the same functionality is highly desirable for both biomedical and technological applications (Nelson et al. Reference Nelson, Kaliakatsos and Abbott2010; Popescu, Tasinkevych & Dietrich Reference Popescu, Tasinkevych and Dietrich2011).

In this work we first quantify the propulsion and reorientation mechanisms associated with the canonical problem of a spherical axisymmetric swimmer placed in a uniform background gradient of solute concentration. We approach this problem using the classical continuum framework of diffusiophoresis (Golestanian et al. Reference Golestanian, Liverpool and Ajdari2007; Jülicher & Prost Reference Jülicher and Prost2009; Sabass & Seifert Reference Sabass and Seifert2012) as opposed to the osmotic framework proposed by Brady and coworkers (Córdova-Figueroa & Brady Reference Córdova-Figueroa and Brady2008; Brady Reference Brady2011; Córdova-Figueroa, Brady & Shklyaev Reference Córdova-Figueroa, Brady and Shklyaev2013). To the best of our knowledge, we use the same set-up as Saha et al. (Reference Saha, Golestanian and Ramaswamy2014), but we generalise and correct their results. Specifically, we derive a fully analytical law for the instantaneous propulsion and orientation of a phoretic swimmer with general axisymmetric surface properties, in the limit of zero Péclet number for both substrate and product. We compute the solution for a weakly reactive swimmer as an expansion in small Damköhler number, including the first-order effects from the substrate which had been neglected by Saha et al. (Reference Saha, Golestanian and Ramaswamy2014). The advantage of having a complete analytical law is that we are able to apply our results to real-life designs of phoretic swimmers, such as the Janus sphere, which could not be achieved using the truncated results previously published.

We next develop a continuum model for the long-time behaviour of a population of non-interacting phoretic swimmers. Most practical applications rely on the synchronised behaviour of a large number of phoretic swimmers (Palacci et al. Reference Palacci, Sacanna, Steinberg, Pine and Chaikin2013, Reference Palacci, Sacanna, Kim, Yi, Pine and Chaikin2014; Wang et al. Reference Wang, Duan, Sen and Mallouk2013) and macroscopic models are essential in understanding how these micro-swimmers would disperse under the competing action of active propulsion and thermal fluctuations. Our approach is inspired by work done on gyrotactic micro-organisms which successfully explained the formation of bioconvection patterns in cells such as Chalmydomonas nivalis. Specifically, we use the framework of generalised Taylor dispersion theory developed by Frankel & Brenner (Reference Frankel and Brenner1991, Reference Frankel and Brenner1993) and apply it to the artificial chemotaxis of phoretic swimmers. We finally compare the predictions of our continuum model with numerical simulations of large samples of swimmers and obtain excellent agreement. Notably, our continuum model demonstrates that the interplay between active propulsion and reorientation leads to the anisotropic diffusion of phoretic swimmers, a piece of physical understanding overlooked by previous studies.

2 Instantaneous behaviour: artificial chemotaxis of an individual particle

In this first part of the paper, we derive the instantaneous law controlling the linear and angular transport of a phoretically active sphere in a solute (or ‘substrate’) gradient.

2.1 Set-up: chemical problem

A spherical particle of radius $R$ is immersed in a uniform gradient of substrate concentration and induces a disturbance to the equilibrium concentration of substrate not only due to the volume that it displaces, but also due to the chemical reaction that it facilitates at its surface. Before introducing the particle, the medium is characterised by a known linear background concentration of substrate, $s_{b}(\boldsymbol{r}^{\prime })=\boldsymbol{r}^{\prime }\boldsymbol{\cdot }\unicode[STIX]{x1D735}s^{\infty }$ , where $\boldsymbol{r}^{\prime }$ is the position vector relative to some location of vanishing substrate concentration. The presence of the device leads to a perturbed concentration of substrate, $s$ , which must be calculated.

Due to a catalytic coating, we assume that the swimmer promotes the conversion of substrate molecules $S$ into product molecules $P$ , after an intermediary stage of binding to the enzyme molecules $E$ . The classical model for this type of chemical behaviour is referred to as Michaelis–Menten kinetics, and starts from the reaction formula

(2.1)

where $k_{1},k_{-1}$ and $k_{2}$ are reaction rates (for a review, see Nelson Reference Nelson2008).

When the supply of substrate is much larger than the supply of enzyme, the reaction is limited by the availability of enzymatic sites, and the system classically reaches a steady non-equilibrium state. The rate of conversion of substrate molecules $S$ into product molecules $P$ is given by the Michaelis–Menten rule

(2.2) $$\begin{eqnarray}\unicode[STIX]{x1D705}(s)=\frac{\unicode[STIX]{x1D705}_{1}\unicode[STIX]{x1D705}_{2}s}{\unicode[STIX]{x1D705}_{2}+\unicode[STIX]{x1D705}_{1}s},\end{eqnarray}$$

where $s$ is the substrate concentration, $\unicode[STIX]{x1D705}_{1}=k_{1}k_{2}/(k_{-1}+k_{2})$ and $\unicode[STIX]{x1D705}_{2}=k_{2}$ . Note that $\unicode[STIX]{x1D705}(s)$ is the rate of reaction per enzyme molecule, so the total flux at the surface of the swimmer will be the reaction rate $\unicode[STIX]{x1D705}(s)$ multiplied by a measure of surface activity $\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})$ , which quantifies the abundance of enzyme molecules at a particular point on the surface of the particle, parameterised by the spherical coordinate angles $(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})$ in a body-fixed frame of reference.

Figure 1. Sketch of the problem. (a) Axisymmetric active particle with axis of symmetry along $z$ . The surface properties (surface activity $\unicode[STIX]{x1D70E}$ , substrate mobility $\unicode[STIX]{x1D707}_{s}$ and product mobility $\unicode[STIX]{x1D707}_{p}$ ) are functions of the polar angle $\unicode[STIX]{x1D703}$ alone. (b) Active particle placed in a uniform substrate gradient, $\unicode[STIX]{x1D735}s^{\infty }$ , which lies in the $yz$ -plane. The presence of the background gradient, together with the catalytic coating of the swimmer, induce local tangential chemical gradients on the surface of the sphere (light grey arrows) which give rise to an apparent fluid slip velocity (white arrows) leading to translation (linear velocity $\boldsymbol{V}$ ) and rotation (angular velocity $\unicode[STIX]{x1D74E}$ ) of the particle.

We next assume that the spherical swimmer has axisymmetric surface properties, and define the $z$ -axis to be its axis of symmetry (see figure 1 a). We choose the origin at the centre of the sphere, and define our spherical coordinate system such that the polar angle $\unicode[STIX]{x1D703}$ is measured relative to the axis of symmetry of the swimmer. In this coordinate system, any surface property of the swimmer can be expressed as series of Legendre polynomials of the polar angle $P_{l}(\cos \unicode[STIX]{x1D703})$ . In particular, the axisymmetric surface activity $\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D703})$ can be written as

(2.3) $$\begin{eqnarray}\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D703})=\mathop{\sum }_{l=0}^{\infty }\unicode[STIX]{x1D70E}_{l}P_{l}(\cos \unicode[STIX]{x1D703}).\end{eqnarray}$$

We note that the typical scale of substrate concentration surrounding a swimmer centred at position $\boldsymbol{r}_{c}$ is given by the background concentration of substrate at that position, i.e. $s\sim s_{b}(\boldsymbol{r}_{c})$ . It is this scale that determines the magnitude of the reaction rate on the surface of the swimmer. Using this scale, the nonlinear reaction rate from (2.2) has two distinguished limits in which analytical progress can be made.

In the case when $s_{b}(\boldsymbol{r}_{c})\gg \unicode[STIX]{x1D705}_{2}/\unicode[STIX]{x1D705}_{1}$ , the reaction saturates at a constant rate $\unicode[STIX]{x1D705}_{2}$ . This is the limit studied by Golestanian et al. (Reference Golestanian, Liverpool and Ajdari2005, Reference Golestanian, Liverpool and Ajdari2007), Popescu et al. (Reference Popescu, Dietrich, Tasinkevych and Ralston2010) and Michelin & Lauga (Reference Michelin and Lauga2015), amongst others. The swimmer could still perform chemotaxis in this regime but solely due to the substrate, since the product problem has become independent of the substrate dynamics. If the product is released at a constant rate $\unicode[STIX]{x1D705}_{2}\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D703})$ from the surface of the swimmer, it can only lead to translation along the axis of symmetry of the swimmer. Therefore, there is no chemotactic effect from the product in this regime.

In contrast, when $s_{b}(\boldsymbol{r}_{c})\ll \unicode[STIX]{x1D705}_{2}/\unicode[STIX]{x1D705}_{1}$ , the reaction rate becomes proportional to the substrate concentration,

(2.4) $$\begin{eqnarray}\unicode[STIX]{x1D705}(s)\sim \unicode[STIX]{x1D705}_{1}s,\end{eqnarray}$$

which means that the release of product molecules at the surface of the swimmer also depends on the direction of the background gradient in substrate concentration. In this limit, both substrate and product can contribute to chemotaxis in a non-trivial way.

In what follows we shall focus on the linear regime (2.4) which provides a mathematically tractable problem, yet one that leads to interesting chemotactic effects. Traditional diffusiophoresis (whereby a passive colloid is propelled due to a background gradient in the concentration of a chemical species) will also be captured in our model as the special case where the surface activity $\unicode[STIX]{x1D70E}$ is zero, but the presence of chemical reactions will add further complexity to the chemotactic behaviour of the swimmer.

In the limit of zero Péclet number for both substrate and product, we can neglect the advection of molecules and solve a diffusion-reaction problem with the appropriate boundary conditions. In the bulk of the fluid, we must solve the steady diffusion equation for each component

(2.5a,b ) $$\begin{eqnarray}D_{s}\unicode[STIX]{x1D6FB}^{2}s=0,\quad D_{p}\unicode[STIX]{x1D6FB}^{2}p=0,\end{eqnarray}$$

where $s$ and $p$ are the volume concentrations of substrate and product molecules, while $D_{s}$ and $D_{p}$ are their respective diffusivities. The diffusive flux normal to the particle surface is given by the consumption and production of substrate and product molecules, respectively. Given the definition of the surface activity, $\unicode[STIX]{x1D70E}$ , we thus write these boundary conditions as

(2.6a,b ) $$\begin{eqnarray}\left.-D_{s}\frac{\unicode[STIX]{x2202}s}{\unicode[STIX]{x2202}n}\right|_{|\boldsymbol{r}|=R}=\left.-\unicode[STIX]{x1D705}_{1}\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D703})s\right|_{|\boldsymbol{r}|=R},\quad \left.-D_{p}\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}n}\right|_{|\boldsymbol{r}|=R}=\left.\unicode[STIX]{x1D705}_{1}\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D703})s\right|_{|\boldsymbol{r}|=R}.\end{eqnarray}$$

Finally, we impose that the perturbation due to the presence of the particle decays at infinity, such that we recover the background concentrations far away from the particle. Therefore, as $|\boldsymbol{r}|\rightarrow \infty$ , we impose that

(2.7a,b ) $$\begin{eqnarray}s\rightarrow s_{b}(\boldsymbol{r}_{c})+\boldsymbol{r}\boldsymbol{\cdot }\unicode[STIX]{x1D735}s^{\infty },\quad p\rightarrow 0,\end{eqnarray}$$

where $\boldsymbol{r}_{c}$ is the instantaneous position of the centre of the swimmer relative to the reference location where $s_{b}(\boldsymbol{r}^{\prime })=0$ , and $\boldsymbol{r}=\boldsymbol{r}^{\prime }-\boldsymbol{r}_{c}$ is the position relative to the centre of the swimmer.

2.2 Set-up: hydrodynamics

The physical basis for the self-propulsion of phoretic swimmers lies in the interactions between the chemical molecules present in the fluid and the surface of the particle (see sketch in figure 1 b with notation). Following the classical framework (Anderson Reference Anderson1989), interactions are restricted to a thin diffuse layer around the particle where the chemical molecules obey a Boltzmann distribution. The energy of the particle is given by a total interaction potential $\unicode[STIX]{x1D6F7}$ which includes effects such as van der Waals forces and excluded volume effects. We assume that the size of the boundary layer is much smaller than the inverse of the local curvature everywhere on the surface of the particle, so that we can use a locally flat approximation.

If there is a gradient in the concentration of the chemical species outside the diffuse boundary layer, it will give rise generically to an osmotic flow near the surface of the swimmer. This flow will lead to an apparent slip velocity at the edge of the diffuse boundary layer proportional to the tangential gradient in concentration of the chemical species. The constant of proportionality is called the diffusio–osmotic mobility, $\unicode[STIX]{x1D707}$ , and it is a property of the surface interactions. A detailed description of diffusiophoresis can be found in Anderson (Reference Anderson1989), accompanied by the classical integral expression for the mobility $\unicode[STIX]{x1D707}$ in terms of the total interaction potential $\unicode[STIX]{x1D6F7}$ . For the purpose of the present paper, we may ignore the details of the interaction potential and simply work with the mobility as a given surface property.

Because our swimmer is axisymmetric, we can also express the diffusio–osmotic mobility as a series of Legendre polynomials, as we did for the surface activity. In general there are two relevant mobilities, one for the substrate molecules ( $\unicode[STIX]{x1D707}_{s}$ ) and the other for the product molecules ( $\unicode[STIX]{x1D707}_{p}$ ), both of which we decompose as

(2.8a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D707}_{s}(\unicode[STIX]{x1D703})=\mathop{\sum }_{l=0}^{\infty }\unicode[STIX]{x1D707}_{sl}P_{l}(\cos \unicode[STIX]{x1D703}),\quad \unicode[STIX]{x1D707}_{p}(\unicode[STIX]{x1D703})=\mathop{\sum }_{l=0}^{\infty }\unicode[STIX]{x1D707}_{pl}P_{l}(\cos \unicode[STIX]{x1D703}).\end{eqnarray}$$

In the case of our spherical particle, the slip velocity due to either one of the chemical species at a point $(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})$ on the surface of the swimmer can then be written as

(2.9) $$\begin{eqnarray}\boldsymbol{v}_{slip}^{\unicode[STIX]{x1D6FC}}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})=\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D6FC}}(\unicode[STIX]{x1D703})(\boldsymbol{I}-\hat{\boldsymbol{r}}\hat{\boldsymbol{r}})\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}|_{r=R}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719}),\end{eqnarray}$$

where $\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D6FC}}$ is the diffusio–osmotic mobility of the surface with respect to the chemical species in question, $(\boldsymbol{I}-\hat{\boldsymbol{r}}\hat{\boldsymbol{r}})\boldsymbol{\cdot }\unicode[STIX]{x1D735}$ is the tangential gradient operator on the surface of the sphere, and $\unicode[STIX]{x1D6FC}$ is the concentration of the chemical species of interest (substrate or product).

The final step in the set-up of the problem is to translate the local slip velocity into the equations of motion of the particle. Within the low Reynolds number regime, inertia can be neglected and therefore it is relevant to talk about the instantaneous linear and angular velocities of the particle. A classical calculation using the reciprocal theorem on a force-free, torque-free sphere (see Stone & Samuel Reference Stone and Samuel1996) establishes that the linear velocity is given by

(2.10) $$\begin{eqnarray}\boldsymbol{V}^{\unicode[STIX]{x1D6FC}}=-\frac{1}{4\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}}\int _{0}^{\unicode[STIX]{x03C0}}\boldsymbol{v}_{slip}^{\unicode[STIX]{x1D6FC}}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})\sin \unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D719},\end{eqnarray}$$

and the angular velocity by

(2.11) $$\begin{eqnarray}\unicode[STIX]{x1D74E}^{\unicode[STIX]{x1D6FC}}=-\frac{3}{8\unicode[STIX]{x03C0}R}\int _{0}^{2\unicode[STIX]{x03C0}}\int _{0}^{\unicode[STIX]{x03C0}}\hat{\boldsymbol{r}}\wedge \boldsymbol{v}_{slip}^{\unicode[STIX]{x1D6FC}}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})\sin \unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D719}.\end{eqnarray}$$

Using the short-hand $\langle \cdots \rangle$ to denote the surface average, we can write these as

(2.12a,b ) $$\begin{eqnarray}\boldsymbol{V}^{\unicode[STIX]{x1D6FC}}=-\langle \boldsymbol{v}_{slip}^{\unicode[STIX]{x1D6FC}}\rangle ,\quad \unicode[STIX]{x1D74E}^{\unicode[STIX]{x1D6FC}}=-\frac{3}{2R}\langle \hat{\boldsymbol{r}}\wedge \boldsymbol{v}_{slip}^{\unicode[STIX]{x1D6FC}}\rangle .\end{eqnarray}$$

Due to the linearity of Stokes flow and of the chemical problem, the contributions from substrate ( $s$ ) and product ( $p$ ) can be calculated separately and added up at the end to give the total linear and angular velocities as

(2.13a,b ) $$\begin{eqnarray}\boldsymbol{V}=\boldsymbol{V}^{s}+\boldsymbol{V}^{p},\quad \unicode[STIX]{x1D74E}=\unicode[STIX]{x1D74E}^{s}+\unicode[STIX]{x1D74E}^{p}.\end{eqnarray}$$

In the simple case of a particle with uniform surface mobility and no surface activity, the only variations in concentration are due to the linear background gradient. If the interaction potential between the surface of the particle and the chemical species were attractive (respectively repulsive), then this would give rise to a negative (respectively positive) mobility (see Anderson Reference Anderson1989). Thus, the osmotic flows induced near the surface would go against (or up) the tangential chemical gradient. According to (2.10) this means that the particle in this case would move in the direction of increasing (respectively decreasing) chemical concentration, as expected from a colloid which displays attractive (respectively repulsive) interactions with the chemical species. This is the physical basis for the chemotaxis of phoretic swimmers, and the following subsections are concerned with formulating a precise mathematical description of the chemotactic mechanism.

2.3 Small Damköhler number approximation

In its current form, the boundary condition (2.6) at the surface of the sphere does not admit a closed-form solution for the substrate. In order to make further progress, we must make a final simplifying assumption about the relative importance of reaction and diffusion for the substrate molecules. This balance is quantified by the Damköhler number, $Da$ , which we define as

(2.14) $$\begin{eqnarray}Da\equiv \frac{\unicode[STIX]{x1D705}_{1}\unicode[STIX]{x1D6F4}R}{D_{s}}\sim \frac{|\unicode[STIX]{x1D705}_{1}\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D703})s|}{|D_{s}\unicode[STIX]{x2202}s/\unicode[STIX]{x2202}n|},\end{eqnarray}$$

where $\unicode[STIX]{x1D705}_{1}$ is the volumetric rate or reaction (in units $\text{L}^{3}\text{T}^{-1}$ ) and $\unicode[STIX]{x1D6F4}$ is the typical scale of surface activity (in units $\text{L}^{-2}$ ). In the limit of small Damköhler number, $Da\ll 1$ , the particle behaves as a passive sphere at leading order ( $Da^{0}$ ) but, at first order ( $Da^{1}$ ), it produces small disturbances in the substrate and product distributions due to the chemical reactivity of its surface.

To make this statement more rigorous, we proceed with non-dimensionalising the equations of the chemical problem. For the chemical concentrations we use the scale imposed by the given background concentration of substrate, such that our dimensionless variables are $\tilde{s}=s/s_{b}(\boldsymbol{r}_{c})$ and $\tilde{p}=p/s_{b}(\boldsymbol{r}_{c})$ . We scale lengths by the radius of the swimmer, $\tilde{r}=r/R$ , and the surface activity by its typical scale, $\tilde{\unicode[STIX]{x1D70E}}=\unicode[STIX]{x1D70E}/\unicode[STIX]{x1D6F4}$ . Under these scalings, the dimensionless versions of (2.5)–(2.7) are

(2.15a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D6FB}^{2}\tilde{s}=0,\quad \unicode[STIX]{x1D6FB}^{2}\tilde{p}=0,\end{eqnarray}$$

in the bulk of the fluid, with boundary conditions on the surface of the swimmer

(2.16a,b ) $$\begin{eqnarray}\left.\frac{\unicode[STIX]{x2202}\tilde{s}}{\unicode[STIX]{x2202}n}\right|_{\tilde{r}=1}=Da\tilde{\unicode[STIX]{x1D70E}}(\unicode[STIX]{x1D703})\tilde{s}|_{\tilde{r}=1},\quad \left.\frac{\unicode[STIX]{x2202}\tilde{p}}{\unicode[STIX]{x2202}n}\right|_{\tilde{r}=1}=-\unicode[STIX]{x1D6FF}Da\tilde{\unicode[STIX]{x1D70E}}(\unicode[STIX]{x1D703})\tilde{s}|_{\tilde{r}=1},\end{eqnarray}$$

and away from the swimmer

(2.17a,b ) $$\begin{eqnarray}\tilde{s}\rightarrow 1+\tilde{\boldsymbol{r}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\tilde{s}^{\infty },\quad \tilde{p}\rightarrow 0\quad \text{as}~\tilde{r}\rightarrow \infty .\end{eqnarray}$$

In addition to the Damköhler number, a second dimensionless parameter appears, $\unicode[STIX]{x1D6FF}=D_{s}/D_{p}$ , which is the ratio of the two diffusivities.

Due to the existence of a background substrate profile, the non-dimensional substrate distribution around the swimmer will include an $O(1)$ term, such that $\tilde{s}=\tilde{s}_{0}+Da\tilde{s}_{1}+O(Da^{2})$ . This can be seen from the boundary conditions (2.16) and (2.17). However, since there was no product in the ambient medium to begin with, the scale of the product distribution is determined by the chemical reactions at the surface of the swimmer. In the small Damköhler number approximation the swimmer is only weakly reactive, which means that the leading-order product distribution appears only at $O(Da)$ , such that $\tilde{p}=Da\tilde{p}_{1}+O(Da^{2})$ .

Due to the linearity of the Stokes equations, the contributions from substrate and product add up to give the total linear and angular velocities of the swimmer. Therefore, if we want to include all terms of equal importance for the motion of the swimmer, we must calculate the substrate and product concentrations to the same level of accuracy in our Damköhler number expansion. In the current paper we work out the solution up to $O(Da)$ , which brings the first non-trivial contributions to the final answer and also captures the expected phenomenology.

Returning to dimensional variables, we will approximate the substrate concentration to be $s\simeq s_{0}+s_{1}$ , where $s_{0}$ is the solution of the leading-order problem

(2.18a-c ) $$\begin{eqnarray}\unicode[STIX]{x1D6FB}^{2}s_{0}=0,\quad \left.\frac{\unicode[STIX]{x2202}s_{0}}{\unicode[STIX]{x2202}n}\right|_{r=R}=0,\quad s_{0}\rightarrow s_{b}(\boldsymbol{r}_{c})+\boldsymbol{r}\boldsymbol{\cdot }\unicode[STIX]{x1D735}s^{\infty }~\text{as}~r\rightarrow \infty ,\end{eqnarray}$$

while $s_{1}$ is the solution of the first-order problem

(2.19a-c ) $$\begin{eqnarray}\unicode[STIX]{x1D6FB}^{2}s_{1}=0,\quad \left.\frac{\unicode[STIX]{x2202}s_{1}}{\unicode[STIX]{x2202}n}\right|_{r=R}=\frac{\unicode[STIX]{x1D705}_{1}\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D703})}{D_{s}}\left.s_{0}\right|_{r=R},\quad s_{1}\rightarrow 0~\text{as}~r\rightarrow \infty .\end{eqnarray}$$

Similarly, we approximate the product concentration to be $p\simeq p_{1}$ , where $p_{1}$ is the solution of the first-order problem

(2.20a-c ) $$\begin{eqnarray}\unicode[STIX]{x1D6FB}^{2}p_{1}=0,\quad \left.\frac{\unicode[STIX]{x2202}p_{1}}{\unicode[STIX]{x2202}n}\right|_{r=R}=-\frac{\unicode[STIX]{x1D705}_{1}\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D703})}{D_{p}}\left.s_{0}\right|_{r=R},\quad p_{1}\rightarrow 0~\text{as}~r\rightarrow \infty .\end{eqnarray}$$

Notice that the first-order problems for the substrate and product are the same, up to a factor of $-D_{s}/D_{p}$ , so we only need to solve for one of them.

2.4 Leading-order effects

Recall that our $z$ -axis is the axis of symmetry of the swimmer. Without loss of generality, we may then choose our $y$ -axis such that the substrate gradient lies in the $yz$ -plane, that is $\unicode[STIX]{x1D735}s^{\infty }=\unicode[STIX]{x1D735}s_{y}^{\infty }\hat{\boldsymbol{y}}+\unicode[STIX]{x1D735}s_{z}^{\infty }\hat{\boldsymbol{z}}$ .

The leading-order substrate concentration can be written as the sum of the background concentration plus a perturbation which satisfies Laplace’s equation and decays at infinity. The general solution to Laplace’s equation can be written classically as a sum of spherical harmonics $Y_{l}^{m}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})\propto \text{e}^{\text{i}m\unicode[STIX]{x1D719}}P_{l}^{m}(\cos \unicode[STIX]{x1D703})$ , where $P_{l}^{m}$ is an associated Legendre polynomial, multiplied by an appropriate radial component. For our perturbation to decay at infinity this radial component has to be a multiple of $r^{-(l+1)}$ , so we can write

(2.21) $$\begin{eqnarray}s_{0}(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})=s_{b}(\boldsymbol{r}_{c})+\boldsymbol{r}\boldsymbol{\cdot }\unicode[STIX]{x1D735}s^{\infty }+\mathop{\sum }_{l=0}^{\infty }\mathop{\sum }_{m=-l}^{l}f_{l}^{m}r^{-(l+1)}Y_{l}^{m}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719}).\end{eqnarray}$$

This expression can be simplified further if we take into account that the zero-flux boundary condition from (2.18b ) can only excite the spherical harmonics of order $m=0$ and $m=\pm 1$ , which are already present in the solution due to the form of the concentration at infinity. Due to our choice of $y$ -axis, the $m=\pm 1$ components must only lead to dependence on $\sin \unicode[STIX]{x1D719}$ and not $\cos \unicode[STIX]{x1D719}$ , so the substrate concentration has at most the form

(2.22) $$\begin{eqnarray}s_{0}(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})=s_{b}(\boldsymbol{r}_{c})+\boldsymbol{r}\boldsymbol{\cdot }\unicode[STIX]{x1D735}s^{\infty }+\mathop{\sum }_{l=0}^{\infty }\left(\frac{R}{r}\right)^{l+1}[a_{l}P_{l}(\cos \unicode[STIX]{x1D703})+b_{l}P_{l}^{1}(\cos \unicode[STIX]{x1D703})\sin \unicode[STIX]{x1D719}].\end{eqnarray}$$

Imposing the zero-flux boundary condition from (2.18b ), we find the coefficients

(2.23a,b ) $$\begin{eqnarray}a_{1}=\frac{\unicode[STIX]{x1D735}s_{z}^{\infty }R}{2},\quad b_{1}=-\frac{\unicode[STIX]{x1D735}s_{y}^{\infty }R}{2},\end{eqnarray}$$

and $a_{l},b_{l}=0$ for $l\neq 1$ , which lead to the following substrate distribution

(2.24) $$\begin{eqnarray}s_{0}=s_{b}(\boldsymbol{r}_{c})+\boldsymbol{r}\boldsymbol{\cdot }\unicode[STIX]{x1D735}s^{\infty }+\frac{R^{3}}{2r^{2}}(\unicode[STIX]{x1D735}s_{z}^{\infty }\cos \unicode[STIX]{x1D703}+\unicode[STIX]{x1D735}s_{y}^{\infty }\sin \unicode[STIX]{x1D703}\sin \unicode[STIX]{x1D719}).\end{eqnarray}$$

Notably, this expression can be simplified into

(2.25) $$\begin{eqnarray}s_{0}=s_{b}(\boldsymbol{r}_{c})+\boldsymbol{r}\boldsymbol{\cdot }\unicode[STIX]{x1D735}s^{\infty }\left(1+\frac{1}{2}\left(\frac{R}{r}\right)^{3}\right).\end{eqnarray}$$

In order to compute the slip velocity, $\boldsymbol{v}_{slip}^{s,0}$ , due to the leading-order substrate profile, we have to evaluate

(2.26) $$\begin{eqnarray}(\boldsymbol{I}-\hat{\boldsymbol{r}}\hat{\boldsymbol{r}})\boldsymbol{\cdot }\unicode[STIX]{x1D735}s_{0}=\frac{1}{r}\frac{\unicode[STIX]{x2202}s_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}\hat{\unicode[STIX]{x1D73D}}+\frac{1}{r\sin \unicode[STIX]{x1D703}}\frac{\unicode[STIX]{x2202}s_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}\hat{\unicode[STIX]{x1D753}}\end{eqnarray}$$

on the surface of the swimmer. Using (2.9), we obtain

(2.27) $$\begin{eqnarray}\boldsymbol{v}_{slip}^{s,0}=\frac{3\unicode[STIX]{x1D707}_{s}(\unicode[STIX]{x1D703})}{2}[(-\unicode[STIX]{x1D735}s_{z}^{\infty }\sin \unicode[STIX]{x1D703}+\unicode[STIX]{x1D735}s_{y}^{\infty }\cos \unicode[STIX]{x1D703}\sin \unicode[STIX]{x1D719})\hat{\unicode[STIX]{x1D73D}}+\unicode[STIX]{x1D735}s_{y}^{\infty }\cos \unicode[STIX]{x1D719}\,\hat{\unicode[STIX]{x1D753}}].\end{eqnarray}$$

Then, by averaging this slip velocity over the surface of the spherical swimmer, we determine the contribution to linear velocity, $\boldsymbol{V}_{0}^{s}$ , coming from the leading-order substrate problem as

(2.28) $$\begin{eqnarray}\boldsymbol{V}_{0}^{s}=-\left(\unicode[STIX]{x1D707}_{s0}+{\textstyle \frac{1}{10}}\unicode[STIX]{x1D707}_{s2}\right)\unicode[STIX]{x1D735}s^{\infty }+{\textstyle \frac{3}{10}}\unicode[STIX]{x1D707}_{s2}\,\hat{\boldsymbol{z}}\hat{\boldsymbol{z}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}s^{\infty }.\end{eqnarray}$$

The calculations leading to this result are detailed in the first subsection of appendix A.

Although this result is valid for a general surface mobility, $\unicode[STIX]{x1D707}_{s}(\unicode[STIX]{x1D703})$ , only two Legendre modes of the mobility contribute to the final results, namely those adjacent to the forcing coming from the linear gradient. Furthermore, our result is in agreement with the linear velocity calculated by Anderson for a passive sphere placed in a non-uniform field (equation (37a) from Anderson Reference Anderson1989). In our case that field is the concentration of substrate molecules, and the sphere is passive, at leading order, due to our small Damköhler number approximation.

Similarly, one has to compute the average of $\hat{\boldsymbol{r}}\wedge \boldsymbol{v}_{slip}^{s,0}$ over the surface of the spherical swimmer in order to determine the angular velocity contribution, $\unicode[STIX]{x1D74E}_{0}^{s}$ , resulting from the leading-order substrate problem

(2.29) $$\begin{eqnarray}\unicode[STIX]{x1D74E}_{0}^{s}=-\frac{3\unicode[STIX]{x1D707}_{s1}}{4R}\,\hat{\boldsymbol{z}}\wedge \unicode[STIX]{x1D735}s^{\infty }.\end{eqnarray}$$

The calculations leading to this are detailed in the second subsection of appendix A.

This result agrees with the angular velocity calculated by Anderson for a passive sphere placed in a non-uniform field (equation (37b) from Anderson Reference Anderson1989). Note again that this result is valid for a general surface mobility $\unicode[STIX]{x1D707}_{s}(\unicode[STIX]{x1D703})$ , but only one mode contributes to the reorientation of a passive swimmer. This mode corresponds to a linear dependence of mobility on the cosine of the polar angle, with opposite and equal in magnitude mobilities at the two poles. By symmetry we would expect this type of swimmer to align with the chemical gradient, but we can also predict the direction of rotation using a simple physical argument. Following the negative and positive chemotactic response of swimmers with uniform positive and negative mobility, respectively, as discussed at the end of § 2.2, we would expect the pole of negative mobility to pull towards the higher concentration of substrate and the pole of positive mobility to pull towards the lower concentration. This will then make the swimmer rotate until its axis of symmetry is aligned with the substrate gradient and the pole of negative mobility faces higher concentration, which is precisely the dynamics captured by (2.29). Since the two poles are equally strong, the tug-of-war cannot be won by either opponent and this mode has zero linear velocity, in agreement with (2.28).

2.5 First-order effects

We now turn our attention to finding the corrections to first order in Damköhler number for the substrate problem. In order to satisfy Laplace’s equation and the decay at infinity, the first-order substrate concentration must have the general form

(2.30) $$\begin{eqnarray}s_{1}(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})=\mathop{\sum }_{l=0}^{\infty }\mathop{\sum }_{m=-l}^{l}F_{l}^{m}r^{-(l+1)}Y_{l}^{m}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719}),\end{eqnarray}$$

but the normal-flux boundary condition from (2.19b ) can only excite modes with either no $\unicode[STIX]{x1D719}$ -dependence or $\sin \unicode[STIX]{x1D719}$ -dependence which appear in the leading-order substrate concentration from (2.25). This means that the first-order substrate concentration can be simplified at most to the sum

(2.31) $$\begin{eqnarray}s_{1}(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})=\mathop{\sum }_{l=0}^{\infty }\left(\frac{R}{r}\right)^{l+1}[A_{l}P_{l}(\cos \unicode[STIX]{x1D703})-B_{l}\sin \unicode[STIX]{x1D703}P_{l}^{\prime }(\cos \unicode[STIX]{x1D703})\sin \unicode[STIX]{x1D719}],\end{eqnarray}$$

where we have used the fact that $P_{l}^{1}(\cos \unicode[STIX]{x1D703})=-\sin \unicode[STIX]{x1D703}P_{l}^{\prime }(\cos \unicode[STIX]{x1D703})$ . The coefficients $A_{l}$ and $B_{l}$ can be determined from the normal-flux boundary condition in (2.19). These calculations are detailed in the third subsection of appendix A, and lead to the following expressions

(2.32) $$\begin{eqnarray}\displaystyle & \displaystyle A_{l}=-\frac{\unicode[STIX]{x1D705}_{1}R}{(l+1)D_{s}}\left[s_{b}(\boldsymbol{r}_{c})\unicode[STIX]{x1D70E}_{l}+\frac{3}{2}\unicode[STIX]{x1D735}s_{z}^{\infty }R\left(\frac{l+1}{2l+3}\unicode[STIX]{x1D70E}_{l+1}+\frac{l}{2l-1}\unicode[STIX]{x1D70E}_{l-1}\right)\right], & \displaystyle\end{eqnarray}$$
(2.33) $$\begin{eqnarray}\displaystyle & \displaystyle B_{l}=-\frac{3\unicode[STIX]{x1D705}_{1}R^{2}\unicode[STIX]{x1D735}s_{y}^{\infty }}{2(l+1)D_{s}}\left(\frac{\unicode[STIX]{x1D70E}_{l+1}}{2l+3}-\frac{\unicode[STIX]{x1D70E}_{l-1}}{2l-1}\right). & \displaystyle\end{eqnarray}$$

Since the first-order product problem (2.20) is the same as the first-order substrate problem (2.19), except for a factor of $-D_{s}/D_{p}$ , we can use our results for $s_{1}$ and state that the first-order product concentration $p_{1}$ is given by

(2.34) $$\begin{eqnarray}p_{1}=\mathop{\sum }_{l=0}^{\infty }\left(\frac{R}{r}\right)^{l+1}[C_{l}P_{l}(\cos \unicode[STIX]{x1D703})-D_{l}\sin \unicode[STIX]{x1D703}P_{l}^{\prime }(\cos \unicode[STIX]{x1D703})\sin \unicode[STIX]{x1D719}],\end{eqnarray}$$

where the coefficients $C_{l}$ and $D_{l}$ can be written as

(2.35) $$\begin{eqnarray}\displaystyle & \displaystyle C_{l}=\frac{\unicode[STIX]{x1D705}_{1}R}{(l+1)D_{p}}\left[s_{b}(\boldsymbol{r}_{c})\unicode[STIX]{x1D70E}_{l}+\frac{3}{2}\unicode[STIX]{x1D735}s_{z}^{\infty }R\left(\frac{l+1}{2l+3}\unicode[STIX]{x1D70E}_{l+1}+\frac{l}{2l-1}\unicode[STIX]{x1D70E}_{l-1}\right)\right], & \displaystyle\end{eqnarray}$$
(2.36) $$\begin{eqnarray}\displaystyle & \displaystyle D_{l}=\frac{3\unicode[STIX]{x1D705}_{1}R^{2}\unicode[STIX]{x1D735}s_{y}^{\infty }}{2(l+1)D_{p}}\left(\frac{\unicode[STIX]{x1D70E}_{l+1}}{2l+3}-\frac{\unicode[STIX]{x1D70E}_{l-1}}{2l-1}\right). & \displaystyle\end{eqnarray}$$

We next need to calculate the contribution of the first-order substrate concentration to the slip velocity. Using (2.9) and (2.31), we find that

(2.37) $$\begin{eqnarray}\displaystyle \boldsymbol{v}_{slip}^{s,1} & = & \displaystyle \frac{\unicode[STIX]{x1D707}_{s}(\unicode[STIX]{x1D703})}{R}\mathop{\sum }_{l=0}^{\infty }\left(-A_{l}\sin \unicode[STIX]{x1D703}P_{l}^{\prime }(\cos \unicode[STIX]{x1D703})-B_{l}\cos \unicode[STIX]{x1D703}P_{l}^{\prime }(\cos \unicode[STIX]{x1D703})\sin \unicode[STIX]{x1D719}\right.\nonumber\\ \displaystyle & & \displaystyle +\,\left.B_{l}\sin ^{2}\unicode[STIX]{x1D703}P_{l}^{\prime \prime }(\cos \unicode[STIX]{x1D703})\sin \unicode[STIX]{x1D719}\right)\hat{\unicode[STIX]{x1D73D}}+\frac{\unicode[STIX]{x1D707}_{s}(\unicode[STIX]{x1D703})}{R}\mathop{\sum }_{l=0}^{\infty }(-B_{l}P_{l}^{\prime }(\cos \unicode[STIX]{x1D703})\cos \unicode[STIX]{x1D719})\hat{\unicode[STIX]{x1D753}}.\end{eqnarray}$$

The last two subsections in appendix A present the detailed calculations involved in averaging $\boldsymbol{v}_{slip}^{s,1}$ and $\hat{\boldsymbol{r}}\wedge \boldsymbol{v}_{slip}^{s,1}$ over the surface of the spherical swimmer, in order to calculate the contributions of the first-order substrate problem to the linear and angular velocities of the swimmer. The final results are that the first-order substrate problem contributes

(2.38) $$\begin{eqnarray}\displaystyle \boldsymbol{V}_{1}^{s} & = & \displaystyle -\frac{\unicode[STIX]{x1D705}_{1}s_{b}(\boldsymbol{r}_{c})}{D_{s}}\mathop{\sum }_{l=1}^{\infty }\left(\frac{l}{2l+1}\right)\unicode[STIX]{x1D70E}_{l}\left(\frac{\unicode[STIX]{x1D707}_{s,l+1}}{2l+3}-\frac{\unicode[STIX]{x1D707}_{s,l-1}}{2l-1}\right)\hat{\boldsymbol{z}}\nonumber\\ \displaystyle & & \displaystyle -\,\frac{3\unicode[STIX]{x1D705}_{1}R}{4D_{s}}\mathop{\sum }_{l=1}^{\infty }\left(\frac{l}{2l+1}\right)\left(\frac{\unicode[STIX]{x1D70E}_{l+1}}{2l+3}-\frac{\unicode[STIX]{x1D70E}_{l-1}}{2l-1}\right)\left(\frac{l+1}{2l-1}\unicode[STIX]{x1D707}_{s,l-1}+\frac{l}{2l+3}\unicode[STIX]{x1D707}_{s,l+1}\right)\unicode[STIX]{x1D735}s^{\infty }\nonumber\\ \displaystyle & & \displaystyle +\,\frac{3\unicode[STIX]{x1D705}_{1}R}{4D_{s}}\mathop{\sum }_{l=1}^{\infty }\left(\frac{l}{2l+1}\right)\left(\frac{3(l+1)\unicode[STIX]{x1D70E}_{l+1}\unicode[STIX]{x1D707}_{s,l-1}}{(2l+3)(2l-1)}-\frac{(l+2)\unicode[STIX]{x1D70E}_{l+1}\unicode[STIX]{x1D707}_{s,l+1}}{(2l+3)^{2}}\right.\nonumber\\ \displaystyle & & \displaystyle +\,\left.\frac{(l-1)\unicode[STIX]{x1D70E}_{l-1}\unicode[STIX]{x1D707}_{s,l-1}}{(2l-1)^{2}}-\frac{3l\unicode[STIX]{x1D70E}_{l-1}\unicode[STIX]{x1D707}_{s,l+1}}{(2l-1)(2l+3)}\right)\hat{\boldsymbol{z}}\hat{\boldsymbol{z}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}s^{\infty }\end{eqnarray}$$

to the linear velocity of the swimmer, and

(2.39) $$\begin{eqnarray}\unicode[STIX]{x1D74E}_{1}^{s}=-\frac{9\unicode[STIX]{x1D705}_{1}}{8D_{s}}\mathop{\sum }_{l=1}^{\infty }\left(\frac{l}{2l+1}\right)\unicode[STIX]{x1D707}_{sl}\left(\frac{\unicode[STIX]{x1D70E}_{l+1}}{2l+3}-\frac{\unicode[STIX]{x1D70E}_{l-1}}{2l-1}\right)\,\hat{\boldsymbol{z}}\wedge \unicode[STIX]{x1D735}s^{\infty }\end{eqnarray}$$

to its angular velocity.

The contributions from the first-order product problem will be the same as the ones above, with an additional factor of $-D_{s}/D_{p}$ which carries through from the solution of the chemical problem, and a change of subscript on the mobility coefficients $\unicode[STIX]{x1D707}_{sl}\mapsto \unicode[STIX]{x1D707}_{pl}$ .

2.6 Summary of results

In summary, we may write the instantaneous linear and angular velocities of the phoretic swimmer, when its centre is at position $\boldsymbol{r}_{c}$ , in the compact form

(2.40) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{V}=U(\boldsymbol{r}_{c})\hat{\boldsymbol{z}}+\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D735}s^{\infty }+\unicode[STIX]{x1D6FD}\,\hat{\boldsymbol{z}}\hat{\boldsymbol{z}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}s^{\infty }, & \displaystyle\end{eqnarray}$$
(2.41) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D74E}=\unicode[STIX]{x1D6F7}\,\hat{\boldsymbol{z}}\wedge \unicode[STIX]{x1D735}s^{\infty }. & \displaystyle\end{eqnarray}$$

Due to the linearity of the Stokes equations, the total linear and angular velocities of the phoretic swimmer are simply the sum of contributions from the substrate and the product, at leading and first order. Using the results from the previous subsections, we deduce that the coefficients appearing in (2.40) and (2.41) are given by

(2.42) $$\begin{eqnarray}\displaystyle & \hspace{-160.00024pt}\displaystyle U(\boldsymbol{r}_{c})=\unicode[STIX]{x1D705}_{1}s_{b}(\boldsymbol{r}_{c})\mathop{\sum }_{l=1}^{\infty }\left(\frac{l}{2l+1}\right)\qquad \qquad & \displaystyle \nonumber\\ \displaystyle & \displaystyle \times \,\unicode[STIX]{x1D70E}_{l}\left[\frac{1}{2l+3}\left(\frac{\unicode[STIX]{x1D707}_{p,l+1}}{D_{p}}-\frac{\unicode[STIX]{x1D707}_{s,l+1}}{D_{s}}\right)-\frac{1}{2l-1}\left(\frac{\unicode[STIX]{x1D707}_{p,l-1}}{D_{p}}-\frac{\unicode[STIX]{x1D707}_{s,l-1}}{D_{s}}\right)\right], & \displaystyle\end{eqnarray}$$
(2.43) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FC}=-\left(\unicode[STIX]{x1D707}_{s0}+\frac{1}{10}\unicode[STIX]{x1D707}_{s2}\right)+\frac{3\unicode[STIX]{x1D705}_{1}R}{4}\mathop{\sum }_{l=1}^{\infty }\left(\frac{l}{2l+1}\right)\left(\frac{\unicode[STIX]{x1D70E}_{l+1}}{2l+3}-\frac{\unicode[STIX]{x1D70E}_{l-1}}{2l-1}\right)\qquad & \displaystyle \nonumber\\ \displaystyle & \displaystyle \times \,\left[\frac{l+1}{2l-1}\left(\frac{\unicode[STIX]{x1D707}_{p,l-1}}{D_{p}}-\frac{\unicode[STIX]{x1D707}_{s,l-1}}{D_{s}}\right)\right.+\left.\frac{l}{2l+3}\left(\frac{\unicode[STIX]{x1D707}_{p,l+1}}{D_{p}}-\frac{\unicode[STIX]{x1D707}_{s,l+1}}{D_{s}}\right)\right], & \displaystyle\end{eqnarray}$$
(2.44) $$\begin{eqnarray}\displaystyle & \hspace{-170.00026pt}\displaystyle \unicode[STIX]{x1D6FD}=\frac{3\unicode[STIX]{x1D707}_{s2}}{10}-\frac{3\unicode[STIX]{x1D705}_{1}R}{4}\mathop{\sum }_{l=1}^{\infty }\left(\frac{l}{2l+1}\right) & \displaystyle \nonumber\\ \displaystyle & \hspace{-30.00005pt}\displaystyle \times \,\left[\left(\frac{3(l+1)\unicode[STIX]{x1D70E}_{l+1}}{(2l+3)(2l-1)}+\frac{(l-1)\unicode[STIX]{x1D70E}_{l-1}}{(2l-1)^{2}}\right)\left(\frac{\unicode[STIX]{x1D707}_{p,l-1}}{D_{p}}-\frac{\unicode[STIX]{x1D707}_{s,l-1}}{D_{s}}\right)\right. & \displaystyle \nonumber\\ \displaystyle & \hspace{-20.00003pt}\displaystyle \left.-\,\left(\frac{(l+2)\unicode[STIX]{x1D70E}_{l+1}}{(2l+3)^{2}}+\frac{3l\unicode[STIX]{x1D70E}_{l-1}}{(2l-1)(2l+3)}\right)\left(\frac{\unicode[STIX]{x1D707}_{p,l+1}}{D_{p}}-\frac{\unicode[STIX]{x1D707}_{s,l+1}}{D_{s}}\right)\right], & \displaystyle\end{eqnarray}$$
(2.45) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6F7}=-\frac{3\unicode[STIX]{x1D707}_{s1}}{4R}+\frac{9\unicode[STIX]{x1D705}_{1}}{8}\mathop{\sum }_{l=1}^{\infty }\left(\frac{l}{2l+1}\right)\left(\frac{\unicode[STIX]{x1D707}_{pl}}{D_{p}}-\frac{\unicode[STIX]{x1D707}_{sl}}{D_{s}}\right)\left(\frac{\unicode[STIX]{x1D70E}_{l+1}}{2l+3}-\frac{\unicode[STIX]{x1D70E}_{l-1}}{2l-1}\right). & \displaystyle\end{eqnarray}$$

To the best of our knowledge, the set-up and modelling assumptions used in deriving these results follow those used by Saha et al. (Reference Saha, Golestanian and Ramaswamy2014). Our final equations, (2.42)–(2.45), represent an extension and correction of the instantaneous laws of motion published in Saha et al. (Reference Saha, Golestanian and Ramaswamy2014) to a phoretic swimmer with the most general axisymmetric surface activity $\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D703})$ and surface mobilities $\unicode[STIX]{x1D707}_{s}(\unicode[STIX]{x1D703})$ , $\unicode[STIX]{x1D707}_{p}(\unicode[STIX]{x1D703})$ . It is also important to note that our solution includes the first-order contributions from the substrate, an effect which had been neglected by Saha et al. (Reference Saha, Golestanian and Ramaswamy2014).

2.7 Janus spherical swimmer

Because our extended laws of motion allow for arbitrary axisymmetric surface properties, we can apply our results to one of the most common designs of phoretic swimmers, namely that of a Janus sphere. This type of swimmer, depicted in figure 2, has uniform surface properties on each of its two hemispheres, such that the surface activity is given by

(2.46) $$\begin{eqnarray}\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D703})=\left\{\begin{array}{@{}ll@{}}\unicode[STIX]{x1D70E}^{U},\quad & 0\leqslant \unicode[STIX]{x1D703}<\unicode[STIX]{x03C0}/2,\\[2.0pt] \unicode[STIX]{x1D70E}^{L},\quad & \unicode[STIX]{x03C0}/2<\unicode[STIX]{x1D703}\leqslant \unicode[STIX]{x03C0},\end{array}\right.\end{eqnarray}$$

where the superscript denotes the lower ( $L$ ) and upper ( $U$ ) halves of the swimmer. We also assume similar spatial distributions of the two mobilities $\unicode[STIX]{x1D707}_{s}(\unicode[STIX]{x1D703})$ and $\unicode[STIX]{x1D707}_{p}(\unicode[STIX]{x1D703})$ (see full notation in figure 2).

Figure 2. Schematic representation of a Janus spherical swimmer with surface properties $\{\unicode[STIX]{x1D70E}^{U},\unicode[STIX]{x1D707}_{s}^{U},\unicode[STIX]{x1D707}_{p}^{U}\}$ on the upper hemisphere and $\{\unicode[STIX]{x1D70E}^{L},\unicode[STIX]{x1D707}_{s}^{L},\unicode[STIX]{x1D707}_{p}^{L}\}$ on the lower hemisphere. The linear and angular velocities of this phoretic swimmer are given by (2.50)–(2.53).

With the help of the integral identity (B 12), we can project (2.46) onto the space of Legendre polynomials, and we find that

(2.47) $$\begin{eqnarray}\unicode[STIX]{x1D70E}_{0}={\textstyle \frac{1}{2}}(\unicode[STIX]{x1D70E}^{U}+\unicode[STIX]{x1D70E}^{L}),\end{eqnarray}$$

but that all other even modes vanish

(2.48) $$\begin{eqnarray}\unicode[STIX]{x1D70E}_{2k}=0,\quad k\geqslant 1.\end{eqnarray}$$

For the odd modes, we have

(2.49) $$\begin{eqnarray}\unicode[STIX]{x1D70E}_{2k+1}=(-1)^{k}(\unicode[STIX]{x1D70E}^{U}-\unicode[STIX]{x1D70E}^{L})\frac{(4k+3)(2k+1)!!}{(4k+2)(2k+2)!!},\end{eqnarray}$$

and equivalent identities for the mobility coefficients $\unicode[STIX]{x1D707}_{sl}$ and $\unicode[STIX]{x1D707}_{pl}$ .

Upon substituting these expressions into our results for the linear and angular velocity of a general axisymmetric swimmer, (2.42)–(2.45), we find that the parameters characterising the motion of a Janus sphere in a uniform chemical gradient are given by

(2.50) $$\begin{eqnarray}\displaystyle & \displaystyle U_{Janus}=\frac{\unicode[STIX]{x1D705}_{1}s_{b}}{8}(\unicode[STIX]{x1D70E}^{L}-\unicode[STIX]{x1D70E}^{U})\left(\frac{\unicode[STIX]{x1D707}_{p}^{L}+\unicode[STIX]{x1D707}_{p}^{U}}{D_{p}}-\frac{\unicode[STIX]{x1D707}_{s}^{L}+\unicode[STIX]{x1D707}_{s}^{U}}{D_{s}}\right), & \displaystyle\end{eqnarray}$$
(2.51) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FC}_{Janus}=-\frac{1}{2}(\unicode[STIX]{x1D707}_{s}^{U}+\unicode[STIX]{x1D707}_{s}^{L})-\frac{\unicode[STIX]{x1D705}_{1}R}{8}(\unicode[STIX]{x1D70E}^{U}+\unicode[STIX]{x1D70E}^{L})\left(\frac{\unicode[STIX]{x1D707}_{p}^{L}+\unicode[STIX]{x1D707}_{p}^{U}}{D_{p}}-\frac{\unicode[STIX]{x1D707}_{s}^{L}+\unicode[STIX]{x1D707}_{s}^{U}}{D_{s}}\right) & \displaystyle \nonumber\\ \displaystyle & \displaystyle -\,\frac{3\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D705}_{1}R}{4}(\unicode[STIX]{x1D70E}^{U}-\unicode[STIX]{x1D70E}^{L})\left(\frac{\unicode[STIX]{x1D707}_{p}^{U}-\unicode[STIX]{x1D707}_{p}^{L}}{D_{p}}-\frac{\unicode[STIX]{x1D707}_{s}^{U}-\unicode[STIX]{x1D707}_{s}^{L}}{D_{s}}\right),\ & \displaystyle\end{eqnarray}$$
(2.52) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FD}_{Janus}=-\frac{3\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D705}_{1}R}{4}(\unicode[STIX]{x1D70E}^{U}-\unicode[STIX]{x1D70E}^{L})\left(\frac{\unicode[STIX]{x1D707}_{p}^{U}-\unicode[STIX]{x1D707}_{p}^{L}}{D_{p}}-\frac{\unicode[STIX]{x1D707}_{s}^{U}-\unicode[STIX]{x1D707}_{s}^{L}}{D_{s}}\right), & \displaystyle\end{eqnarray}$$
(2.53) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6F7}_{Janus}=\frac{9}{16R}(\unicode[STIX]{x1D707}_{s}^{L}-\unicode[STIX]{x1D707}_{s}^{U})+\frac{9\unicode[STIX]{x1D705}_{1}}{64}(\unicode[STIX]{x1D70E}^{L}+\unicode[STIX]{x1D70E}^{U})\left(\frac{\unicode[STIX]{x1D707}_{p}^{L}-\unicode[STIX]{x1D707}_{p}^{U}}{D_{p}}-\frac{\unicode[STIX]{x1D707}_{s}^{L}-\unicode[STIX]{x1D707}_{s}^{U}}{D_{s}}\right), & \displaystyle\end{eqnarray}$$

with swimming kinematics as in (2.40) and (2.41). Note that the coefficients $\unicode[STIX]{x1D6FC}_{Janus}$ and $\unicode[STIX]{x1D6FD}_{Janus}$ are more easily calculated from (A 38) than from (2.43)–(2.44), and that they involve a numerical factor

(2.54) $$\begin{eqnarray}\unicode[STIX]{x1D6FE}=\mathop{\sum }_{k=1}^{\infty }k\left(\frac{(2k-1)!!}{(2k-1)(2k+2)!!}\right)^{2}\approx 0.0336.\end{eqnarray}$$

Our result for $U_{Janus}$ agrees with Golestanian et al. (Reference Golestanian, Liverpool and Ajdari2007) if we equate their concentration-independent surface activity $\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D703})$ to the component of our surface activity which is independent of the background gradient, i.e.  $\unicode[STIX]{x1D705}_{1}s_{b}\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D703})$ . Instead of a single chemical mobility to diffusivity ratio $\unicode[STIX]{x1D707}/D$ , we also have an effective value, $(\unicode[STIX]{x1D707}/D)_{eff}=\unicode[STIX]{x1D707}_{p}/D_{p}-\unicode[STIX]{x1D707}_{s}/D_{s}$ , due to the interplay of the two chemical species in our problem.

The simple expressions in (2.50)–(2.53) provide an estimate for the linear and angular velocity of the swimmer based on information about its surface properties alone, and could thus be used in the design and fabrication of phoretic Janus particles.

Note that, at leading order, the only non-zero parameters are $\unicode[STIX]{x1D6FC}_{Janus}$ and $\unicode[STIX]{x1D6F7}_{Janus}$ . Thus, the dominant behaviour of a spherical Janus swimmer consists in a constant linear velocity along the direction of the chemical gradient, with the orientation of the swimmer being irrelevant. In order to fabricate an efficient chemotactic swimmer of this type, one must simply ensure that $\unicode[STIX]{x1D707}_{s}^{U}+\unicode[STIX]{x1D707}_{s}^{L}$ has the desired sign for positive or negative chemotaxis. The simplest realisation of this is a swimmer with uniform mobility. However, if it happens that $\unicode[STIX]{x1D707}_{s}^{U}+\unicode[STIX]{x1D707}_{s}^{L}=0$ then the chemotactic behaviour of the Janus sphere will depend on first-order effects coming from the chemical reactions. At the end of this subsection we discuss the different possible strategies for chemotaxis.

For a swimmer that responds to both substrate and product molecules it is possible, in theory, to choose the surface properties $\{\unicode[STIX]{x1D70E}^{L},\unicode[STIX]{x1D70E}^{U},\unicode[STIX]{x1D707}_{s}^{L},\unicode[STIX]{x1D707}_{s}^{U},\unicode[STIX]{x1D707}_{p}^{L},\unicode[STIX]{x1D707}_{p}^{U}\}$ in order to obtain any set of prescribed values for $\{U_{Janus},\unicode[STIX]{x1D6FC}_{Janus},\unicode[STIX]{x1D6FD}_{Janus},\unicode[STIX]{x1D6F7}_{Janus}\}$ , although this may be hard to achieve in practice. Surprisingly, the same cannot be said for a swimmer that responds only to product molecules and has $\unicode[STIX]{x1D707}_{s}=0$ , even though there are four degrees of freedom and four target parameters. This is due to the specific combinations in which the four degrees of freedom appear in the expressions for $\{U_{Janus},\unicode[STIX]{x1D6FC}_{Janus},\unicode[STIX]{x1D6FD}_{Janus},\unicode[STIX]{x1D6F7}_{Janus}\}$ . We find that a Janus sphere with $\unicode[STIX]{x1D707}_{s}=0$ is subject to the following dependence relation between the four parameters describing its motion

(2.55) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}_{Janus}=\frac{16\unicode[STIX]{x1D6FE}R^{2}U_{Janus}\unicode[STIX]{x1D6F7}_{Janus}}{3s_{b}\unicode[STIX]{x1D6FD}_{Janus}}+\unicode[STIX]{x1D6FD}_{Janus}.\end{eqnarray}$$

This represents a constraint on the range of behaviour that can be obtained with a swimmer that has zero mobility with respect to the substrate molecules.

For a swimmer that responds to the product molecules, and possibly the substrate molecules as well, we can distinguish different strategies for chemotaxis depending on the quantitative details of the problem. If the chemical gradient is sufficiently strong, the linear velocity of the swimmer will be dominated by the chemotactic ‘sedimentation’ terms, $\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D735}s^{\infty }+\unicode[STIX]{x1D6FD}\,\hat{\boldsymbol{z}}\hat{\boldsymbol{z}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}s^{\infty }$ . In this case it is better to opt for a simpler swimmer with uniform surface properties, which ensures that $\unicode[STIX]{x1D6FD}_{Janus}=0$ and the swimmer undergoes direct chemotactic sedimentation along the gradient, with the sign of $\unicode[STIX]{x1D6FC}_{Janus}$ determining the sense of chemotaxis.

If, on the other hand, the chemical gradient is sufficiently weak, then the linear velocity of the swimmer will be dominated by propulsion along its axis of symmetry, i.e.  $U\hat{\boldsymbol{z}}$ . In this case, an effective chemotactic swimmer must be able to align with the chemical gradient and then propel along its axis, which corresponds to non-zero values for both $U_{Janus}$ and $\unicode[STIX]{x1D6F7}_{Janus}$ . In this situation, it is crucial that the two hemispheres of the Janus sphere have different surface activities as well as different surface mobilities. If $U_{Janus}$ and $\unicode[STIX]{x1D6F7}_{Janus}$ have the same sign then the swimmer will perform positive chemotaxis, whereas if they have opposite signs the swimmer will perform negative chemotaxis. For intermediate values of the chemical gradient, the optimal strategy for chemotaxis will consist of a combination of the two mechanisms.

3 Long-time behaviour: artificial chemotaxis of a population of non-interacting particles

3.1 Motivation

Having established the physical law which describes the instantaneous movement of an individual swimmer, we now want to develop an understanding of the long-time behaviour of a phoretic swimmer, or equivalently that of a large number of non-interacting swimmers. Although we neglect the interaction between particles, we focus on understanding the effects of another essential ingredient – stochasticity. Thermal noise plays an important role on the scale of phoretic swimmers, which are typically a few microns in diameter. The orientation of the swimmers will fluctuate as a result of thermal energy in the surrounding fluid, leading to diffusive behaviour characterised by a rotational diffusivity $D_{r}$ given by the Einstein–Smoluchowski relation (Einstein Reference Einstein1905; von Smoluchowski Reference von Smoluchowski1906)

(3.1) $$\begin{eqnarray}D_{r}=\frac{k_{B}T}{8\unicode[STIX]{x03C0}\unicode[STIX]{x1D702}R^{3}},\end{eqnarray}$$

where $k_{B}$ is Boltzmann’s constant, $T$ is the temperature and $\unicode[STIX]{x1D702}$ is the dynamic viscosity of the fluid.

The physical basis for a macroscopic continuum model lies in the competition between the stochastic effect of rotational diffusion and the deterministic processes by which the swimmer propels and aligns with the external chemical gradient. If we consider the system on a time scale much larger than the characteristic time scale for rotational diffusion, such that we are essentially averaging over fluctuations in the swimmer orientation, we would expect to see the spatial distribution of swimmers evolve in time according to an effective advection–diffusion equation. The purpose of this section is to quantify this process rigorously.

The first modelling approach we considered was based on the classical continuum model for gyrotactic swimming micro-organisms proposed by Pedley & Kessler (Reference Pedley and Kessler1990). This is a simple and intuitive model, but its simplicity comes at a cost. The authors use a phenomenological definition of the diffusivity tensor which requires the introduction of a direction correlation time, also called relaxation time by some authors (see Batchelor Reference Batchelor1976). This quantity is not calculated explicitly in their paper, and would have to be estimated from experiments or numerical simulations. Furthermore, the authors assume that the direction correlation time is isotropic, an assumption probably not suitable for a swimmer that has a preferential orientation.

Therefore, we turned instead our attention to the theory of generalised Taylor dispersion (GTD) that Frankel and Brenner used for the study of orientable Brownian particles (Frankel & Brenner Reference Frankel and Brenner1991, Reference Frankel and Brenner1993) and which has been successfully applied to gyrotactic swimming micro-organisms in more recent years (Hill & Bees Reference Hill and Bees2002; Manela & Frankel Reference Manela and Frankel2003; Bearon, Bees & Croze Reference Bearon, Bees and Croze2012). By addressing the coupling between dynamics in the orientational space and in the physical space, GTD theory provides a rigorous approach to modelling collections of orientable particles.

Fundamentally, GTD theory uses a moment expansion for the probability density of swimmers, in a similar way to Golestanian (Reference Golestanian2012), Pohl & Stark (Reference Pohl and Stark2014) and Bickel et al. (Reference Bickel, Zecua and Würger2014). What is novel in the present paper is that the chemical reactions considered in § 2 lead to a more complex instantaneous behaviour than previously investigated, which brings a further degree of freedom to our continuum model (what we later call the ‘indirect chemotactic index’) and new phenomenology along with it. Furthermore, the diffusivity tensor is isotropic in all of the above references, because none of them possess the necessary combination of active reorientation and orientation-dependent velocity which promotes the emergence of anisotropic diffusion, as is the case with the phoretic swimmers in our paper.

3.2 Set-up

We consider a uniform substrate gradient, as in the derivation of the instantaneous behaviour, and no background flow. Since the orientation of the swimmer changes in time, it no longer makes sense to define one of the principal Cartesian axes along the axis of symmetry of the swimmer. Instead, we take the positive $z$ -axis to be aligned with the fixed chemical gradient such that $\unicode[STIX]{x1D735}s^{\infty }=|\unicode[STIX]{x1D735}s^{\infty }|\boldsymbol{k}$ . We denote the position of the swimmer by $\boldsymbol{x}$ and its axis of symmetry by the unit vector $\boldsymbol{n}=(\sin \unicode[STIX]{x1D703}\cos \unicode[STIX]{x1D719},\sin \unicode[STIX]{x1D703}\sin \unicode[STIX]{x1D719},\cos \unicode[STIX]{x1D703})$ , where $\unicode[STIX]{x1D703}$ and $\unicode[STIX]{x1D719}$ are the usual polar and azimuthal angles from spherical coordinates.

In this notation, the linear velocity of the phoretic swimmer in (2.40) can be written as

(3.2) $$\begin{eqnarray}\boldsymbol{V}(\boldsymbol{x},\boldsymbol{n})=U(\boldsymbol{x})(\boldsymbol{n}+\unicode[STIX]{x1D708}(\boldsymbol{x})\boldsymbol{k}+\unicode[STIX]{x1D707}(\boldsymbol{x})\boldsymbol{n}\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{k}),\end{eqnarray}$$

where we have introduced two dimensionless parameters

(3.3a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D708}(\boldsymbol{x})\equiv \frac{\unicode[STIX]{x1D6FC}|\unicode[STIX]{x1D735}s^{\infty }|}{U(\boldsymbol{x})}\quad \text{and}\quad \unicode[STIX]{x1D707}(\boldsymbol{x})\equiv \frac{\unicode[STIX]{x1D6FD}|\unicode[STIX]{x1D735}s^{\infty }|}{U(\boldsymbol{x})}.\end{eqnarray}$$

We call $\unicode[STIX]{x1D708}$ and $\unicode[STIX]{x1D707}$ the direct and indirect chemotactic indices, respectively, since they measure the relative effectiveness of direct chemotactic sedimentation, $\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D735}s^{\infty }$ , and indirect chemotactic sedimentation, $\unicode[STIX]{x1D6FD}\boldsymbol{n}\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}s^{\infty }$ , against active propulsion, $U(\boldsymbol{x})\boldsymbol{n}$ .

Note that the relative magnitude $\unicode[STIX]{x1D708}/\unicode[STIX]{x1D707}$ of direct to indirect chemotactic sedimentation determines the angle at which the phoretic swimmer sediments relative to the chemical gradient. This is similar to the Stokes flow sedimentation of an elongated object, such as a rod, under the action of gravity. In Stokes flow the sedimentation angle is given by the geometry (and the drag coefficients) of the elongated body, while in our case it is a function of the chemical properties of the spherical swimmer.

Lastly, we need to consider the mechanism by which the swimmer aligns with the chemical gradient. The orientation of the swimmer evolves in time according to

(3.4) $$\begin{eqnarray}\dot{\boldsymbol{n}}=\unicode[STIX]{x1D74E}\wedge \boldsymbol{n},\end{eqnarray}$$

and using (2.41) we arrive at the following reorientation law

(3.5) $$\begin{eqnarray}\dot{\boldsymbol{n}}=\unicode[STIX]{x1D6FA}[\boldsymbol{k}-(\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{n})\boldsymbol{n}],\end{eqnarray}$$

where $\unicode[STIX]{x1D6FA}=\unicode[STIX]{x1D6F7}|\unicode[STIX]{x1D735}s^{\infty }|$ . This prompts us to define a third dimensionless parameter

(3.6) $$\begin{eqnarray}\unicode[STIX]{x1D706}\equiv \frac{\unicode[STIX]{x1D6FA}}{D_{r}},\end{eqnarray}$$

which is a rotational Péclet number measuring the relative importance of active reorientation to rotational diffusion.

3.3 GTD theory

Although GTD theory is now a classical tool, it is useful to introduce it clearly and define the differential operator involved in it. We start from the conservation equation for $P(\boldsymbol{n},\boldsymbol{x},t)$ , the probability density function of finding a swimmer with orientation $\boldsymbol{n}$ at position $\boldsymbol{x}$ at time $t$ , which satisfies the conservation equation

(3.7) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}P}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}_{r}\boldsymbol{\cdot }(\boldsymbol{V}(\boldsymbol{x},\boldsymbol{n})P-D_{t}\unicode[STIX]{x1D735}_{r}P)+\unicode[STIX]{x1D735}_{n}\boldsymbol{\cdot }(\dot{\boldsymbol{n}}P-D_{r}\unicode[STIX]{x1D735}_{n}P)=0,\end{eqnarray}$$

where the differential operators for physical space and orientational space are

(3.8a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D735}_{r}\equiv \boldsymbol{i}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{1}}+\boldsymbol{j}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{2}}+\boldsymbol{k}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{3}},\quad \unicode[STIX]{x1D735}_{n}\equiv \hat{\unicode[STIX]{x1D73D}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}+\hat{\unicode[STIX]{x1D753}}\frac{1}{\sin \unicode[STIX]{x1D703}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}},\end{eqnarray}$$

respectively. The law in (3.7) says that the local rate of change in the distribution of swimmers is due to the flux of swimmers in physical space (advection and translational diffusion) and the flux of swimmers in orientational space (advection and rotational diffusion).

Since we are interested in the spatial distribution of the particles, we define the particle density

(3.9) $$\begin{eqnarray}\unicode[STIX]{x1D70C}(\boldsymbol{x},t)\equiv \int _{S_{2}}P(\boldsymbol{n},\boldsymbol{x},t)\,\text{d}^{2}\boldsymbol{n}.\end{eqnarray}$$

Upon careful consideration of the moments of $P(\boldsymbol{n},\boldsymbol{x},t)$ , Frankel and Brenner showed that at the macroscopic level, and on time scales $t\gg D_{r}^{-1}$ , a population of non-interacting orientable particles obeys the advection–diffusion equation

(3.10) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}_{\boldsymbol{r}}\boldsymbol{\cdot }(\boldsymbol{V}_{s}\unicode[STIX]{x1D70C}-\unicode[STIX]{x1D63F}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\boldsymbol{r}}\unicode[STIX]{x1D70C})=0,\end{eqnarray}$$

where $\boldsymbol{V}_{s}$ is the mean swimming velocity and $\unicode[STIX]{x1D63F}$ is the diffusivity tensor. The mean swimming velocity is given by

(3.11) $$\begin{eqnarray}\boldsymbol{V}_{s}(\boldsymbol{x})=\int _{S_{2}}f(\boldsymbol{n})\boldsymbol{V}(\boldsymbol{x},\boldsymbol{n})\,\text{d}^{2}\boldsymbol{n},\end{eqnarray}$$

whereas the diffusivity tensor is the sum of contributions from passive translational diffusion ( $\unicode[STIX]{x1D63F}_{pass}$ ) and active phoretic motion ( $\unicode[STIX]{x1D63F}_{act}$ ) as

(3.12) $$\begin{eqnarray}\unicode[STIX]{x1D63F}(\boldsymbol{x})=\unicode[STIX]{x1D63F}_{pass}+\unicode[STIX]{x1D63F}_{act}(\boldsymbol{x}).\end{eqnarray}$$

The passive part of the diffusivity tensor, $\unicode[STIX]{x1D63F}_{pass}$ , is simply the average of the translational diffusivity tensor, $\unicode[STIX]{x1D63F}_{t}(\boldsymbol{n})$ , over the space of orientations (see Frankel & Brenner Reference Frankel and Brenner1993). For a spherical swimmer it is reasonable to assume that Brownian noise acts isotropically, and so $\unicode[STIX]{x1D63F}_{pass}=D_{t}\unicode[STIX]{x1D644}$ with $D_{t}$ being the translational diffusion coefficient. Recently it has been shown by Agudo-Canalejo, Illien & Golestanian (Reference Agudo-Canalejo, Illien and Golestanian2018) that, at the nano-scale, spatial variations in the translational diffusion coefficient of enzymes can compete with phoretic effects, and it may well be possible that the catalytic coating of phoretic swimmers could lead to enhanced diffusion coefficients as the swimmer moves through the chemical gradient. However, due to a lack of documentation on the enhanced diffusion capabilities of a microscopic phoretic swimmer, and due to the limitations of GTD theory, to consider spatial variations of $D_{t}$ would go beyond the scope of this paper.

With the assumption that $\unicode[STIX]{x1D63F}_{pass}$ is constant and independent of the phoretic mechanisms for motion, we will focus from here onwards on quantifying the active diffusivity tensor given by

(3.13) $$\begin{eqnarray}\unicode[STIX]{x1D63F}_{act}(\boldsymbol{x})=\int _{S_{2}}[\boldsymbol{b}(\boldsymbol{x},\boldsymbol{n})\boldsymbol{V}(\boldsymbol{x},\boldsymbol{n})]^{sym}\,\text{d}^{2}\boldsymbol{n},\end{eqnarray}$$

where $[\cdots ]^{sym}$ denotes the symmetric part of the tensor.

The expressions for $\boldsymbol{V}_{s}(\boldsymbol{x})$ and $\unicode[STIX]{x1D63F}_{act}(\boldsymbol{x})$ contain two new fields that have to be computed: $f(\boldsymbol{n})$ is the long-time steady distribution of the particle orientation, and $\boldsymbol{b}(\boldsymbol{x},\boldsymbol{n})$ represents the long-time relative displacement of a particle given that its instantaneous orientation is $\boldsymbol{n}$ , compared to its expected position over all possible orientations. For a formal definition of these quantities we need to introduce the linear operator

(3.14) $$\begin{eqnarray}\text{L}(\star )=\unicode[STIX]{x1D735}_{\boldsymbol{n}}\boldsymbol{\cdot }[\dot{\boldsymbol{n}}(\star )-D_{r}\unicode[STIX]{x1D735}_{\boldsymbol{n}}(\star )],\end{eqnarray}$$

which describes the flux of a quantity in orientational space. Then $f$ and $\boldsymbol{b}$ are the solutions to the following equations and boundary conditions

(3.15a,b ) $$\begin{eqnarray}\text{L}f=0,\quad \int _{S_{2}}f(\boldsymbol{n})\,\text{d}^{2}\boldsymbol{n}=1,\end{eqnarray}$$
(3.16a,b ) $$\begin{eqnarray}\text{L}\boldsymbol{b}=f(\boldsymbol{n})(\boldsymbol{V}(\boldsymbol{x},\boldsymbol{n})-\boldsymbol{V}_{s}(\boldsymbol{x})),\quad \int _{S_{2}}\boldsymbol{b}\,\text{d}^{2}\boldsymbol{n}=0.\end{eqnarray}$$

For a full derivation of these results see Frankel & Brenner (Reference Frankel and Brenner1993) or, for a more succinct account, Hill & Bees (Reference Hill and Bees2002).

To provide some physical intuition for these two equations, we observe that equation (3.15) corresponds to solving $\unicode[STIX]{x2202}f/\unicode[STIX]{x2202}t=0$ and therefore finding the steady-state distribution of orientations subject to a global normalisation condition. Similarly, equation (3.16) can be thought of as solving $\unicode[STIX]{x2202}\boldsymbol{b}/\unicode[STIX]{x2202}t=f(\boldsymbol{n})(\boldsymbol{V}(\boldsymbol{x},\boldsymbol{n})-\boldsymbol{V}_{s}(\boldsymbol{x}))$ which means that the rate of change of the relative displacement is due to the relative velocity, weighted by the probability of finding a swimmer with that specific orientation. The obvious boundary condition to impose here is that, averaged over all possible orientations, the relative displacement must vanish.

One important difference between our application of GTD theory to phoretic swimmers, compared to previous applications to orientable Brownian particles and gyrotactic micro-swimmers, is that the swimming velocity $\boldsymbol{V}(\boldsymbol{x},\boldsymbol{n})$ of phoretic swimmers depends on position as well. The classical theory still applies provided that the velocity is slowly varying in space, such that the process of averaging over particle orientations can be carried out on time scales much larger than the diffusive time scale, $D_{r}^{-1}$ , but much smaller than the time scale of travel to a region of significantly different swimming speed. So there must exist a time scale $\unicode[STIX]{x1D70F}$ , satisfying

(3.17) $$\begin{eqnarray}D_{r}^{-1}\ll \unicode[STIX]{x1D70F}\ll |\unicode[STIX]{x1D735}U|^{-1},\end{eqnarray}$$

on which we can average (3.7) over the space of swimmer orientations. In particular, we must have

(3.18) $$\begin{eqnarray}|\unicode[STIX]{x1D735}U|\ll D_{r}.\end{eqnarray}$$

In that case, the advection–diffusion equation (3.10) will be valid on time scales greater than or equal to $\unicode[STIX]{x1D70F}$ and will involve a mean swimming velocity, $\boldsymbol{V}_{s}(\boldsymbol{x})$ , and a diffusivity tensor, $\unicode[STIX]{x1D63F}(\boldsymbol{x})$ , which are slowly varying in space. Under this approximation, quantities such as $\unicode[STIX]{x1D708}(\boldsymbol{x}),\unicode[STIX]{x1D707}(\boldsymbol{x}),\boldsymbol{V}(\boldsymbol{x},\boldsymbol{n})$ and $\boldsymbol{b}(\boldsymbol{x},\boldsymbol{n})$ are also slowly varying in space, and we will drop the explicit mention of spatial dependence from our notation.

3.4 Distribution of swimmer orientations

Since the orientation of the swimmers fluctuates over time from thermal noise, the only fixed meaningful direction in our problem is that of the chemical gradient. Therefore, quantities such as $f(\boldsymbol{n})$ and $\boldsymbol{b}(\boldsymbol{n})$ will inherit rotational symmetry around the direction of the chemical gradient (the $\boldsymbol{k}$ axis). For such an axisymmetric problem, the linear operator simplifies to

(3.19) $$\begin{eqnarray}\text{L}(\star )=-\frac{\unicode[STIX]{x1D6FA}}{\sin \unicode[STIX]{x1D703}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}(\sin ^{2}\unicode[STIX]{x1D703}(\star ))-\frac{D_{r}}{\sin \unicode[STIX]{x1D703}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}\left(\sin \unicode[STIX]{x1D703}\frac{\unicode[STIX]{x2202}(\star )}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}\right).\end{eqnarray}$$

Then (3.15a ) has a trivial first integral

(3.20) $$\begin{eqnarray}\unicode[STIX]{x1D706}\sin ^{2}\unicode[STIX]{x1D703}f+\sin \unicode[STIX]{x1D703}\frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}=0,\end{eqnarray}$$

where the integrating constant can be set to zero assuming that $f$ is well behaved at the poles. The solution to this, subject to the appropriate normalisation condition, is the same as that found by Pedley & Kessler (Reference Pedley and Kessler1990). The long-time steady distribution of the swimmer orientation is

(3.21) $$\begin{eqnarray}f(\boldsymbol{n})=\hat{f}(\cos \unicode[STIX]{x1D703})=\frac{\unicode[STIX]{x1D706}\text{e}^{\unicode[STIX]{x1D706}\cos \unicode[STIX]{x1D703}}}{4\unicode[STIX]{x03C0}\sinh \unicode[STIX]{x1D706}},\end{eqnarray}$$

where $\unicode[STIX]{x1D706}$ is the rotational Péclet number defined in (3.6).

3.5 Mean swimming velocity

From (3.2) and (3.11), the mean swimming velocity is obtained as

(3.22) $$\begin{eqnarray}\boldsymbol{V}_{s}=U\left(\int _{S_{2}}f(\boldsymbol{n})\boldsymbol{n}\,\text{d}^{2}\boldsymbol{n}+\unicode[STIX]{x1D708}\int _{S_{2}}f(\boldsymbol{n})\boldsymbol{k}\,\text{d}^{2}\boldsymbol{n}+\unicode[STIX]{x1D707}\int _{S_{2}}f(\boldsymbol{n})\boldsymbol{n}\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{k}\,\text{d}^{2}\boldsymbol{n}\right).\end{eqnarray}$$

We use the normalisation of $f(\boldsymbol{n})$ to simplify the second integral, and average over the azimuthal angle in the other two integrals in order to get

(3.23) $$\begin{eqnarray}\boldsymbol{V}_{s}=U\left(\int _{0}^{\unicode[STIX]{x03C0}}\hat{f}(\cos \unicode[STIX]{x1D703})\cos \unicode[STIX]{x1D703}\sin \unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D703}+\unicode[STIX]{x1D708}+\unicode[STIX]{x1D707}\int _{0}^{\unicode[STIX]{x03C0}}\hat{f}(\cos \unicode[STIX]{x1D703})\cos ^{2}\unicode[STIX]{x1D703}\sin \unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D703}\right)\boldsymbol{k}.\end{eqnarray}$$

Using our previous expression for $\hat{f}(\cos \unicode[STIX]{x1D703})$ , we obtain the following analytical expression for the mean swimming velocity

(3.24) $$\begin{eqnarray}\boldsymbol{V}_{s}=U\left[\frac{-1+\unicode[STIX]{x1D706}\coth \unicode[STIX]{x1D706}}{\unicode[STIX]{x1D706}}+\unicode[STIX]{x1D708}+\unicode[STIX]{x1D707}\left(\frac{2+\unicode[STIX]{x1D706}^{2}-2\unicode[STIX]{x1D706}\coth \unicode[STIX]{x1D706}}{\unicode[STIX]{x1D706}^{2}}\right)\right]\boldsymbol{k}.\end{eqnarray}$$

Note that the linear velocity of the swimmer relative to the mean swimming velocity is a function of $\unicode[STIX]{x1D706}$ and $\unicode[STIX]{x1D707}$ only since

(3.25) $$\begin{eqnarray}\boldsymbol{V}(\boldsymbol{n})-\boldsymbol{V}_{s}=U[\boldsymbol{n}+\unicode[STIX]{x1D707}\boldsymbol{n}\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{k}-\tilde{u} (\unicode[STIX]{x1D706},\unicode[STIX]{x1D707})\boldsymbol{k}],\end{eqnarray}$$

where we have defined for convenience

(3.26) $$\begin{eqnarray}\tilde{u} (\unicode[STIX]{x1D706},\unicode[STIX]{x1D707})\equiv \frac{-1+\unicode[STIX]{x1D706}\coth \unicode[STIX]{x1D706}}{\unicode[STIX]{x1D706}}+\unicode[STIX]{x1D707}\left(\frac{2+\unicode[STIX]{x1D706}^{2}-2\unicode[STIX]{x1D706}\coth \unicode[STIX]{x1D706}}{\unicode[STIX]{x1D706}^{2}}\right).\end{eqnarray}$$

As a result, both $\boldsymbol{b}(\boldsymbol{n})$ and the diffusivity tensor $\unicode[STIX]{x1D63F}_{act}$ will only depend on the value of $\unicode[STIX]{x1D706}$ and $\unicode[STIX]{x1D707}$ . Physically, the direct chemotactic index $\unicode[STIX]{x1D708}$ corresponds to a component of velocity that is constant in time, which imparts a constant drift to the phoretic swimmer without affecting the way in which the swimmer diffuses through space.

3.6 Active diffusivity tensor

The first step in calculating the active diffusivity tensor is to solve for the vector field $\boldsymbol{b}(\boldsymbol{n})$ . The expression for the relative swimming velocity from (3.25) together with the definition of $\boldsymbol{b}(\boldsymbol{n})$ from (3.16) suggest that the only dependence on the azimuthal angle that can be excited in $\boldsymbol{b}(\boldsymbol{n})$ is a $\cos \unicode[STIX]{x1D719}$ dependence in its $x$ -component and a $\sin \unicode[STIX]{x1D719}$ dependence in its $y$ -component. Therefore, we look to solve for $\boldsymbol{b}(\boldsymbol{n})$ in the form

(3.27) $$\begin{eqnarray}\boldsymbol{b}(\boldsymbol{n})=(b_{\bot }(\unicode[STIX]{x1D703})\cos \unicode[STIX]{x1D719},b_{\bot }(\unicode[STIX]{x1D703})\sin \unicode[STIX]{x1D719},b_{\Vert }(\unicode[STIX]{x1D703})),\end{eqnarray}$$

where the function $b_{\bot }(\unicode[STIX]{x1D703})$ is the same in the $x$ and $y$ -components due to rotational symmetry around the direction of the chemical gradient. We then decompose

(3.28) $$\begin{eqnarray}\displaystyle & \displaystyle b_{\bot }(\unicode[STIX]{x1D703})=\frac{U}{D_{r}}\mathop{\sum }_{l=1}^{\infty }-b_{l}^{\bot }P_{l}^{1}(\cos \unicode[STIX]{x1D703})=\frac{U}{D_{r}}\mathop{\sum }_{l=1}^{\infty }b_{l}^{\bot }\sin \unicode[STIX]{x1D703}P_{l}^{\prime }(\cos \unicode[STIX]{x1D703}), & \displaystyle\end{eqnarray}$$
(3.29) $$\begin{eqnarray}\displaystyle & \displaystyle b_{\Vert }(\unicode[STIX]{x1D703})=\frac{U}{D_{r}}\mathop{\sum }_{l=0}^{\infty }b_{l}^{\Vert }P_{l}(\cos \unicode[STIX]{x1D703}), & \displaystyle\end{eqnarray}$$

and substitute into (3.16) to obtain, after some simplification, the equalities

(3.30) $$\begin{eqnarray}\displaystyle & \displaystyle \mathop{\sum }_{l=1}^{\infty }b_{l}^{\bot }[(l(l+1)-3\unicode[STIX]{x1D706}\cos \unicode[STIX]{x1D703})P_{l}^{\prime }+\unicode[STIX]{x1D706}\sin ^{2}\unicode[STIX]{x1D703}P_{l}^{\prime \prime }]=(1+\unicode[STIX]{x1D707}\cos \unicode[STIX]{x1D703})\hat{f}(\cos \unicode[STIX]{x1D703}), & \displaystyle\end{eqnarray}$$
(3.31) $$\begin{eqnarray}\displaystyle & \displaystyle \mathop{\sum }_{l=1}^{\infty }b_{l}^{\Vert }[(l(l+1)-2\unicode[STIX]{x1D706}\cos \unicode[STIX]{x1D703})P_{l}+\unicode[STIX]{x1D706}\sin ^{2}\unicode[STIX]{x1D703}P_{l}^{\prime }]=(\cos \unicode[STIX]{x1D703}(1+\unicode[STIX]{x1D707}\cos \unicode[STIX]{x1D703})-\tilde{u} )\hat{f}(\cos \unicode[STIX]{x1D703}).\qquad & \displaystyle\end{eqnarray}$$

Note that $b_{0}^{\Vert }=0$ results from imposing $\int _{S_{2}}\boldsymbol{b}\,\text{d}^{2}\boldsymbol{n}=0$ .

The most important step in the derivation is choosing the appropriate inner product for these two equations, so that the integrals on the left-hand side may be evaluated exactly. In (3.30) we make the substitution $\unicode[STIX]{x1D709}=\cos \unicode[STIX]{x1D703}$ and take $\int \ldots (1-\unicode[STIX]{x1D709}^{2})P_{m}(\unicode[STIX]{x1D709})\,\text{d}\unicode[STIX]{x1D709}$ of both sides. The left-hand side can then be evaluated using the identities (B 8), (B 10) and (B 11). We obtain a set of linear equations for the coefficients $b_{l}^{\bot }$ as

(3.32) $$\begin{eqnarray}\unicode[STIX]{x1D617}_{ml}b_{l}^{\bot }=\int _{-1}^{+1}(1+\unicode[STIX]{x1D707}\unicode[STIX]{x1D709})\hat{f}(\unicode[STIX]{x1D709})(1-\unicode[STIX]{x1D709}^{2})P_{m}(\unicode[STIX]{x1D709})\,\text{d}\unicode[STIX]{x1D709},\quad (m\geqslant 0)\end{eqnarray}$$

where the pentadiagonal matrix $\unicode[STIX]{x1D64B}$ has components

(3.33) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D617}_{ml} & = & \displaystyle \frac{2\unicode[STIX]{x1D706}(m+3)(m+2)(m+1)^{2}}{(2m+5)(2m+3)(2m+1)}\unicode[STIX]{x1D6FF}_{l,m+2}+\frac{2(m+2)^{2}(m+1)^{2}}{(2m+3)(2m+1)}\unicode[STIX]{x1D6FF}_{l,m+1}\nonumber\\ \displaystyle & & \displaystyle -\,\frac{2\unicode[STIX]{x1D706}m(m+1)(2m^{2}+2m-1)}{(2m+3)(2m+1)(2m-1)}\unicode[STIX]{x1D6FF}_{lm}+\frac{2m^{2}(m-1)^{2}}{(2m+1)(2m-1)}\unicode[STIX]{x1D6FF}_{l,m-1}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{2\unicode[STIX]{x1D706}m^{2}(m-1)(m-2)}{(2m+1)(2m-1)(2m-3)}\unicode[STIX]{x1D6FF}_{l,m-2}.\end{eqnarray}$$

In (3.31) we make the same substitution and take $\int \ldots P_{m}(\unicode[STIX]{x1D709})\,\text{d}\unicode[STIX]{x1D709}$ of both sides. Using identities (B 2), (B 6) and (B 8) we obtain another set of linear equations, this time for the coefficients $b_{l}^{\Vert }$

(3.34) $$\begin{eqnarray}\unicode[STIX]{x1D61B}_{ml}b_{l}^{\Vert }=\int _{-1}^{+1}\left[\unicode[STIX]{x1D709}(1+\unicode[STIX]{x1D707}\unicode[STIX]{x1D709})-\tilde{u} \right]\hat{f}(\unicode[STIX]{x1D709})P_{m}(\unicode[STIX]{x1D709})\,\text{d}\unicode[STIX]{x1D709},\quad (m\geqslant 0),\end{eqnarray}$$

where the tridiagonal matrix $\unicode[STIX]{x1D64F}$ has components

(3.35) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D61B}_{ml}=\frac{2\unicode[STIX]{x1D706}m(m+1)}{(2m+1)(2m+3)}\unicode[STIX]{x1D6FF}_{l,m+1}+\frac{2m(m+1)}{2m+1}\unicode[STIX]{x1D6FF}_{lm}-\frac{2\unicode[STIX]{x1D706}m(m+1)}{(2m-1)(2m+1)}\unicode[STIX]{x1D6FF}_{l,m-1}. & & \displaystyle\end{eqnarray}$$

Equations (3.32) and (3.34) represent infinite systems of linear equations for the coefficients $b_{m}^{\bot }$ and $b_{m}^{\Vert }$ which cannot be inverted analytically. Progress can be made by truncating the series for $b_{\bot }(\unicode[STIX]{x1D703})$ and $b_{\Vert }(\unicode[STIX]{x1D703})$ to $N$ terms and inverting the resulting $N$ -by- $N$ systems numerically. Rather fortunately, only the first two coefficients are needed to compute the active diffusivity tensor, as we shall now see.

If we substitute our expressions for $\boldsymbol{V}(\boldsymbol{n})$ and $\boldsymbol{b}(\boldsymbol{n})$ into the definition of the active diffusivity tensor, (3.13), and average over the azimuthal angle, we find that all the off-diagonal terms vanish. Furthermore, we have $(\unicode[STIX]{x1D63F}_{act})_{11}=(\unicode[STIX]{x1D63F}_{act})_{22}$ due to rotational symmetry about the $\boldsymbol{k}$ axis, which is the direction of the chemical gradient. We denote $D_{\Vert }\equiv (\unicode[STIX]{x1D63F}_{act})_{33}$ to be the active diffusivity parallel to the chemical gradient, and $D_{\bot }\equiv (\unicode[STIX]{x1D63F}_{act})_{11}=(\unicode[STIX]{x1D63F}_{act})_{22}$ to be the active diffusivity in the plane normal to it. From (3.13) these quantities are

(3.36) $$\begin{eqnarray}\displaystyle & \displaystyle D_{\bot }=\unicode[STIX]{x03C0}\int _{0}^{\unicode[STIX]{x03C0}}b_{\bot }(\unicode[STIX]{x1D703})U\sin \unicode[STIX]{x1D703}(1+\unicode[STIX]{x1D707}\cos \unicode[STIX]{x1D703})\sin \unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D703}, & \displaystyle\end{eqnarray}$$
(3.37) $$\begin{eqnarray}\displaystyle & \displaystyle D_{\Vert }=2\unicode[STIX]{x03C0}\int _{0}^{\unicode[STIX]{x03C0}}b_{\Vert }(\unicode[STIX]{x1D703})U\cos \unicode[STIX]{x1D703}(1+\unicode[STIX]{x1D707}\cos \unicode[STIX]{x1D703})\sin \unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D703}. & \displaystyle\end{eqnarray}$$

Using identities (B 5), (B 6), (B 8) we can simplify these expression to finally obtain

(3.38) $$\begin{eqnarray}\displaystyle & \displaystyle D_{\bot }=\frac{4\unicode[STIX]{x03C0}U^{2}}{D_{r}}\left(\frac{b_{1}^{\bot }}{3}+\frac{\unicode[STIX]{x1D707}b_{2}^{\bot }}{5}\right), & \displaystyle\end{eqnarray}$$
(3.39) $$\begin{eqnarray}\displaystyle & \displaystyle D_{\Vert }=\frac{4\unicode[STIX]{x03C0}U^{2}}{D_{r}}\left(\frac{b_{1}^{\Vert }}{3}+\frac{2\unicode[STIX]{x1D707}b_{2}^{\Vert }}{15}\right). & \displaystyle\end{eqnarray}$$

Since only the first two coefficients $b_{m}^{\bot }$ and $b_{m}^{\Vert }$ ( $m=1,2$ ) are needed in the calculation of the diffusivity tensor, the truncated systems converge very rapidly and we find that it is usually sufficient to truncate to $N$ of order $O(\unicode[STIX]{x1D706},\unicode[STIX]{x1D707})$ or $N=2$ , whichever is larger.

3.7 Numerical simulations

In order to validate the continuum model we performed numerical simulations of phoretic swimmers using a simple Euler–Maruyama scheme. To capture the correct dynamics we use a time step much smaller than the time scale $D_{r}^{-1}$ for rotational diffusion, and calculate all relevant macroscopic quantities as averages over a time period larger than $D_{r}^{-1}$ . At each time step, the position and orientation of the swimmers are updated with the deterministic contribution from the linear and angular velocities of the swimmer, as in the classical Euler method. In addition to this, the orientation of the swimmer experiences a stochastic contribution due to thermal noise. During the $n$ th time step the axis of symmetry of the swimmer is deflected in a random direction by an angle $\unicode[STIX]{x0394}\unicode[STIX]{x1D703}_{n}=2\sqrt{D_{r}\unicode[STIX]{x0394}t}\unicode[STIX]{x0394}W_{n}$ , where $\unicode[STIX]{x0394}W_{n}$ are independent random variables taken from the standard normal distribution (Saragosti, Silberzan & Buguin Reference Saragosti, Silberzan and Buguin2012).

Figure 3. Comparison of the distribution of swimmer orientations obtained from numerical simulations with the predictions of the continuum model. (a) Histograms of final orientations for 10 000 swimmers whose initial orientations were chosen from a uniform distribution, after being evolved in time for 500 time steps of size $\unicode[STIX]{x0394}t=10^{-2}D_{r}^{-1}$ . The red solid line corresponds to the analytical expression from (3.21) for the steady-state distribution of swimmer orientations. (b) Sample of 1000 swimmer orientations chosen from the steady-state distribution. Each swimmer orientation corresponds to a unit vector which is represented here as a blue dot on the unit sphere. From top to bottom the rotational Péclet number is $\unicode[STIX]{x1D706}=0.1,3,10$ .

The only parameters that enter the simulation are the phenomenological constants $\unicode[STIX]{x1D707}$ and $\unicode[STIX]{x1D706}$ which describe deterministic translation and reorientation, while the stochastic reorientation time scale $D_{r}^{-1}$ is normalised to one unit of time in our simulations. Since the direct chemotactic index $\unicode[STIX]{x1D708}$ contributes to a constant linear velocity of the particle, we can take it to be zero in our simulations, which means that we are effectively moving to a frame of reference that is sedimenting with the particle, translating at a constant speed $\unicode[STIX]{x1D6FC}|\unicode[STIX]{x1D735}s^{\infty }|$ in the direction of the chemical gradient. We only consider a non-zero value for $\unicode[STIX]{x1D708}$ in figure 6 when we normalise the active diffusivity coefficients (derived from theory only, no simulations) using the mean swimming velocity. In this case we take $\unicode[STIX]{x1D708}$ to be an arbitrary non-zero value in order to avoid the singularity at $\unicode[STIX]{x1D706},\unicode[STIX]{x1D707}=0$ where the mean swimming velocity would also go to zero otherwise.

To validate our numerical method we first compared the steady-state distribution of orientations for our simulated swimmers against the analytical expression in (3.21), which is a well-established result in the literature (Pedley & Kessler Reference Pedley and Kessler1990; Bearon et al. Reference Bearon, Bees and Croze2012). This validation is illustrated in figure 3, which shows perfect agreement between numerical simulations and theory.

Figure 4. Time evolution of 200 independent trajectories, all starting from the origin at $t=0$ . Each snapshot is taken after a time interval of $3D_{r}^{-1}$ , and the chemical gradient is along the horizontal axis. The large dots are placed at regular distances of $3V_{s}D_{r}^{-1}$ from the origin which confirms that the cloud of phoretic swimmers drifts along the direction of the chemical gradient at the mean swimming velocity ( $V_{s}$ ), as predicted by the continuum model. This simulation corresponds to parameters $\unicode[STIX]{x1D706}=5$ and $\unicode[STIX]{x1D707}=3$ , and the time step used was $\unicode[STIX]{x0394}t=10^{-2}D_{r}^{-1}$ .

We then investigated qualitatively the behaviour of a cloud of swimmers initially starting from the origin. These results are depicted in figure 4. The cloud of swimmers is seen to diffuse anisotropically and to drift along the direction of the chemical gradient at the mean swimming velocity, $V_{s}$ , as predicted by the continuum model.

Figure 5. Comparison of the values obtained from numerical simulations against those predicted by GTD theory for the dependence of component $D_{\bot }$ of the active diffusivity tensor on $\unicode[STIX]{x1D706}$ . The simulations were run with 100 samples of 1000 swimmers using a time step $\unicode[STIX]{x0394}t=10^{-2}D_{r}^{-1}$ over a period $2D_{r}^{-1}$ . Each sample was used to calculate one value of the diffusivity $D_{\bot }$ , and then averaged to obtain one data point. The error bars represent one standard deviation amongst the 100 values obtained for the diffusivity. The solid lines represent the predictions of the continuum model from (3.32)–(3.35) and (3.38)–(3.39). In all plots the diffusivity is non-dimensionalised by $U^{2}D_{r}^{-1}$ . The four graphs correspond to parameter values $\unicode[STIX]{x1D707}=0$ , $\unicode[STIX]{x1D707}=1$ , $\unicode[STIX]{x1D707}=5$ and $\unicode[STIX]{x1D707}=10$ .

The final comparison focuses on the values of the active diffusivity tensor computed from stochastic simulations against those given by GTD theory. We plot in figure 5 four comparative graphs for the dependence of $D_{\bot }$ on $\unicode[STIX]{x1D706}$ at different values of the indirect chemotactic index $\unicode[STIX]{x1D707}$ (results for $D_{\Vert }$ , not shown, are very similar). In all cases, the agreement between the simulations and the theory is excellent and the data points obtained numerically are always within one standard deviation of the theoretical prediction.

Figure 6. Iso-values of the diffusivities $D_{\bot }$ and $D_{\Vert }$ computed using GTD theory in the $(\unicode[STIX]{x1D706},\unicode[STIX]{x1D707})$ plane. Both diffusivities are normalised using the mean swimming velocity and the time scale of rotational diffusion, i.e.  $D\mapsto D/(V_{s}^{2}D_{r}^{-1})$ . These contour plots are obtained with a value of $\unicode[STIX]{x1D708}=1.5$ for the direct chemotactic index, but there are no qualitative changes as $\unicode[STIX]{x1D708}$ is varied (not shown). We observe consistent decay of both diffusivities as $\unicode[STIX]{x1D706}\rightarrow \infty$ at a fixed value of $\unicode[STIX]{x1D707}$ , as well as a saddle point in diffusivity for $\unicode[STIX]{x1D706},\unicode[STIX]{x1D707}=O(1)$ .

We can also use the continuum model to understand how the active diffusivity tensor changes as we vary the values of the rotational Péclet number, $\unicode[STIX]{x1D706}$ , and the indirect chemotactic index, $\unicode[STIX]{x1D707}$ . We plot in figure 6 iso-values of the two components of the active diffusivity tensor normalised by $V_{s}^{2}D_{r}^{-1}$ , with $V_{s}$ being the magnitude of the mean swimming velocity. The active diffusivity decays as the rotational Péclet number $\unicode[STIX]{x1D706}\rightarrow \infty$ for any fixed value of $\unicode[STIX]{x1D707}$ , as expected, because the active reorientation of the particle enhances directed motion along the chemical gradient and suppresses diffusion. We also observe a local maximum in diffusion at $\unicode[STIX]{x1D706},\unicode[STIX]{x1D707}=0$ and a saddle point for $\unicode[STIX]{x1D706},\unicode[STIX]{x1D707}$ of $O(1)$ . The surprising finding is that diffusivity increases as $\unicode[STIX]{x1D707}\rightarrow \infty$ for a fixed value of $\unicode[STIX]{x1D706}$ , meaning that indirect chemotactic sedimentation tends to augment diffusion relative to the mean drift velocity of the swimmer.

Finally, we investigate how the anisotropy of the diffusion depends on the parameters of the problem; specifically we measure the ratio $D_{\bot }/D_{\Vert }$ of the diffusivity in the plane normal to the chemical gradient to the diffusivity parallel to the gradient. Iso-values of this ratio in the $(\unicode[STIX]{x1D706},\unicode[STIX]{x1D707})$ plane are shown in figure 7. We find that the contour $D_{\bot }/D_{\Vert }=1$ passes through the origin, in agreement with our expectation to recover isotropic diffusion in the limit where both chemotactic alignment and indirect chemotactic sedimentation are weak. We observe that the ratio $D_{\bot }/D_{\Vert }$ falls below unity in the region of the $(\unicode[STIX]{x1D706},\unicode[STIX]{x1D707})$ plane above this contour and rises above it in the region to the right of the contour. This indicates that chemotactic alignment, whose strength is given by $\unicode[STIX]{x1D706}$ , favours the reduction of diffusion along the chemical gradient. On the other hand, the indirect chemotactic sedimentation quantified by $\unicode[STIX]{x1D707}$ favours the reduction of diffusion in the plane normal to the chemical gradient.

From figure 7 we also note that $D_{\bot }/D_{\Vert }>1$ for the parameters $\unicode[STIX]{x1D706}=5$ and $\unicode[STIX]{x1D707}=3$ which were used in our analysis of the time evolution of a cloud of swimmers. This is apparent in figure 4 where the cloud of swimmers, initially at the origin, grows anisotropically over time, with more spreading being observed in the vertical than in the horizontal direction, the latter being the direction of the chemical gradient.

3.8 Asymptotics

For small $\unicode[STIX]{x1D706}$ and $\unicode[STIX]{x1D707}$ we can invert the 2-by-2 systems obtained from the truncation of the linear systems of equations (3.32) and (3.34) exactly, in order to obtain the approximate analytical expressions for $D_{\bot }$ and $D_{\Vert }$ that are given in appendix C. We then take the limit $\unicode[STIX]{x1D706},\unicode[STIX]{x1D707}\rightarrow 0$ in these expressions and find that

(3.40a,b ) $$\begin{eqnarray}D_{\bot }\rightarrow \frac{U^{2}}{6D_{r}},\quad D_{\Vert }\rightarrow \frac{U^{2}}{6D_{r}},\end{eqnarray}$$

which is precisely what we would expect in the purely diffusive limit where both chemotactic alignment and indirect chemotactic sedimentation are weak. We recover isotropy of the active diffusivity tensor and the correct factor of $1/6$ for the dispersal of active swimmers in three dimensions (Berg Reference Berg1975).

Furthermore, we discover that the linear terms vanish and that up to quadratic order

(3.41) $$\begin{eqnarray}\displaystyle & \displaystyle D_{\bot }\sim \frac{U^{2}}{D_{r}}\left[\frac{1}{6}-\frac{\unicode[STIX]{x1D706}^{2}}{40}-\frac{\unicode[STIX]{x1D707}^{2}}{90}+\text{h.o.t.}\right], & \displaystyle\end{eqnarray}$$
(3.42) $$\begin{eqnarray}\displaystyle & \displaystyle D_{\Vert }\sim \frac{U^{2}}{D_{r}}\left[\frac{1}{6}-\frac{7\unicode[STIX]{x1D706}^{2}}{135}-\frac{\unicode[STIX]{x1D706}\unicode[STIX]{x1D707}}{9}-\frac{2\unicode[STIX]{x1D707}^{2}}{135}+\text{h.o.t.}\right]. & \displaystyle\end{eqnarray}$$

This explains our previous observation that the diffusivity has a local maximum at $\unicode[STIX]{x1D706}=0,\unicode[STIX]{x1D707}=0$ , which can be seen in figure 6. Thus, for small values of $\unicode[STIX]{x1D706}$ and $\unicode[STIX]{x1D707}$ both mechanisms of chemotactic alignment and indirect chemotactic sedimentation have the effect of suppressing diffusion relative to the mean drift velocity of the swimmers.

Figure 7. The ratio $D_{\bot }/D_{\Vert }$ computed using GTD theory in the $(\unicode[STIX]{x1D706},\unicode[STIX]{x1D707})$ plane. As expected, the contour $D_{\bot }/D_{\Vert }=1$ passes through the origin, because the system reverts to isotropic diffusion in the limit $\unicode[STIX]{x1D706},\unicode[STIX]{x1D707}\rightarrow 0$ . Increasing $\unicode[STIX]{x1D707}$ above this contour leads to a ratio $D_{\bot }/D_{\Vert }<1$ , meaning that indirect chemotactic sedimentation enhances diffusion along the chemical gradient relative to diffusion in the normal plane. On the other hand, increasing $\unicode[STIX]{x1D706}$ to the right of this contour leads to a ratio $D_{\bot }/D_{\Vert }>1$ , meaning that chemotactic alignment suppresses diffusion along the chemical gradient relative to diffusion in the normal plane.

4 Conclusions

The idea of using autophoresis to design particles capable of independent propulsion and reorientation in the presence of chemical stimuli has gained much attention in recent years, not least because it eliminates the need for active steering or fine control of external (e.g. electric or magnetic) fields. In this paper we derive a general law for the instantaneous behaviour of a spherical axisymmetric swimmer placed in a uniform chemical gradient, which extends and corrects results published in Saha et al. (Reference Saha, Golestanian and Ramaswamy2014). We also use our framework to calculate the linear and angular velocity of a Janus sphere, which has great relevance for experimental studies.

The main contribution of the present paper is to obtain a fully analytical solution for a non-trivial transport problem involving a chemically active phoretic swimmer placed in a uniform chemical gradient, a canonical set-up which had yet to be solved in a general form and presented in a pedagogical manner. Furthermore, our systematic analysis of the different sources contributing to chemotaxis could help inform the design of phoretic swimmers in future experiments, as we discussed in the subsection on Janus spheres.

We reinforce the rationale for our first modelling assumption by noting that the regime in which the surface reaction rate is linear in the substrate concentration is the most interesting and relevant case to artificial chemotaxis, because the product problem is sensitive to the substrate gradient in this limit. Our second modelling assumption of small Damköhler number could be relaxed, but in this case the substrate problem would no longer have a closed analytical solution and we would have to resort to numerical methods for solving a truncated system. The instantaneous model could be further extended through the inclusion of advective effects, which have been highlighted by several studies (Khair Reference Khair2013; Michelin, Lauga & Bartolo Reference Michelin, Lauga and Bartolo2013; Michelin & Lauga Reference Michelin and Lauga2014), or by considering the possible ionic effects recently pointed out by Brown & Poon (Reference Brown and Poon2014).

In the second half of the paper we present the calculations involved in applying GTD theory (Frankel & Brenner Reference Frankel and Brenner1991, Reference Frankel and Brenner1993) to the artificial chemotaxis of phoretic swimmers, and we obtain very good agreement between the continuum model and numerical simulations. We observe non-trivial variations of the active diffusivity in the two-dimensional parameter space of rotational Péclet number and indirect chemotactic index, which creates the possibility of novel pattern formation in systems where these parameters vary with position. The indirect chemotactic index is an added degree of complexity that is new to our continuum model compared to previous studies since it is a direct consequence of the complex chemical dynamics considered in our derivation of the instantaneous behaviour.

The limitations of the present continuum model must however be acknowledged. The main difference between the problem of autophoretic swimmers and that of gyrotactic micro-organisms is that the swimming velocity of the former depends on position, whereas for the latter it is constant. This imposes a further condition on the applicability of the model, namely that the swimming velocity is sufficiently slowly varying in space, which may not be the case in certain practical applications. Furthermore, our model is three-dimensional, whereas many experimental projects are concerned with two-dimensional distributions of autophoretic swimmers above a flat plate.

The model in this paper describes the long-time behaviour of a single phoretic swimmer, and therefore it is valid only for a population of non-interacting autophoretic swimmers. This makes our continuum model unable to explain the complex collective behaviour observed in dense populations of interacting phoretic swimmers (Palacci et al. Reference Palacci, Sacanna, Steinberg, Pine and Chaikin2013; Saha et al. Reference Saha, Golestanian and Ramaswamy2014). Nevertheless, it represents an appropriate starting point for investigating the emergent behaviour of dilute suspensions of chemically active colloidal particles, such as the recent work done on phoretic swimmers ‘riding’ active density waves (see Geiseler et al. Reference Geiseler, Hanggi, Marchesoni, Mulhern and Savel’ev2016).

In the phoretic swimming literature, many studies assume that active swimming leads to an enhanced, but still isotropic, spatial diffusion coefficient, whereas active reorientation leads to a modified rotational diffusion coefficient, and that the two mechanisms are independent of each other. We clearly show in this paper the importance of having an anisotropic diffusivity tensor which encapsulates the effect of both mechanisms, and which represents the correct way to think about the dispersion of autophoretic swimmers. On a phenomenological level, this represents an important contribution of our continuum model to the general understanding of the long-time behaviour of phoretic swimmers placed in a chemical gradient.

Acknowledgements

We thank the anonymous referees for their useful comments on the early version of the manuscript. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement 682754 to E.L.). This work was also funded by a Summer Research Studentship from Trinity College and a George and Lilian Schiff Studentship (M.T.-C.).

Appendix A

This appendix contains detailed calculations that complement §§ 2.4 and 2.5 from the derivation of the instantaneous behaviour.

A.1 Leading order: linear velocity due to substrate

To compute the linear velocity due to the leading-order substrate concentration, we must average the slip velocity

(A 1) $$\begin{eqnarray}\boldsymbol{v}_{slip}^{s,0}=\frac{3\unicode[STIX]{x1D707}_{s}(\unicode[STIX]{x1D703})}{2}[(-\unicode[STIX]{x1D735}s_{z}^{\infty }\sin \unicode[STIX]{x1D703}+\unicode[STIX]{x1D735}s_{y}^{\infty }\cos \unicode[STIX]{x1D703}\sin \unicode[STIX]{x1D719})\hat{\unicode[STIX]{x1D73D}}+\unicode[STIX]{x1D735}s_{y}^{\infty }\cos \unicode[STIX]{x1D719}\,\hat{\unicode[STIX]{x1D753}}]\end{eqnarray}$$

over the surface of the sphere. The contribution due to the $\hat{\unicode[STIX]{x1D73D}}$ term in (2.27) is given by

(A 2) $$\begin{eqnarray}\frac{1}{4\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}}\int _{0}^{\unicode[STIX]{x03C0}}\unicode[STIX]{x1D707}_{s}(\unicode[STIX]{x1D703})\frac{3}{2}(-\unicode[STIX]{x1D735}s_{z}^{\infty }\sin \unicode[STIX]{x1D703}+\unicode[STIX]{x1D735}s_{y}^{\infty }\cos \unicode[STIX]{x1D703}\sin \unicode[STIX]{x1D719})\left(\begin{array}{@{}c@{}}\cos \unicode[STIX]{x1D703}\cos \unicode[STIX]{x1D719}\\ \cos \unicode[STIX]{x1D703}\sin \unicode[STIX]{x1D719}\\ -\text{sin}\unicode[STIX]{x1D703}\\ \end{array}\right)\sin \unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D719},\end{eqnarray}$$

where we have expanded the vector $\hat{\unicode[STIX]{x1D73D}}$ into Cartesian coordinates. The only terms with a non-zero average over azimuthal angle are

(A 3) $$\begin{eqnarray}\frac{3}{8\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}}\int _{0}^{\unicode[STIX]{x03C0}}\mathop{\sum }_{l=0}^{\infty }\unicode[STIX]{x1D707}_{sl}P_{l}(\cos \unicode[STIX]{x1D703})(\unicode[STIX]{x1D735}s_{z}^{\infty }\sin ^{2}\unicode[STIX]{x1D703}\hat{\boldsymbol{z}}+\unicode[STIX]{x1D735}s_{y}^{\infty }\cos ^{2}\unicode[STIX]{x1D703}\sin ^{2}\unicode[STIX]{x1D719}\hat{\boldsymbol{y}})\sin \unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D719},\end{eqnarray}$$

where we have replaced the mobility $\unicode[STIX]{x1D707}_{s}(\unicode[STIX]{x1D703})$ by its series representation. Next, we carry out the integral over the azimuthal angle and use a standard change of variable $u=\cos \unicode[STIX]{x1D703}$ for the polar angle to obtain

(A 4) $$\begin{eqnarray}\frac{3}{4}\mathop{\sum }_{l=0}^{\infty }\unicode[STIX]{x1D707}_{sl}\int _{-1}^{+1}P_{l}(u)\left(\unicode[STIX]{x1D735}s_{z}^{\infty }(1-u^{2})\hat{\boldsymbol{z}}+\frac{1}{2}\unicode[STIX]{x1D735}s_{y}^{\infty }u^{2}\hat{\boldsymbol{y}}\right)\,\text{d}u.\end{eqnarray}$$

Using the classical identities (B 2) and (B 5) we obtain the final result for the $\hat{\unicode[STIX]{x1D73D}}$ contribution as

(A 5) $$\begin{eqnarray}\left(\unicode[STIX]{x1D707}_{s0}-{\textstyle \frac{1}{5}}\unicode[STIX]{x1D707}_{s2}\right)\unicode[STIX]{x1D735}s_{z}^{\infty }\hat{\boldsymbol{z}}+\left({\textstyle \frac{1}{4}}\unicode[STIX]{x1D707}_{s0}+{\textstyle \frac{1}{10}}\unicode[STIX]{x1D707}_{s2}\right)\unicode[STIX]{x1D735}s_{y}^{\infty }\hat{\boldsymbol{y}}.\end{eqnarray}$$

Following the same principles, the contribution due to the $\hat{\unicode[STIX]{x1D753}}$ term in (2.27) is

(A 6) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{1}{4\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}}\int _{0}^{\unicode[STIX]{x03C0}}\,\unicode[STIX]{x1D707}_{s}(\unicode[STIX]{x1D703})\frac{3}{2}\unicode[STIX]{x1D735}s_{y}^{\infty }\cos \unicode[STIX]{x1D719}\left(\begin{array}{@{}c@{}}-\text{sin}\unicode[STIX]{x1D719}\\ \cos \unicode[STIX]{x1D719}\\ 0\\ \end{array}\right)\sin \unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D719}, & \displaystyle\end{eqnarray}$$
(A 7) $$\begin{eqnarray}\displaystyle & \displaystyle =\frac{3}{8\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}}\int _{0}^{\unicode[STIX]{x03C0}}\mathop{\sum }_{l=0}^{\infty }\unicode[STIX]{x1D707}_{sl}P_{l}(\cos \unicode[STIX]{x1D703})\unicode[STIX]{x1D735}s_{y}^{\infty }\cos ^{2}\unicode[STIX]{x1D719}\,\hat{\boldsymbol{y}}\sin \unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D719}, & \displaystyle\end{eqnarray}$$
(A 8) $$\begin{eqnarray}\displaystyle & \displaystyle =\frac{3}{8}\mathop{\sum }_{l=0}^{\infty }\unicode[STIX]{x1D707}_{sl}\int _{-1}^{+1}P_{l}(u)\unicode[STIX]{x1D735}s_{y}^{\infty }\hat{\boldsymbol{y}}\,\text{d}u, & \displaystyle\end{eqnarray}$$
(A 9) $$\begin{eqnarray}\displaystyle & \displaystyle =\frac{3}{4}\unicode[STIX]{x1D707}_{s0}\unicode[STIX]{x1D735}s_{y}^{\infty }\hat{\boldsymbol{y}}. & \displaystyle\end{eqnarray}$$

As a consequence, the slip velocity averaged over the surface is finally

(A 10) $$\begin{eqnarray}\langle \boldsymbol{v}_{slip}^{s,0}\rangle =\left(\unicode[STIX]{x1D707}_{s0}-{\textstyle \frac{1}{5}}\unicode[STIX]{x1D707}_{s2}\right)\unicode[STIX]{x1D735}s_{z}^{\infty }\hat{\boldsymbol{z}}+\left(\unicode[STIX]{x1D707}_{s0}+{\textstyle \frac{1}{10}}\unicode[STIX]{x1D707}_{s2}\right)\unicode[STIX]{x1D735}s_{y}^{\infty }\hat{\boldsymbol{y}}.\end{eqnarray}$$

From (2.10) we deduce that the contribution of the substrate to the linear velocity of the particle, at leading order, is

(A 11) $$\begin{eqnarray}\boldsymbol{V}_{0}^{s}=-\left(\unicode[STIX]{x1D707}_{s0}+{\textstyle \frac{1}{10}}\unicode[STIX]{x1D707}_{s2}\right)\unicode[STIX]{x1D735}s^{\infty }+{\textstyle \frac{3}{10}}\unicode[STIX]{x1D707}_{s2}\,\hat{\boldsymbol{z}}\hat{\boldsymbol{z}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}s^{\infty }.\end{eqnarray}$$

A.2 Leading order: angular velocity due to substrate

In order to calculate the angular velocity, we must average the following quantity

(A 12) $$\begin{eqnarray}\hat{\boldsymbol{r}}\wedge \boldsymbol{v}_{slip}^{s,0}=\frac{3\unicode[STIX]{x1D707}_{s}(\unicode[STIX]{x1D703})}{2}[(\unicode[STIX]{x1D735}s_{y}^{\infty }\cos \unicode[STIX]{x1D703}\sin \unicode[STIX]{x1D719}-\unicode[STIX]{x1D735}s_{z}^{\infty }\sin \unicode[STIX]{x1D703})\hat{\unicode[STIX]{x1D753}}-\unicode[STIX]{x1D735}s_{y}^{\infty }\cos \unicode[STIX]{x1D719}\,\hat{\unicode[STIX]{x1D73D}}]\end{eqnarray}$$

over the surface of the sphere. The contribution due to the $\hat{\unicode[STIX]{x1D753}}$ term in (A 12) is

(A 13) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{1}{4\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}}\int _{0}^{\unicode[STIX]{x03C0}}\unicode[STIX]{x1D707}_{s}(\unicode[STIX]{x1D703})\frac{3}{2}(\unicode[STIX]{x1D735}s_{y}^{\infty }\cos \unicode[STIX]{x1D703}\sin \unicode[STIX]{x1D719}-\unicode[STIX]{x1D735}s_{z}^{\infty }\sin \unicode[STIX]{x1D703})\left(\begin{array}{@{}c@{}}-\text{sin}\unicode[STIX]{x1D719}\\ \cos \unicode[STIX]{x1D719}\\ 0\\ \end{array}\right)\sin \unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D719}, & \displaystyle\end{eqnarray}$$
(A 14) $$\begin{eqnarray}\displaystyle & \displaystyle =\frac{3}{8\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}}\int _{0}^{\unicode[STIX]{x03C0}}\mathop{\sum }_{l=0}^{\infty }\unicode[STIX]{x1D707}_{sl}P_{l}(\cos \unicode[STIX]{x1D703})(-\unicode[STIX]{x1D735}s_{y}^{\infty }\cos \unicode[STIX]{x1D703}\sin ^{2}\unicode[STIX]{x1D719}\hat{\boldsymbol{x}})\sin \unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D719}, & \displaystyle\end{eqnarray}$$
(A 15) $$\begin{eqnarray}\displaystyle & \displaystyle =-\frac{3}{8}\mathop{\sum }_{l=0}^{\infty }\unicode[STIX]{x1D707}_{sl}\int _{-1}^{+1}P_{l}(u)\unicode[STIX]{x1D735}s_{y}^{\infty }u\,\hat{\boldsymbol{x}}\,\text{d}u, & \displaystyle\end{eqnarray}$$
(A 16) $$\begin{eqnarray}\displaystyle & \displaystyle =-\frac{1}{4}\unicode[STIX]{x1D707}_{s1}\unicode[STIX]{x1D735}s_{y}^{\infty }\hat{\boldsymbol{x}}. & \displaystyle\end{eqnarray}$$

Next, we evaluate the contribution due to the $\hat{\unicode[STIX]{x1D73D}}$ term in (A 12) as

(A 17) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{1}{4\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}}\int _{0}^{\unicode[STIX]{x03C0}}\unicode[STIX]{x1D707}_{s}(\unicode[STIX]{x1D703})\frac{3}{2}(-\unicode[STIX]{x1D735}s_{y}^{\infty }\cos \unicode[STIX]{x1D719})\left(\begin{array}{@{}c@{}}\cos \unicode[STIX]{x1D703}\cos \unicode[STIX]{x1D719}\\ \cos \unicode[STIX]{x1D703}\sin \unicode[STIX]{x1D719}\\ -\text{sin}\unicode[STIX]{x1D703}\\ \end{array}\right)\sin \unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D719}, & \displaystyle\end{eqnarray}$$
(A 18) $$\begin{eqnarray}\displaystyle & \displaystyle =\frac{3}{8\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}}\int _{0}^{\unicode[STIX]{x03C0}}\mathop{\sum }_{l=0}^{\infty }\unicode[STIX]{x1D707}_{sl}P_{l}(\cos \unicode[STIX]{x1D703})(-\unicode[STIX]{x1D735}s_{y}^{\infty }\cos \unicode[STIX]{x1D703}\cos ^{2}\unicode[STIX]{x1D719}\hat{\boldsymbol{x}})\sin \unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D719}, & \displaystyle\end{eqnarray}$$
(A 19) $$\begin{eqnarray}\displaystyle & \displaystyle =-\frac{3}{8}\mathop{\sum }_{l=0}^{\infty }\unicode[STIX]{x1D707}_{sl}\int _{-1}^{+1}P_{l}(u)\unicode[STIX]{x1D735}s_{y}^{\infty }u\,\hat{\boldsymbol{x}}\,\text{d}u, & \displaystyle\end{eqnarray}$$
(A 20) $$\begin{eqnarray}\displaystyle & \displaystyle =-\frac{1}{4}\unicode[STIX]{x1D707}_{s1}\unicode[STIX]{x1D735}s_{y}^{\infty }\hat{\boldsymbol{x}}. & \displaystyle\end{eqnarray}$$

Adding these two results gives us

(A 21) $$\begin{eqnarray}\langle \hat{\boldsymbol{r}}\wedge \boldsymbol{v}_{slip}^{s,0}\rangle =-{\textstyle \frac{1}{2}}\unicode[STIX]{x1D707}_{s1}\unicode[STIX]{x1D735}s_{y}^{\infty }\hat{\boldsymbol{x}},\end{eqnarray}$$

and from (2.11) we deduce that the angular velocity due to the substrate is

(A 22) $$\begin{eqnarray}\unicode[STIX]{x1D74E}_{0}^{s}=-\frac{3\unicode[STIX]{x1D707}_{s1}}{4R}\,\hat{\boldsymbol{z}}\wedge \unicode[STIX]{x1D735}s^{\infty }.\end{eqnarray}$$

A.3 First order: substrate concentration

The normal gradient of the substrate concentration from (2.31) is

(A 23) $$\begin{eqnarray}\left.\frac{\unicode[STIX]{x2202}s_{1}}{\unicode[STIX]{x2202}n}\right|_{r=R}=\mathop{\sum }_{l=0}^{\infty }\frac{-(l+1)}{R}[A_{l}P_{l}(\cos \unicode[STIX]{x1D703})+B_{l}P_{l}^{1}(\cos \unicode[STIX]{x1D703})\sin \unicode[STIX]{x1D719}].\end{eqnarray}$$

We can substitute this result, together with the series representation of $\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D703})$ and the expression for $s_{0}$ from (2.25), into the normal-flux boundary condition from (2.19b ) to get

(A 24) $$\begin{eqnarray}\displaystyle & & \displaystyle -\frac{D_{s}}{R}\mathop{\sum }_{l=0}^{\infty }(l+1)\left[A_{l}P_{l}(\cos \unicode[STIX]{x1D703})+B_{l}P_{l}^{1}(\cos \unicode[STIX]{x1D703})\sin \unicode[STIX]{x1D719}\right]\nonumber\\ \displaystyle & & \displaystyle \quad =\unicode[STIX]{x1D705}_{1}\left(s_{b}(\boldsymbol{r}_{c})+\frac{3}{2}\unicode[STIX]{x1D735}s_{z}^{\infty }R\cos \unicode[STIX]{x1D703}+\frac{3}{2}\unicode[STIX]{x1D735}s_{y}^{\infty }R\sin \unicode[STIX]{x1D703}\sin \unicode[STIX]{x1D719}\right)\mathop{\sum }_{j=0}^{\infty }\unicode[STIX]{x1D70E}_{j}P_{j}(\cos \unicode[STIX]{x1D703}).\qquad\end{eqnarray}$$

We next separate the components with no $\unicode[STIX]{x1D719}$ -dependence by averaging over the azimuthal angle and then take $\int _{0}^{\unicode[STIX]{x03C0}}\ldots P_{k}(\cos \unicode[STIX]{x1D703})\sin \unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D703}$ of both sides to get

(A 25) $$\begin{eqnarray}\displaystyle & & \displaystyle -\frac{D_{s}}{R}\mathop{\sum }_{l=0}^{\infty }(l+1)A_{l}\int _{-1}^{+1}P_{l}(u)P_{k}(u)\,\text{d}u\nonumber\\ \displaystyle & & \displaystyle \quad =\unicode[STIX]{x1D705}_{1}\mathop{\sum }_{j=0}^{\infty }\unicode[STIX]{x1D70E}_{j}\left(s_{b}(\boldsymbol{r}_{c})\int _{-1}^{+1}P_{j}(u)P_{k}(u)\,\text{d}u+\frac{3}{2}\unicode[STIX]{x1D735}s_{z}^{\infty }R\int _{-1}^{+1}uP_{j}(u)P_{k}(u)\,\text{d}u\right).\qquad\end{eqnarray}$$

Using identities (B 2) and (B 6), we obtain the $A_{k}$ set of coefficients as

(A 26) $$\begin{eqnarray}A_{k}=-\frac{\unicode[STIX]{x1D705}_{1}R}{(k+1)D_{s}}\left[s_{b}(\boldsymbol{r}_{c})\unicode[STIX]{x1D70E}_{k}+\frac{3}{2}\unicode[STIX]{x1D735}s_{z}^{\infty }R\left(\frac{k+1}{2k+3}\unicode[STIX]{x1D70E}_{k+1}+\frac{k}{2k-1}\unicode[STIX]{x1D70E}_{k-1}\right)\right].\end{eqnarray}$$

We can also separate the $\sin \unicode[STIX]{x1D719}$ components if we multiply (A 24) by $\sin \unicode[STIX]{x1D719}$ and average over the azimuthal angle. Then we exploit the property that $P_{l}^{1}(\cos \unicode[STIX]{x1D703})=-\sin \unicode[STIX]{x1D703}~P_{l}^{\prime }(\cos \unicode[STIX]{x1D703})$ and change variables to $u=\cos \unicode[STIX]{x1D703}$ to obtain

(A 27) $$\begin{eqnarray}\frac{D_{s}}{R}\mathop{\sum }_{l=0}^{\infty }(l+1)B_{l}P_{l}^{\prime }(u)=\frac{3}{2}\unicode[STIX]{x1D705}_{1}\unicode[STIX]{x1D735}s_{y}^{\infty }R\mathop{\sum }_{j=0}^{\infty }\unicode[STIX]{x1D70E}_{j}P_{j}(u).\end{eqnarray}$$

Taking $\int _{-1}^{+1}\ldots (1-u^{2})P_{k}^{\prime }(u)\,\text{d}u$ of both sides of this equation, we obtain

(A 28) $$\begin{eqnarray}\displaystyle & & \displaystyle \frac{D_{s}}{R}\mathop{\sum }_{l=0}^{\infty }(l+1)B_{l}\int _{-1}^{+1}(1-u^{2})P_{l}^{\prime }(u)P_{k}^{\prime }(u)\,\text{d}u\nonumber\\ \displaystyle & & \displaystyle \quad =\frac{3}{2}\unicode[STIX]{x1D705}_{1}\unicode[STIX]{x1D735}s_{y}^{\infty }R\mathop{\sum }_{j=0}^{\infty }\unicode[STIX]{x1D70E}_{j}\int _{-1}^{+1}(1-u^{2})P_{j}(u)P_{k}^{\prime }(u)\,\text{d}u.\end{eqnarray}$$

Finally, using the identities (B 7) and (B 8) we obtain the $B_{k}$ set of coefficients as

(A 29) $$\begin{eqnarray}B_{k}=-\frac{3\unicode[STIX]{x1D705}_{1}R^{2}\unicode[STIX]{x1D735}s_{y}^{\infty }}{2(k+1)D_{s}}\left(\frac{\unicode[STIX]{x1D70E}_{k+1}}{2k+3}-\frac{\unicode[STIX]{x1D70E}_{k-1}}{2k-1}\right).\end{eqnarray}$$

Although the first-order product distribution is now fully determined (since both the $A_{k}$ and $B_{k}$ are known), for the sake of simplicity we will continue the calculations below using the expression

(A 30) $$\begin{eqnarray}s_{1}=\mathop{\sum }_{l=0}^{\infty }\left(\frac{R}{r}\right)^{l+1}[A_{l}P_{l}(\cos \unicode[STIX]{x1D703})-B_{l}\sin \unicode[STIX]{x1D703}P_{l}^{\prime }(\cos \unicode[STIX]{x1D703})\sin \unicode[STIX]{x1D719}]\end{eqnarray}$$

and will substitute the values for $A_{l}$ and $B_{l}$ only when necessary.

A.4 First order: linear velocity due to substrate

To calculate the first-order contribution of the substrate problem to the linear velocity, we must average the slip velocity

(A 31) $$\begin{eqnarray}\displaystyle \boldsymbol{v}_{slip}^{s,1} & = & \displaystyle \frac{\unicode[STIX]{x1D707}_{s}(\unicode[STIX]{x1D703})}{R}\mathop{\sum }_{l=0}^{\infty }\left(-A_{l}\sin \unicode[STIX]{x1D703}P_{l}^{\prime }(\cos \unicode[STIX]{x1D703})-B_{l}\cos \unicode[STIX]{x1D703}P_{l}^{\prime }(\cos \unicode[STIX]{x1D703})\sin \unicode[STIX]{x1D719}\right.\nonumber\\ \displaystyle & & \displaystyle +\,\left.B_{l}\sin ^{2}\unicode[STIX]{x1D703}P_{l}^{\prime \prime }(\cos \unicode[STIX]{x1D703})\sin \unicode[STIX]{x1D719}\right)\hat{\unicode[STIX]{x1D73D}}+\frac{\unicode[STIX]{x1D707}_{s}(\unicode[STIX]{x1D703})}{R}\mathop{\sum }_{l=0}^{\infty }\left(-B_{l}P_{l}^{\prime }(\cos \unicode[STIX]{x1D703})\cos \unicode[STIX]{x1D719}\right)\hat{\unicode[STIX]{x1D753}}\qquad \quad\end{eqnarray}$$

over the surface of the swimmer.

The contribution to $\langle \boldsymbol{v}_{slip}^{s,1}\rangle$ due to the $\hat{\unicode[STIX]{x1D73D}}$ terms in (A 31) is given by

(A 32) $$\begin{eqnarray}\displaystyle & & \displaystyle \frac{1}{4\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}}\int _{0}^{\unicode[STIX]{x03C0}}\unicode[STIX]{x1D707}_{s}(\unicode[STIX]{x1D703})\frac{1}{R}\mathop{\sum }_{l=0}^{\infty }\left[-A_{l}\sin \unicode[STIX]{x1D703}P_{l}^{\prime }(\cos \unicode[STIX]{x1D703})-B_{l}\cos \unicode[STIX]{x1D703}P_{l}^{\prime }(\cos \unicode[STIX]{x1D703})\sin \unicode[STIX]{x1D719}\right.\nonumber\\ \displaystyle & & \displaystyle \quad \left.+\,B_{l}\sin ^{2}\unicode[STIX]{x1D703}P_{l}^{\prime \prime }(\cos \unicode[STIX]{x1D703})\sin \unicode[STIX]{x1D719}\right]\left(\begin{array}{@{}c@{}}\cos \unicode[STIX]{x1D703}\cos \unicode[STIX]{x1D719}\\ \cos \unicode[STIX]{x1D703}\sin \unicode[STIX]{x1D719}\\ -\text{sin}\unicode[STIX]{x1D703}\\ \end{array}\right)\sin \unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D719}.\end{eqnarray}$$

We can first average over the azimuthal angle and perform the usual change of variable for the polar angle in order to obtain the average

(A 33) $$\begin{eqnarray}\displaystyle & & \displaystyle \frac{1}{2R}\mathop{\sum }_{m,l=0}^{\infty }\unicode[STIX]{x1D707}_{sm}A_{l}\int _{-1}^{+1}(1-u^{2})P_{m}(u)P_{l}^{\prime }(u)\,\text{d}u\,\hat{\boldsymbol{z}}\nonumber\\ \displaystyle & & \displaystyle \quad +\,\frac{1}{4R}\mathop{\sum }_{m,l=0}^{\infty }\unicode[STIX]{x1D707}_{sm}B_{l}\left[\int _{-1}^{+1}u(1-u^{2})P_{m}(u)P_{l}^{\prime \prime }(u)\,\text{d}u-\int _{-1}^{+1}u^{2}P_{m}(u)P_{l}^{\prime }(u)\,\text{d}u\right]\hat{\boldsymbol{y}},\qquad \quad\end{eqnarray}$$

where we have replaced $\unicode[STIX]{x1D707}_{s}(\unicode[STIX]{x1D703})$ by its series representation.

The $\hat{\boldsymbol{z}}$ component of the integral in (A 33) is straightforward, but for the $\hat{\boldsymbol{y}}$ component we need to do an integration by parts on the first integral to bring it to the form

(A 34) $$\begin{eqnarray}-\int _{-1}^{+1}\left[u(1-u^{2})~P_{m}^{\prime }(u)+(1-3u^{2})P_{m}(u)\right]P_{l}^{\prime }(u)\,\text{d}u.\end{eqnarray}$$

As a result, the contribution to linear velocity due to the $\hat{\unicode[STIX]{x1D73D}}$ terms in (A 31) is given by

(A 35) $$\begin{eqnarray}\displaystyle & & \displaystyle \frac{1}{2R}\mathop{\sum }_{m,l=0}^{\infty }\unicode[STIX]{x1D707}_{sm}A_{l}\int _{-1}^{+1}(1-u^{2})P_{m}(u)P_{l}^{\prime }(u)\,\text{d}u~\hat{\boldsymbol{z}}\nonumber\\ \displaystyle & & \displaystyle \quad +\,\frac{1}{4R}\mathop{\sum }_{m,l=0}^{\infty }\unicode[STIX]{x1D707}_{sm}B_{l}\left[-\int _{-1}^{+1}u(1-u^{2})P_{m}^{\prime }(u)P_{l}^{\prime }(u)\,\text{d}u-2\int _{-1}^{+1}(1-u^{2})P_{m}(u)P_{l}^{\prime }(u)\,\text{d}u\right.\nonumber\\ \displaystyle & & \displaystyle \quad \left.+\,\int _{-1}^{+1}P_{m}(u)P_{l}^{\prime }(u)\,\text{d}u\right]\hat{\boldsymbol{y}}.\end{eqnarray}$$

We can evaluate explicitly all but the last of those integrals, but fortunately the contribution to $\langle \boldsymbol{v}_{slip}^{s,1}\rangle$ due to the $\hat{\unicode[STIX]{x1D753}}$ term in (A 31) is exactly

(A 36) $$\begin{eqnarray}\displaystyle & & \displaystyle \frac{1}{4\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}}\int _{0}^{\unicode[STIX]{x03C0}}\unicode[STIX]{x1D707}_{s}(\unicode[STIX]{x1D703})\frac{1}{R}\mathop{\sum }_{l=0}^{\infty }(-B_{l}P_{l}^{\prime }(\cos \unicode[STIX]{x1D703})\cos \unicode[STIX]{x1D719})\left(\begin{array}{@{}c@{}}-\text{sin}\unicode[STIX]{x1D719}\\ \cos \unicode[STIX]{x1D719}\\ 0\\ \end{array}\right)\sin \unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D719}\nonumber\\ \displaystyle & & \displaystyle \quad =-\frac{1}{4R}\mathop{\sum }_{m,l=0}^{\infty }\unicode[STIX]{x1D707}_{sm}B_{l}\int _{-1}^{+1}P_{m}(u)P_{l}^{\prime }(u)\,\text{d}u\,\hat{\boldsymbol{y}},\end{eqnarray}$$

so the undetermined integrals cancel out exactly. Using identities (B 8), (B 9) and (2.10), we deduce that the first-order substrate contribution to linear velocity is given by

(A 37) $$\begin{eqnarray}\displaystyle \boldsymbol{V}_{1}^{s} & = & \displaystyle \frac{1}{R}\mathop{\sum }_{l=0}^{\infty }\frac{l(l+1)}{2l+1}A_{l}\left(\frac{\unicode[STIX]{x1D707}_{s,l+1}}{2l+3}-\frac{\unicode[STIX]{x1D707}_{s,l-1}}{2l-1}\right)\hat{\boldsymbol{z}}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{1}{2R}\mathop{\sum }_{l=0}^{\infty }\frac{l(l+1)}{2l+1}B_{l}\left(\frac{l\unicode[STIX]{x1D707}_{s,l+1}}{2l+3}+\frac{(l+1)\unicode[STIX]{x1D707}_{s,l-1}}{2l-1}\right)\hat{\boldsymbol{y}}.\end{eqnarray}$$

If we substitute the values of the coefficients $A_{l}$ and $B_{l}$ from (A 26) and (A 29), we get

(A 38) $$\begin{eqnarray}\displaystyle \boldsymbol{V}_{1}^{s} & = & \displaystyle -\frac{\unicode[STIX]{x1D705}_{1}s_{b}(\boldsymbol{r}_{c})}{D_{s}}\mathop{\sum }_{l=0}^{\infty }\left(\frac{l}{2l+1}\right)~\unicode[STIX]{x1D70E}_{l}\left(\frac{\unicode[STIX]{x1D707}_{s,l+1}}{2l+3}-\frac{\unicode[STIX]{x1D707}_{s,l-1}}{2l-1}\right)\hat{\boldsymbol{z}}\nonumber\\ \displaystyle & & \displaystyle -\,\frac{3\unicode[STIX]{x1D705}_{1}R}{2D_{s}}\mathop{\sum }_{l=0}^{\infty }\left(\frac{l}{2l+1}\right)\left(\frac{l+1}{2l+3}\unicode[STIX]{x1D70E}_{l+1}+\frac{l}{2l-1}\unicode[STIX]{x1D70E}_{l-1}\right)\left(\frac{\unicode[STIX]{x1D707}_{s,l+1}}{2l+3}-\frac{\unicode[STIX]{x1D707}_{s,l-1}}{2l-1}\right)\unicode[STIX]{x1D735}s_{z}^{\infty }\hat{\boldsymbol{z}}\nonumber\\ \displaystyle & & \displaystyle -\,\frac{3\unicode[STIX]{x1D705}_{1}R}{4D_{s}}\mathop{\sum }_{l=0}^{\infty }\left(\frac{l}{2l+1}\right)\left(\frac{\unicode[STIX]{x1D70E}_{l+1}}{2l+3}-\frac{\unicode[STIX]{x1D70E}_{l-1}}{2l-1}\right)\left(\frac{l+1}{2l-1}\unicode[STIX]{x1D707}_{s,l-1}+\frac{l}{2l+3}\unicode[STIX]{x1D707}_{s,l+1}\right)\unicode[STIX]{x1D735}s_{y}^{\infty }\hat{\boldsymbol{y}}.\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

We can further manipulate the last two terms in order to get an expression in terms of $\unicode[STIX]{x1D735}s^{\infty }$ and $\hat{\boldsymbol{z}}\hat{\boldsymbol{z}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}s^{\infty }$ for the linear velocity as

(A 39) $$\begin{eqnarray}\displaystyle \boldsymbol{V}_{1}^{s} & = & \displaystyle -\frac{\unicode[STIX]{x1D705}_{1}s_{b}(\boldsymbol{r}_{c})}{D_{s}}\mathop{\sum }_{l=1}^{\infty }\left(\frac{l}{2l+1}\right)\unicode[STIX]{x1D70E}_{l}\left(\frac{\unicode[STIX]{x1D707}_{s,l+1}}{2l+3}-\frac{\unicode[STIX]{x1D707}_{s,l-1}}{2l-1}\right)\hat{\boldsymbol{z}}\nonumber\\ \displaystyle & & \displaystyle -\,\frac{3\unicode[STIX]{x1D705}_{1}R}{4D_{s}}\mathop{\sum }_{l=1}^{\infty }\left(\frac{l}{2l+1}\right)\left(\frac{\unicode[STIX]{x1D70E}_{l+1}}{2l+3}-\frac{\unicode[STIX]{x1D70E}_{l-1}}{2l-1}\right)\left(\frac{l+1}{2l-1}\unicode[STIX]{x1D707}_{s,l-1}+\frac{l}{2l+3}\unicode[STIX]{x1D707}_{s,l+1}\right)\unicode[STIX]{x1D735}s^{\infty }\nonumber\\ \displaystyle & & \displaystyle +\,\frac{3\unicode[STIX]{x1D705}_{1}R}{4D_{s}}\mathop{\sum }_{l=1}^{\infty }\left(\frac{l}{2l+1}\right)\left(\frac{3(l+1)\unicode[STIX]{x1D70E}_{l+1}\unicode[STIX]{x1D707}_{s,l-1}}{(2l+3)(2l-1)}-\frac{(l+2)\unicode[STIX]{x1D70E}_{l+1}\unicode[STIX]{x1D707}_{s,l+1}}{(2l+3)^{2}}\right.\nonumber\\ \displaystyle & & \displaystyle +\,\left.\frac{(l-1)\unicode[STIX]{x1D70E}_{l-1}\unicode[STIX]{x1D707}_{s,l-1}}{(2l-1)^{2}}-\frac{3l\unicode[STIX]{x1D70E}_{l-1}\unicode[STIX]{x1D707}_{s,l+1}}{(2l-1)(2l+3)}\right)\hat{\boldsymbol{z}}\hat{\boldsymbol{z}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}s^{\infty }.\end{eqnarray}$$

A.5 First order: angular velocity due to substrate

We next need to evaluate $\langle \hat{\boldsymbol{r}}\wedge \boldsymbol{v}_{slip}^{s,1}\rangle$ for the angular velocity. This angular integral is equal to

(A 40) $$\begin{eqnarray}\displaystyle & & \displaystyle \frac{1}{4\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}}\int _{0}^{\unicode[STIX]{x03C0}}\unicode[STIX]{x1D707}_{s}(\unicode[STIX]{x1D703})\left[\frac{1}{R}\mathop{\sum }_{l=0}^{\infty }\left(-A_{l}\sin \unicode[STIX]{x1D703}P_{l}^{\prime }(\cos \unicode[STIX]{x1D703})-B_{l}\cos \unicode[STIX]{x1D703}P_{l}^{\prime }(\cos \unicode[STIX]{x1D703})\sin \unicode[STIX]{x1D719}\right.\right.\nonumber\\ \displaystyle & & \displaystyle \quad +\,\left.\left.B_{l}\sin ^{2}\unicode[STIX]{x1D703}P_{l}^{\prime \prime }(\cos \unicode[STIX]{x1D703})\sin \unicode[STIX]{x1D719}\right)\hat{\unicode[STIX]{x1D753}}+\frac{1}{R}\mathop{\sum }_{l=0}^{\infty }B_{l}P_{l}^{\prime }(\cos \unicode[STIX]{x1D703})\cos \unicode[STIX]{x1D719}\,\hat{\unicode[STIX]{x1D73D}}\right]\sin \unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D703}\,\text{d}\unicode[STIX]{x1D719}.\qquad \quad\end{eqnarray}$$

The $\hat{\unicode[STIX]{x1D753}}$ contribution of this integral is

(A 41) $$\begin{eqnarray}\frac{1}{4R}\mathop{\sum }_{m,l=0}^{\infty }\unicode[STIX]{x1D707}_{sm}B_{l}\left[\int _{-1}^{+1}uP_{m}(u)P_{l}^{\prime }(u)\,\text{d}u-\int _{-1}^{+1}(1-u^{2})P_{m}(u)P_{l}^{\prime \prime }(u)\,\text{d}u\right]\hat{\boldsymbol{x}}.\end{eqnarray}$$

We may integrate by parts the second integral in order to get

(A 42) $$\begin{eqnarray}\frac{1}{4R}\mathop{\sum }_{m,l=0}^{\infty }\unicode[STIX]{x1D707}_{sm}B_{l}\left[-\int _{-1}^{+1}uP_{m}(u)P_{l}^{\prime }(u)\,\text{d}u+\int _{-1}^{+1}(1-u^{2})P_{m}^{\prime }(u)P_{l}^{\prime }(u)\,\text{d}u\right]\hat{\boldsymbol{x}}.\end{eqnarray}$$

The $\hat{\unicode[STIX]{x1D73D}}$ contribution to the integral in (A 40) is simply given by

(A 43) $$\begin{eqnarray}\frac{1}{4R}\mathop{\sum }_{m,l=0}^{\infty }\unicode[STIX]{x1D707}_{sm}B_{l}\int _{-1}^{+1}uP_{m}(u)P_{l}^{\prime }(u)\,\text{d}u\,\hat{\boldsymbol{x}},\end{eqnarray}$$

so by adding the two angular contributions and using identity (B 7) we obtain the surface average

(A 44) $$\begin{eqnarray}\langle \hat{\boldsymbol{r}}\wedge \boldsymbol{v}_{slip}^{s,1}\rangle =\frac{1}{2R}\mathop{\sum }_{l=0}^{\infty }\frac{l(l+1)}{2l+1}\unicode[STIX]{x1D707}_{sl}B_{l}\,\hat{\boldsymbol{x}}.\end{eqnarray}$$

Given (2.11), the angular velocity due to the first-order substrate concentration takes the final form

(A 45) $$\begin{eqnarray}\unicode[STIX]{x1D74E}_{1}^{s}=-\frac{9\unicode[STIX]{x1D705}_{1}}{8D_{s}}\mathop{\sum }_{l=1}^{\infty }\left(\frac{l}{2l+1}\right)\unicode[STIX]{x1D707}_{sl}\left(\frac{\unicode[STIX]{x1D70E}_{l+1}}{2l+3}-\frac{\unicode[STIX]{x1D70E}_{l-1}}{2l-1}\right)\,\hat{\boldsymbol{z}}\wedge \unicode[STIX]{x1D735}s^{\infty },\end{eqnarray}$$

where we have substituted the value of $B_{l}$ from (A 29).

Appendix B

This appendix contains a list of all the properties and integral identities of Legendre polynomials that we have used in the rest of the paper. These are:

  1. (i) The differential equation satisfied by Legendre polynomials,

    (B 1) $$\begin{eqnarray}((1-u^{2})P_{k}^{\prime }(u))^{\prime }+k(k+1)P_{k}(u)=0.\end{eqnarray}$$
  2. (ii) The orthogonality condition for Legendre polynomials,

    (B 2) $$\begin{eqnarray}\int _{-1}^{+1}P_{k}(u)P_{l}(u)\,\text{d}u=\frac{2}{2k+1}\unicode[STIX]{x1D6FF}_{lk}.\end{eqnarray}$$
  3. (iii) Two useful recurrence relations involving Legendre polynomials and their derivatives,

    (B 3) $$\begin{eqnarray}\displaystyle & \displaystyle (2k+1)uP_{k}=(k+1)P_{k+1}+kP_{k-1}, & \displaystyle\end{eqnarray}$$
    (B 4) $$\begin{eqnarray}\displaystyle & \displaystyle (1-u^{2})P_{k}^{\prime }=k\left(P_{k-1}-uP_{k}\right). & \displaystyle\end{eqnarray}$$
  4. (iv) The first three Legendre polynomials,

    (B 5a-c ) $$\begin{eqnarray}P_{0}(u)=1,\quad P_{1}(u)=u,\quad P_{2}(u)={\textstyle \frac{1}{2}}(3u^{2}-1).\end{eqnarray}$$
  5. (v) A series of integral identities which can be derived from the basic properties (B 1)–(B 4):

    (B 6) $$\begin{eqnarray}\displaystyle & \displaystyle \int _{-1}^{+1}uP_{k}(u)P_{l}(u)\,\text{d}u=\frac{2}{2k+1}\left(\frac{k+1}{2k+3}\unicode[STIX]{x1D6FF}_{l,k+1}+\frac{k}{2k-1}\unicode[STIX]{x1D6FF}_{l,k-1}\right), & \displaystyle\end{eqnarray}$$
    (B 7) $$\begin{eqnarray}\displaystyle & \displaystyle \int _{-1}^{+1}(1-u^{2})P_{k}^{\prime }(u)P_{l}^{\prime }(u)\,\text{d}u=\frac{2k(k+1)}{2k+1}\unicode[STIX]{x1D6FF}_{lk}, & \displaystyle\end{eqnarray}$$
    (B 8) $$\begin{eqnarray}\displaystyle & \displaystyle \int _{-1}^{+1}(1-u^{2})P_{k}^{\prime }(u)P_{l}(u)\,\text{d}u=\frac{2k(k+1)}{2k+1}\left(\frac{\unicode[STIX]{x1D6FF}_{l,k-1}}{2k-1}-\frac{\unicode[STIX]{x1D6FF}_{l,k+1}}{2k+3}\right), & \displaystyle\end{eqnarray}$$
    (B 9) $$\begin{eqnarray}\displaystyle & \displaystyle \int _{-1}^{+1}(1-u^{2})uP_{k}^{\prime }(u)P_{l}^{\prime }(u)\,\text{d}u=\frac{2k(k+1)}{2k+1}\left(\frac{k+2}{2k+3}\unicode[STIX]{x1D6FF}_{l,k+1}+\frac{k-1}{2k-1}\unicode[STIX]{x1D6FF}_{l,k-1}\right),\qquad & \displaystyle\end{eqnarray}$$
    (B 10) $$\begin{eqnarray}\displaystyle & \hspace{-80.00012pt}\displaystyle \int _{-1}^{+1}(1-u^{2})uP_{k}(u)P_{l}^{\prime }(u)\,\text{d}u=\frac{2(k+1)(k+2)(k+3)}{(2k+1)(2k+3)(2k+5)}\unicode[STIX]{x1D6FF}_{l,k+2} & \displaystyle \nonumber\\ \displaystyle & \displaystyle +\,\frac{2k(k+1)}{(2k-1)(2k+1)(2k+3)}\unicode[STIX]{x1D6FF}_{lk}-\frac{2(k-2)(k-1)k}{(2k-3)(2k-1)(2k+1)}\unicode[STIX]{x1D6FF}_{l,k-2}, & \displaystyle\end{eqnarray}$$
    (B 11) $$\begin{eqnarray}\displaystyle & \hspace{-80.00012pt}\displaystyle \int _{-1}^{+1}(1-u^{2})^{2}P_{k}(u)P_{l}^{\prime \prime }(u)\,\text{d}u=\frac{2(k+1)(k+2)(k+3)(k+4)}{(2k+1)(2k+3)(2k+5)}\unicode[STIX]{x1D6FF}_{l,k+2} & \displaystyle \nonumber\\ \displaystyle & \displaystyle -\,\frac{4(k-1)k(k+1)(k+2)}{(2k-1)(2k+1)(2k+3)}\unicode[STIX]{x1D6FF}_{lk}+\frac{2(k-3)(k-2)(k-1)k}{(2k-3)(2k-1)(2k+1)}\unicode[STIX]{x1D6FF}_{l,k-2}. & \displaystyle\end{eqnarray}$$
  6. (vi) And a useful half-integral for our calculations on the Janus sphere:

    (B 12) $$\begin{eqnarray}\int _{0}^{1}P_{k}(u)\,\text{d}u=\left\{\begin{array}{@{}ll@{}}1,\quad & k=0\\ 0,\quad & k=2,4,6,\ldots \\ (-1)^{(m-1)/2}{\displaystyle \frac{m!!}{m(m+1)!!}},\quad & k=1,3,5,\ldots .\end{array}\right.\end{eqnarray}$$

Appendix C

Expressions obtained from the GTD model for $D_{\bot }$ and $D_{\Vert }$ by inverting the truncated 2-by-2 systems corresponding to equations (3.32) and (3.34):

(C 1) $$\begin{eqnarray}\displaystyle D_{\bot } & = & \displaystyle \frac{U^{2}D_{r}^{-1}}{3\unicode[STIX]{x1D706}^{4}(20+\unicode[STIX]{x1D706}^{2})}\left[-5\unicode[STIX]{x1D706}^{2}(9+\unicode[STIX]{x1D706}^{2})-24\unicode[STIX]{x1D706}(5+2\unicode[STIX]{x1D706}^{2})\unicode[STIX]{x1D707}+(120+41\unicode[STIX]{x1D706}^{2}-3\unicode[STIX]{x1D706}^{4})\unicode[STIX]{x1D707}^{2}\right.\nonumber\\ \displaystyle & & \displaystyle \left.+\,\left(45\unicode[STIX]{x1D706}^{3}+8\unicode[STIX]{x1D706}^{2}(15+\unicode[STIX]{x1D706}^{2})\unicode[STIX]{x1D707}-\unicode[STIX]{x1D706}(120+\unicode[STIX]{x1D706}^{2})\unicode[STIX]{x1D707}^{2}\right)\coth \unicode[STIX]{x1D706}\right]\!,\end{eqnarray}$$
(C 2) $$\begin{eqnarray}\displaystyle D_{\Vert } & = & \displaystyle \frac{U^{2}D_{r}^{-1}}{6\unicode[STIX]{x1D706}^{4}(15+\unicode[STIX]{x1D706}^{2})}\left[15\unicode[STIX]{x1D706}^{2}(5+4\unicode[STIX]{x1D706}^{2})+2\unicode[STIX]{x1D706}(255\unicode[STIX]{x1D706}+133\unicode[STIX]{x1D706}^{3}+8\unicode[STIX]{x1D706}^{5})\unicode[STIX]{x1D707}\right.\nonumber\\ \displaystyle & & \displaystyle -\,4(105+43\unicode[STIX]{x1D706}^{2}-\unicode[STIX]{x1D706}^{4})\unicode[STIX]{x1D707}^{2}\nonumber\\ \displaystyle & & \displaystyle +\,\left(-15\unicode[STIX]{x1D706}^{3}+20\unicode[STIX]{x1D706}^{2}(-33+5\unicode[STIX]{x1D706}-7\unicode[STIX]{x1D706}^{2})\unicode[STIX]{x1D707}+12\unicode[STIX]{x1D706}(40+5\unicode[STIX]{x1D706}^{2}-\unicode[STIX]{x1D706}^{4})\unicode[STIX]{x1D707}^{2}\right)\coth \unicode[STIX]{x1D706}\nonumber\\ \displaystyle & & \displaystyle \left.+\,6\unicode[STIX]{x1D706}^{2}(\unicode[STIX]{x1D706}-2\unicode[STIX]{x1D707})(5\unicode[STIX]{x1D707}-10\unicode[STIX]{x1D706}-\unicode[STIX]{x1D706}^{2}\unicode[STIX]{x1D707})\coth ^{2}\unicode[STIX]{x1D706}\right]\!.\end{eqnarray}$$

References

Agudo-Canalejo, J., Illien, P. & Golestanian, R. 2018 Phoresis and enhanced diffusion compete in enzyme chemotaxis. Nano Lett. 18, 27112717.Google Scholar
Anderson, J. L. 1989 Colloid transport by interfacial forces. Annu. Rev. Fluid Mech. 21, 6199.Google Scholar
Batchelor, G. K. 1976 Brownian diffusion of particles with hydrodynamic interaction. J. Fluid Mech. 74, 129.Google Scholar
Bearon, R. N., Bees, M. A. & Croze, O. A. 2012 Biased swimming cells do not disperse in pipes as tracers: a population model based on microscale behaviour. Phys. Fluids 24, 121902.Google Scholar
Berg, H. C. 1975 Chemotaxis in bacteria. Annu. Rev. Biophys. Bioengng 4, 119136.Google Scholar
Bickel, T., Majee, A. & Würger, A. 2013 Flow pattern in the vicinity of self-propelling hot Janus particles. Phys. Rev. E 88, 012301.Google Scholar
Bickel, T., Zecua, G. & Würger, A. 2014 Polarization of active Janus particles. Phys. Rev. E 89, 050303(R).Google Scholar
Brady, J. 2011 Particle motion driven by solute gradients with application to autonomous motion: continuum and colloidal perspectives. J. Fluid Mech. 667, 216259.Google Scholar
Brown, A. & Poon, W. 2014 Ionic effects in self-propelled Pt-coated Janus swimmers. Soft Matt. 10, 40164027.Google Scholar
Córdova-Figueroa, U. M. & Brady, J. F. 2008 Osmotic propulsion: the osmotic motor. Phys. Rev. Lett. 100, 158303.Google Scholar
Córdova-Figueroa, U. M., Brady, J. F. & Shklyaev, S. 2013 Osmotic propulsion of colloidal particles via constant surface flux. Soft Matt. 9, 63826390.Google Scholar
Ebbens, S., Gregory, D. A., Dunderdale, G., Howse, J. R., Ibrahim, Y., Liverpool, T. B. & Golestanian, R. 2014 Electrokinetic effects in catalytic platinum-insulator Janus swimmers. Eur. Phys. Lett. 106, 58003.Google Scholar
Ebbens, S., Tu, M.-H., Howse, J. R. & Golestanian, R. 2012 Size dependence of the propulsion velocity for catalytic Janus-sphere swimmers. Phys. Rev. E 85, 020401.Google Scholar
Ebbens, S. J. & Howse, J. R. 2011 Direct observation of the direction of motion for spherical catalytic swimmers. Langmuir 27, 1229312296.Google Scholar
Einstein, A. 1905 Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen. Ann. Phys. 322, 549560.Google Scholar
Frankel, I. & Brenner, H. 1991 Generalized Taylor dispersion phenomena in unbounded homogeneous shear flows. J. Fluid Mech. 230, 147181.Google Scholar
Frankel, I. & Brenner, H. 1993 Taylor dispersion of orientable Brownian particles in unbounded homogeneous shear flows. J. Fluid Mech. 255, 129156.Google Scholar
Geiseler, A., Hanggi, P., Marchesoni, F., Mulhern, C. & Savel’ev, S. 2016 Chemotaxis of artificial microswimmers in active density waves. Phys. Rev. E 94, 012613.Google Scholar
Golestanian, R. 2012 Collective behavior of thermally active colloids. Phys. Rev. Lett. 108, 3.Google Scholar
Golestanian, R., Liverpool, T. B. & Ajdari, A. 2005 Propulsion of a molecular machine by asymmetric distribution of reaction products. Phys. Rev. Lett. 94, 220801.Google Scholar
Golestanian, R., Liverpool, T. B. & Ajdari, A. 2007 Designing phoretic micro- and nanoswimmers. New J. Phys. 9, 126.Google Scholar
Hill, N. A. & Bees, M. A. 2002 Taylor dispersion of gyrotactic swimming micro-organisms in a linear flow. Phys. Fluids. 14, 25982605.Google Scholar
Howse, J. R., Jones, R. A. L., Ryan, A. J., Gough, T., Vafabakhsh, R. & Golestanian, R. 2007 Self-motile colloidal particles: from directed propulsion to random walk. Phys. Rev. Lett. 99, 048102.Google Scholar
Izri, Z., van der Linden, M. N., Michelin, S. & Dauchot, O. 2014 Self-propulsion of pure water droplets by spontaneous Marangoni stress driven motion. Phys. Rev. Lett. 113, 248302.Google Scholar
Jiang, H.-R., Yoshinaga, N. & Sano, M. 2010 Active motion of a Janus particle by self-thermophoresis in a defocused laser beam. Phys. Rev. Lett. 105, 268302.Google Scholar
Jülicher, F. & Prost, J. 2009 Generic theory of colloidal transport. Eur. Phys. J. E 29, 2736.Google Scholar
Khair, A. S. 2013 Diffusiophoresis of colloidal particles in neutral solute gradients at finite Péclet number. J. Fluid Mech. 731, 6494.Google Scholar
Lauga, E. 2016 Bacterial hydrodynamics. Annu. Rev. Fluid Mech. 48, 105130.Google Scholar
Lauga, E. & Powers, T. R. 2009 The hydrodynamics of swimming micro-organisms. Rep. Prog. Phys. 72, 096601.Google Scholar
Manela, A. & Frankel, I. 2003 Generalized Taylor dispersion in suspensions of gyrotactic swimming micro-organisms. J. Fluid Mech. 490, 99127.Google Scholar
Michelin, S. & Lauga, E. 2014 Phoretic self-propulsion at finite Péclet numbers. J. Fluid Mech. 747, 572604.Google Scholar
Michelin, S. & Lauga, E. 2015 Autophoretic locomotion from geometric asymmetry. Eur. Phys. J. E 38, 7.Google Scholar
Michelin, S., Lauga, E. & Bartolo, D. 2013 Spontaneous autophoretic motion of isotropic particles. Phys. Fluids 25, 061701.Google Scholar
Nelson, B. J., Kaliakatsos, I. K. & Abbott, J. J. 2010 Microrobots for minimally invasive medicine. Annu. Rev. Biomed. Engng 12, 5585.Google Scholar
Nelson, P. C. 2008 Biological Physics: Energy, Information, Life. W. H. Freeman.Google Scholar
Palacci, J., Sacanna, S., Kim, S.-H., Yi, G.-R., Pine, D. J. & Chaikin, P. M. 2014 Light-activated self-propelled colloids. Phil. Trans. R. Soc. A 372, 20130372.Google Scholar
Palacci, J., Sacanna, S., Steinberg, A. P., Pine, D. J. & Chaikin, P. M. 2013 Living crystals of light-activated colloidal surfers. Science 339, 936940.Google Scholar
Paxton, W. F., Kistler, K. C., Olmeda, C. C., Sen, A., St. Angelo, S. K., Cao, Y., Mallouk, T. E., Lammert, P. E. & Crespi, V. H. 2004 Catalytic nanomotors: autonomous movement of striped nanorods. J. Am. Chem. Soc. 126, 1342413431.Google Scholar
Pedley, T. J. & Kessler, J. O. 1990 A new continuum model for suspensions of gyrotactic micro-organisms. J. Fluid Mech. 212, 155182.Google Scholar
Pohl, O. & Stark, H. 2014 Dynamic clustering and chemotactic collapse of self-phoretic active particles. Phys. Rev. Lett. 112, 238303.Google Scholar
Popescu, M. N., Dietrich, S., Tasinkevych, M. & Ralston, J. 2010 Phoretic motion of spheroidal particles due to self-generated solute gradients. Eur. Phys J. E 31, 351367.Google Scholar
Popescu, M. N., Tasinkevych, M. & Dietrich, S. 2011 Pulling and pushing a cargo with a catalytically active carrier. Eur. Phys. Lett. 95, 28004.Google Scholar
Sabass, B. & Seifert, U. 2012 Dynamics and efficiency of a self-propelled, diffusiophoretic swimmer. J. Chem. Phys. 136, 064508.Google Scholar
Saha, S., Golestanian, R. & Ramaswamy, S. 2014 Clusters, asters, and collective oscillations in chemotactic colloids. Phys. Rev. E 89, 062316.Google Scholar
Saragosti, J., Silberzan, P. & Buguin, A. 2012 Modeling E. coli tumbles by rotational diffusion. Implications for chemotaxis. PLoS One 7, 16.Google Scholar
Schmitt, M. & Stark, H. 2013 Swimming active droplet: a theoretical analysis. Eur. Phys. Lett. 101, 44008.Google Scholar
Shklyaev, S., Brady, J. F. & Córdova-Figueroa, U. M. 2014 Non-spherical osmotic motor: chemical sailing. J. Fluid Mech. 748, 488520.Google Scholar
von Smoluchowski, M. 1906 Zur kinetischen Theorie der Brownschen Molekularbewegung und der Suspensionen. Ann. Phys. 326, 756780.Google Scholar
Stone, H. A. & Samuel, A. D. T. 1996 Propulsion of microorganisms by surface distortions. Phys. Rev. Lett. 77, 41024104.Google Scholar
Thutupalli, S., Seemann, R. & Herminghaus, S. 2011 Swarming behavior of simple model squirmers. New J. Phys. 13, 073021.Google Scholar
Wang, W., Duan, W., Sen, A. & Mallouk, T. E. 2013 Catalytically powered dynamic assembly of rod-shaped nanomotors and passive tracer particles. Proc. Natl. Acad. Sci. USA 110, 1774417749.Google Scholar
Figure 0

Figure 1. Sketch of the problem. (a) Axisymmetric active particle with axis of symmetry along $z$. The surface properties (surface activity $\unicode[STIX]{x1D70E}$, substrate mobility $\unicode[STIX]{x1D707}_{s}$ and product mobility $\unicode[STIX]{x1D707}_{p}$) are functions of the polar angle $\unicode[STIX]{x1D703}$ alone. (b) Active particle placed in a uniform substrate gradient, $\unicode[STIX]{x1D735}s^{\infty }$, which lies in the $yz$-plane. The presence of the background gradient, together with the catalytic coating of the swimmer, induce local tangential chemical gradients on the surface of the sphere (light grey arrows) which give rise to an apparent fluid slip velocity (white arrows) leading to translation (linear velocity $\boldsymbol{V}$) and rotation (angular velocity $\unicode[STIX]{x1D74E}$) of the particle.

Figure 1

Figure 2. Schematic representation of a Janus spherical swimmer with surface properties $\{\unicode[STIX]{x1D70E}^{U},\unicode[STIX]{x1D707}_{s}^{U},\unicode[STIX]{x1D707}_{p}^{U}\}$ on the upper hemisphere and $\{\unicode[STIX]{x1D70E}^{L},\unicode[STIX]{x1D707}_{s}^{L},\unicode[STIX]{x1D707}_{p}^{L}\}$ on the lower hemisphere. The linear and angular velocities of this phoretic swimmer are given by (2.50)–(2.53).

Figure 2

Figure 3. Comparison of the distribution of swimmer orientations obtained from numerical simulations with the predictions of the continuum model. (a) Histograms of final orientations for 10 000 swimmers whose initial orientations were chosen from a uniform distribution, after being evolved in time for 500 time steps of size $\unicode[STIX]{x0394}t=10^{-2}D_{r}^{-1}$. The red solid line corresponds to the analytical expression from (3.21) for the steady-state distribution of swimmer orientations. (b) Sample of 1000 swimmer orientations chosen from the steady-state distribution. Each swimmer orientation corresponds to a unit vector which is represented here as a blue dot on the unit sphere. From top to bottom the rotational Péclet number is $\unicode[STIX]{x1D706}=0.1,3,10$.

Figure 3

Figure 4. Time evolution of 200 independent trajectories, all starting from the origin at $t=0$. Each snapshot is taken after a time interval of $3D_{r}^{-1}$, and the chemical gradient is along the horizontal axis. The large dots are placed at regular distances of $3V_{s}D_{r}^{-1}$ from the origin which confirms that the cloud of phoretic swimmers drifts along the direction of the chemical gradient at the mean swimming velocity ($V_{s}$), as predicted by the continuum model. This simulation corresponds to parameters $\unicode[STIX]{x1D706}=5$ and $\unicode[STIX]{x1D707}=3$, and the time step used was $\unicode[STIX]{x0394}t=10^{-2}D_{r}^{-1}$.

Figure 4

Figure 5. Comparison of the values obtained from numerical simulations against those predicted by GTD theory for the dependence of component $D_{\bot }$ of the active diffusivity tensor on $\unicode[STIX]{x1D706}$. The simulations were run with 100 samples of 1000 swimmers using a time step $\unicode[STIX]{x0394}t=10^{-2}D_{r}^{-1}$ over a period $2D_{r}^{-1}$. Each sample was used to calculate one value of the diffusivity $D_{\bot }$, and then averaged to obtain one data point. The error bars represent one standard deviation amongst the 100 values obtained for the diffusivity. The solid lines represent the predictions of the continuum model from (3.32)–(3.35) and (3.38)–(3.39). In all plots the diffusivity is non-dimensionalised by $U^{2}D_{r}^{-1}$. The four graphs correspond to parameter values $\unicode[STIX]{x1D707}=0$, $\unicode[STIX]{x1D707}=1$, $\unicode[STIX]{x1D707}=5$ and $\unicode[STIX]{x1D707}=10$.

Figure 5

Figure 6. Iso-values of the diffusivities $D_{\bot }$ and $D_{\Vert }$ computed using GTD theory in the $(\unicode[STIX]{x1D706},\unicode[STIX]{x1D707})$ plane. Both diffusivities are normalised using the mean swimming velocity and the time scale of rotational diffusion, i.e. $D\mapsto D/(V_{s}^{2}D_{r}^{-1})$. These contour plots are obtained with a value of $\unicode[STIX]{x1D708}=1.5$ for the direct chemotactic index, but there are no qualitative changes as $\unicode[STIX]{x1D708}$ is varied (not shown). We observe consistent decay of both diffusivities as $\unicode[STIX]{x1D706}\rightarrow \infty$ at a fixed value of $\unicode[STIX]{x1D707}$, as well as a saddle point in diffusivity for $\unicode[STIX]{x1D706},\unicode[STIX]{x1D707}=O(1)$.

Figure 6

Figure 7. The ratio $D_{\bot }/D_{\Vert }$ computed using GTD theory in the $(\unicode[STIX]{x1D706},\unicode[STIX]{x1D707})$ plane. As expected, the contour $D_{\bot }/D_{\Vert }=1$ passes through the origin, because the system reverts to isotropic diffusion in the limit $\unicode[STIX]{x1D706},\unicode[STIX]{x1D707}\rightarrow 0$. Increasing $\unicode[STIX]{x1D707}$ above this contour leads to a ratio $D_{\bot }/D_{\Vert }<1$, meaning that indirect chemotactic sedimentation enhances diffusion along the chemical gradient relative to diffusion in the normal plane. On the other hand, increasing $\unicode[STIX]{x1D706}$ to the right of this contour leads to a ratio $D_{\bot }/D_{\Vert }>1$, meaning that chemotactic alignment suppresses diffusion along the chemical gradient relative to diffusion in the normal plane.