Pairing, waltzing and scattering of chemotactic active colloids

TL;DR

This paper investigates the complex dynamical interactions of chemotactic active colloids, revealing bound states, scattering behavior, and bifurcations driven by their chemical interactions and propulsion speeds.

Contribution

It introduces a detailed analysis of two-body interactions of chemotactic colloids, highlighting novel bound states, scattering phenomena, and the role of non-reciprocal dynamical forces.

Findings

Chemotactic colloids can form active dimers with translational motion.

Anti-chemotactic colloids always scatter, preventing bound states.

Transitions between bound states depend on propulsion speed and initial conditions.

Abstract

An interacting pair of chemotactic (anti-chemotactic) active colloids, that can rotate their axes of self-propulsion to align {parallel (anti-parallel)} to a chemical gradient, shows dynamical behaviour that varies from bound states to scattering. The underlying two-body interactions are purely dynamical, non-central, non-reciprocal, and controlled by changing the catalytic activity and phoretic mobility. Mutually chemotactic colloids trap each other in a final state of fixed separation; the resulting `active dimer' translates. A second type of bound state is observed where the polar axes undergo periodic cycles leading to phase-synchronised circular motion around a common point. These bound states are formed depending on initial conditions and can unbind on increasing the speed of self propulsion. Mutually anti-chemotactic swimmers always scatter apart. We also classify the fixed points underlying the bound states, and the bifurcations leading to transitions from one type of bound state to another, for the case of a single swimmer in the presence of a localised source of solute.

Figures

Equations

$$A^{(i)}({\theta }_i)=a^{(i)}{\sum }_{\ell }{\alpha }_{\ell }^{(i)}{P}_{\ell }({\theta }_i),M^{(i)}({\theta }_i)$$

$$-D{\nabla }^2\rho =0,-D{\nabla }_{\perp }\rho ({\mathbf{r}}_i,t)=A^{(i)}({\mathbf{r}}_i).$$

$${\mathbf{V}}_i=-\frac{1}{4\pi }\int \text{mbox}ds{\mathbf{v}}_i({\mathbf{r}}_i),{\mathbf{\omega }}_i=-\frac{3}{16\pi {\sigma }_i}\int \text{mbox}ds{\mathbf{r}}_i\times {\mathbf{v}}_i({\mathbf{r}}_i).$$

$${\dot{\theta }}_{n1}=\frac{a^{(1)}{m}^{(1)}}{D\sigma }\omega \left({\Delta }_1,{\Delta }_2,R\right),$$

$$\dot{R_1}=\frac{a^{(1)}{m}^{(1)}}{D}V\left({\Delta }_1,{\Delta }_2,R\right),$$

$$\dot{{\beta }_1}=\frac{a^{(1)}{m}^{(1)}}{D}\Omega \left({\Delta }_1,{\Delta }_2,R\right).$$

$$\omega ({\Delta }_1,{\Delta }_2,R)$$

$$V({\Delta }_1,{\Delta }_2,R)$$

$$\Omega ({\Delta }_1,{\Delta }_2,R)$$

$$\omega$$

$$R_0\approx \sigma \sqrt{\frac{3}{h_0}}{\left(-{\alpha }_1^{(1)}-\frac{2{\alpha }_2^{(1)}{\mu }_1^{(1)}}{5{\mu }_0^{(1)}-{\mu }_2^{(1)}}\right)}^{-\frac{1}{2}}.$$

$$\dot{\delta R}(t)=V(0,0,R=R_0\pm \delta R)\equiv -\nabla \psi .$$

$$\omega ({\Delta }_{10},0,R_0)-\Omega ({\Delta }_{10},0,R_0)=0,V({\Delta }_{10},0,R_0)=0.$$

$$V({\Delta }_{10},{\Delta }_{20},R_0)+V({\Delta }_{20},{\Delta }_{10},R_0)=0,$$

$$\Omega ({\Delta }_{10},{\Delta }_{20},R_0)-\omega ({\Delta }_{10},{\Delta }_{20},R_0)=\Omega ({\Delta }_{20},{\Delta }_{10},R_0)-\omega ({\Delta }_{20},{\Delta }_{10},R_0)=0.$$

$$\cot {\varphi }_{10}=\frac{V({\Delta }_{10},{\Delta }_{20},R_0)}{R_0\Omega ({\Delta }_{10},{\Delta }_{20},R_0)},\cot {\varphi }_{20}=\frac{V({\Delta }_{20},{\Delta }_{10},R_0)}{R_0\Omega ({\Delta }_{20},{\Delta }_{10},R_0)}.$$

$${\nabla }^2\rho$$

$$-D{\nabla }_{\perp }\rho ({\mathbf{r}}_i)$$

$$\bar{\rho }=\sum _n{\rho }_n^{(1)}+\sum _n{\rho }_n^{(2)}.$$

$$-{\nabla }_{\perp }{\rho }_n^{(1)}({\mathbf{r}}_1)$$

$$-{\nabla }_{\perp }{\rho }_n^{(2)}({\mathbf{r}}_2)$$

$$-{\nabla }_{\perp }{\rho }_n^{(i)}({\mathbf{r}}_i)\equiv \sum _{\ell m}{A}_{n,\ell m}^{(i)}{Y}_{\ell m}({\mathbf{r}}_i).$$

$$A_{n,\ell m}^{(1)}=\sum _{k=0}^{n-1}\sum _{j=1}^2{A}_{k,jm}^{(2)}{N}_{\ell jm}^{12}.$$

$${\rho }_n^{(2)}$$

$${\theta }_2({\theta }_1,R_1,R)$$

$$R_2({\theta }_1,R_1,R)$$

$${\varphi }_2$$

$$N_{i\ell m}^{12}=\sigma {\partial }_{R_1}\int Y_{\ell m}^{\ast }(\theta ,\varphi )\frac{Y_{im}\left[{\theta }_2({\theta }_1,R_1,R),\varphi \right]}{R_2({\theta }_1,R_1,R{)}^{\ell +1}}\text{mbox}d{\mathbf{r}}_1.$$

$${\rho }_0+{\rho }_1^{(1)}=\frac{a^{(2)}\sigma }{D}\sum _{\ell m}\frac{4\pi }{2\ell +1}\left[\frac{\sigma R_1^{\ell }}{R^{\ell +1}}+\frac{\ell }{\sigma (\ell +1)}{\left(\frac{{\sigma }^2}{R{R}_1}\right)}^{\ell +1}\right]Y_{\ell m}(\theta ,\varphi )Y_{\ell m}^{\ast }(\beta ,0).$$

$$\omega ({\Delta }_1,{\Delta }_2,r)=\frac{3{\sigma }^2{\alpha }_0^{(2)}{\mu }_1^{(1)}}{8R^2}\sin {\Delta }_1+\frac{3{\sigma }^3{\alpha }_0^{(2)}{\mu }_2^{(1)}}{8R^3}\sin 2{\Delta }_1-\frac{3{\sigma }^3{\alpha }_1^{(2)}{\mu }_1^{(1)}}{32R^3}\left[\sin {\Delta }_1-3\sin ({\Delta }_1+{\Delta }_2)\right]$$

$$-\frac{3{\sigma }^4{\alpha }_1^{(2)}{\mu }_2^{(1)}}{32R^4}\left[5\sin (2{\Delta }_1+{\Delta }_2)-\sin (2{\Delta }_1-{\Delta }_2)\right]+\frac{3{\sigma }^4{\alpha }_2^{(2)}{\mu }_1^{(1)}}{64R^4}\left[5\sin ({\Delta }_1+2{\Delta }_2)+2\sin {\Delta }_1+\sin ({\Delta }_1-2{\Delta }_2)\right]$$

$$+\frac{{\sigma }^5{\alpha }_2^{(2)}{\mu }_2^{(1)}}{64R^5}\left[3\sin (2{\Delta }_1-2{\Delta }_2)+35\sin (2{\Delta }_1+2{\Delta }_2)+10\sin 2{\Delta }_1\right]$$

$$V({\Delta }_1,{\Delta }_2,r)=-\left({\mu }_0^{(1)}-\frac{{\mu }_2^{(1)}}{20}\right)\frac{{\sigma }^2{\alpha }_0^{(2)}}{R^2}+\frac{3{\sigma }^2{\mu }_2^{(1)}{\alpha }_0^{(2)}}{20R^2}\cos 2{\Delta }_1,$$

$$-\frac{2{\sigma }^3{\alpha }_0^{(2)}{\mu }_1^{(1)}}{3R^3}\cos {\Delta }_1+\frac{{\sigma }^3{\alpha }_1^{(2)}{\mu }_{10}}{R^3}\cos {\Delta }_2+\frac{{\sigma }^3{\mu }_2^{(1)}{\alpha }_1^{(2)}}{80R^3}\left[9\cos (2{\Delta }_1+{\Delta }_2)+3\cos (2{\Delta }_1-{\Delta }_2)+4\cos {\Delta }_2\right]$$

$$+\frac{{\sigma }^4{\alpha }_1^{(2)}{\mu }_1^{(1)}}{4R^4}[3\cos ({\Delta }_1+{\Delta }_2)+\cos ({\Delta }_1-{\Delta }_2)]-\frac{{\sigma }^4{\alpha }_2^{(2)}{\mu }_{10}}{4R^4}[1+3\cos 2{\Delta }_2]$$

$$+\frac{{\sigma }^4{\mu }_2^{(1)}{\alpha }_2^{(2)}}{160R^4}[6\cos 2{\Delta }_1+15\cos (2{\Delta }_1+2{\Delta }_2)+3\cos (2{\Delta }_1-2{\Delta }_2)+6\cos 2{\Delta }_2+2]$$

$$-\frac{{\sigma }^5{\mu }_1^{(1)}{\alpha }_2^{(2)}}{6R^5}[2\cos {\Delta }_1+5\cos ({\Delta }_1+2{\Delta }_2)+\cos ({\Delta }_1-2{\Delta }_2)].$$

Full text

1 Introduction

Gradients in the concentration of a solute along the surface of a particle in a fluid produce pressure differences, and thus a slip flow, parallel to the surface. If not anchored, the particle moves through the fluid in the direction opposite to the slip velocity. This phenomenon, known as diffusiophoresis [1, 2], finds dramatic application in the autonomous motility of active Janus colloids [3, 4, 5, 6, 7], that generate the required gradient of solute molecules themselves via a surface patch of catalyst that decomposes a species in the ambient medium. When placed in an externally imposed gradient of solute, these artificial microswimmers rotate their intrinsic polarity, as defined by the catalyst patch, to point – and thus to translate – towards or away from regions of high concentration [8, 10, 11], emulating chemotaxis [12].

Active colloids have emerged as a versatile constituent of synthetic active matter [13], as the surface profiles of catalyst concentration and solute-colloid interaction can be engineered to yield a variety of designed swimmer properties [8, 9]. Recent studies have explored their collective behaviour including pattern formation, motility induced phase separation, chemotactic collapse [8, 14, 15, 16, 17], dynamic swarms [18, 19, 20, 21, 22, 23, 24] and oscillatory clustering [25, 15]. A study of two-particle dynamics with a directory of possible interactions, crucial for a full understanding of the possible collective behaviour in these systems, is the subject of this paper. Earlier studies on pair interactions between motile particles include work on the role of hydrodynamic fields [26, 27, 28, 29] and a general classification scheme [30]. Our focus is on interactions mediated by the diffusing field of the reaction products.

Our pair of active colloids are spheres ($i=1,2$) of radius $\sigma$ characterised by a catalytic coating of surface density $A^{(i)}$ and a mobility $M^{(i)}$. The mobility $M^{(i)}$ is determined by the interaction potential between the surface of the sphere and the molecules of the product. Both vary from point to point on the surface of the sphere. We therefore decompose them into Legendre polynomials, assuming for simplicity that they are axisymmetric with a common polar axis:

[TABLE]

where $a^{(i)}$ is the average rate of production of product molecules per unit area, $m^{(i)}$ sets the scale of the mobility and $\theta_{i}$ is the colatitude on the sphere. The essential phenomenology is conveyed by the simple case where we retain $\ell=0,1,2$ for $\alpha^{(i)}_{\ell}$ and $\mu^{(i)}_{\ell}$. The surface patterns are represented by the sets of numbers $\{\mu^{(i)}_{0},\mu^{(i)}_{1},\mu^{(i)}_{2}\}$ and $\{\alpha^{(i)}_{0},\alpha^{(i)}_{1},\alpha^{(i)}_{2}\}$, with each such assignment defining a distinct swimmer design. As shown in [8], an active colloid with $\mu^{(i)}_{\ell}\neq 0$ for $\ell>0$, generates chemotactic angular velocity that can turn its polar axis to align with a gradient. Parallel and antiparallel alignment are respectively called chemotactic and anti-chemotactic.

Here is a summary of our results. We consider two distinct cases – a swimmer interacting with a fixed source of solute and two swimmers that are both free to move. (i) A colloid which is (positively) chemotactic at every separation from the source, $\mu_{1}^{(1)}<0,\mu_{2}^{(1)}<0$, aligns with the gradient and is trapped by the source when self-propulsion and phoretic repulsion balance. (ii) A colloid that exhibits mixed response to a the source, chemotactic at large separation and antichemotactic when within a few particle radii from source, executes what appears to be periodic motion in closed orbits. (iii) Increasing the speed of the swimmer or increasing the impact parameter drives a transition from bound state to scattering. (iv) An anti-chemotactic colloid is always scattered. Turning to the case of two swimmers we find: (v) Two interacting swimmers can form a stationary or moving dimer stabilised by a balance of self propulsion and phoretic repulsion, with a fixed distance between them. (vi) Two swimmers, at least one of which has a mixed chemotactic response to the other, can form bound states where they revolve around a common point on closed phase synchronised orbits while their polar axes are locked at a finite inclination. (vii) Two anti-chemotactic swimmers always scatter off one another. The steady states have been summarised in state diagrams (see fig. 2) obtained by varying total surface activity $a^{(i)}$ and composition of the motility coat. Finally, we also show that the bound states are robust to the presence of thermal or other fluctuations, so that signatures of these states should be seen in experiments.

The article is organised as follows: In section 2 we present the equations of motion for the swimmers for a given separation and relative orientation and outline steps followed to obtain the dynamics. In section 3 we discuss the bound states obtained when one swimmer interacts with a source of solute. In section 4 we categorize the dynamics shown by two mobile swimmers.

2 Interacting active colloids

Consider two swimmers in a fluid medium constrained to move in a plane that also contains their polar axes. Reactants are converted into products when they come in contact with the enzyme coated colloidal surface, resulting in a spatiotemporal distribution $\rho(\mathbf{r},t)$ of products. The product diffuses freely in the bulk of the fluid with diffusion coefficient $D$. There is a normal flux of the product molecules on surfaces of the swimmers at every point $\mathbf{r}_{i}$

[TABLE]

The activity $a^{(i)}$, can depend on the availability of reactant locally in which case it can be approximated by the Michaelis-Menten [31] factor. In this work we do not consider this aspect and work with constant $a^{(i)}$.

Under conditions of vanishing Reynolds number, the mobility and chemical gradient lead to the establishment of a slip velocity $\mathbf{v}_{i}=M^{(i)}\mathbf{\nabla}_{\parallel}\rho$. The slip velocity produces linear and angular velocities which for a sphere are given by the following surface integrals

[TABLE]

The linear velocity $\mathbf{V}_{i}$ subsumes both self propulsion and phoretic interactions, while $\bm{\omega_{i}}$ is the angular velocity of the swimmer. The sign of $\bm{\omega}_{i}$ determines the ‘chemotactic’ response of the swimmer i.e. whether it aligns to point up or down the gradient [8]. Equations of motion for the positions $\mathbf{R}_{1,2}$ and orientations $\theta_{n1},\theta_{n2}$ of the swimmers in the lab frame are built by first solving eq. 2 for $\rho$. The solution for $\rho$ gives $\mathbf{v}_{i}$ which is inserted into eq. 3 to calculate the velocities. Translation and rotation invariance of space restricts the dependence of these velocities on $R$ and the orientation of the polar axes, $\Delta_{1}$ and $\Delta_{2}$, measured with respect to the line joining the sphere centres (see fig. 1). In Appendix B we show that the field $\rho$ can be calculated to a desired degree of precision as a perturbation series in powers of $\sigma/R$, where $\sigma$ is the colloid size. The calculation uses a method of reflections, analogous to the method of images used in electrostatics. For simplicity, we construct the dynamics for just one set of reflections; corrections from subsequent reflections contribute to progressively larger powers in $\sigma/R$ and can be ignored.

The linear velocity of swimmer 1 in spherical polar coordinates in the frame of reference of swimmer 2 in terms of radial and angular velocities $V$ and $\Omega$ is, $\mathbf{V}_{1}=V\hat{R}+\Omega R\hat{\beta}_{1}$. $\hat{R}$ and $\hat{\beta}_{1}$ are unit vectors in the radial and tangential directions. We assume the swimmer size $\sigma$ to scale distances, while the velocities are scaled by $a^{(1)}m^{(1)}/D$. The equations of motion for swimmer 1 in terms of the separation $R$ and the relative orientations $\Delta_{1}=\theta_{n1}-\beta_{1}$ and $\Delta_{2}=\theta_{n2}-\beta_{1}-\pi$ are

[TABLE]

The velocity of swimmer 1 is a linear superposition of contribution from different harmonics of the catalytic coat of swimmer 2, we first discuss the case when $\alpha_{0}^{(2)}=1$ and $\alpha^{(2)}_{\ell}=0$ for $\ell>0$. The angular velocity for the polar axis is

[TABLE]

At a given position, the polar axis rotates till it aligns parallel ($\Delta_{1}=0$) or anti-parallel ($\Delta_{1}=\pi$) with the local concentration gradient. The nature of the chemotactic response of the swimmer is determined by the sign of the factor $\mu_{1}^{(1)}+\mu_{2}^{(1)}/R$, for small deviations around alignment. When $\mu_{1}^{(1)},\mu_{2}^{(1)}$ carry the same sign, the chemotactic response stays unchanged at all distances from the source, being chemotactic (anti-chemotactic) for negative (positive) sign. For $\mu_{1}^{(1)}<0$ and $\mu_{2}^{(1)}>0$, the response depends on radial distance and reverses when the swimmer approaches the source, leading to orbiting states which we will discuss in the next section. The radial and angular velocities are

[TABLE]

Here the dimensionless parameter $h_{0}={a^{(1)}}/{a^{(2)}}$ is the ratio of the total catalytic activity of swimmers 1 and 2. In eqs. 8 and 9, the terms proportional to $h_{0}$ arise from self propulsion and the rest are due to phoretic response to product field of swimmer 2. As swimmer 1 responds to the product field generated by itself and that produced by swimmer 2 with the same mobility, $h_{0}$ measures the relative strength of self propulsion and interaction. The higher order terms in $\omega$ and $\mathbf{V}$ are always of the form of trigonometric functions of $p\Delta_{1}$, where $p$ is an integer. When the swimmers are sufficiently separated, the polarity relaxes before the position, as it can do so independently of the separation. Assuming this separation of timescale between the radial position and polarity, the mutual radial velocity can be obtained by substituting $\Delta_{1}=0$ in eq. 8. The non monotonic radial velocity, which is a consequence of varying mobility and self propulsion, affords the possibility of it vanishing at finite $R$ producing bound states. The most dominant contribution to the radial velocity is the $1/R^{2}$ term which causes a repulsion of the swimmer from the source for $\mu^{(1)}_{0}>0$. For $\mu^{(1)}_{0}>0$ and $\alpha^{(1)}_{1}<0$, self propulsion is directed opposite to phoretic repulsion thus leading to formation of a bound state. The details of these calculations are presented in the next section.

Finally, we write down a few of the terms in $\omega$ that are present when the anisotropy of chemical field produced by swimmer 2 is considered, and present the full equations in the appendix.

[TABLE]

The colloid is a micron sized particle subject to fluctuations due to thermal as well as non-thermal sources, so that the dynamics should be supplemented with appropriate noises. We present a discussion on the effect of noise on the dynamics in a later section. In the next section we can discuss two possible bound states and explore swimmer designs that lead to one or the other.

3 Single active colloid in a product source

We first consider the case when, $m^{(2)}=0$ and $m^{(1)}\neq 0$; swimmer 1 interacts with a tethered swimmer 2, that we now call source. In this section, an isotropic source is one that is uniformly coated with catalysts, $\alpha_{\ell}^{(2)}=0$ for $\ell>0$; an anisotropic source has at least one nonzero $\alpha_{\ell}^{(2)}$ for $\ell>0$.

3.1 Bound final state

We first discuss a bound state formed by a swimmer with mobility $\mu_{1}^{(1)}<0$ and $\mu_{2}^{(1)}<0$. The swimmer can be trapped on a collection of points where both translational and angular velocities vanish. A swimmer that impinges on a source aligns with the local concentration gradient, provided it spends enough time in its vicinity. It is trapped if the radial velocity vanishes. For an isotropic source, the trapping surface is a circle of radius $R_{0}$ centred at the source. The final orientation is $\Delta_{1}=0$, so that $R_{0}$ satisfies the condition $V(0,0,R_{0})=0$. Using eq. 8, and ignoring terms of order higher than $\sigma^{2}/R_{0}^{2}$, we find

[TABLE]

Swimmers designed with parameters that allow a positive $R_{0}$ are trapped. $R_{0}$ depends on swimmer design and can be tuned by varying the ratio $\alpha^{(1)}_{0}/\alpha^{(1)}_{2}$, for example. Swimmers of distinct designs are thus trapped at different distances from the source. This provides a method for sorting the colloids experimentally.

The motion of the swimmer in a chemical field can be equivalent to colloid motion in an externally applied attractive force field in a dissipative medium. This is seen as follows: with the assumption that $\Delta_{1}$ has relaxed to its stationary value, the equation of motion of the colloid for deviations away from $R_{0}$, denoted by $\delta R$ is

[TABLE]

The source produces a confining potential $\psi(R)=\int_{R_{0}}^{R}\mbox{d}R^{\prime}V(0,0,R^{\prime})$, when $\partial_{R}^{2}\psi(R)>0$. An isotropic source produces a radially symmetric potential, see fig. 4 (c). Near an anisotropic source, $\omega$ vanishes when $\theta_{n1}-\theta_{n2}=\pi$, thus creating an anisotropic confining potential with two traps along the polar axis of the source. The swimmer gets trapped in one or the other depending on the initial conditions. To illustrate the strength of the confinement we plot the anisotropic potential produced along the source axis in fig. 4 (d). The region of state space where the bound state is formed is shown in pink in fig. 2 (a) - (c). Note a swimmer is trapped either when both $\mu^{(1)}_{1,2}<0$ or either one of them is both negative and sufficiently larger in magnitude than the other.

3.2 Orbits

A swimmer with $\mu_{2}^{(1)}>0$ and $\mu_{1}^{(1)}<0$, can execute persistent periodic orbits around the source. This particular choice of parameters ensures that the swimmer is chemotactic at large distance and anti-chemotactic near the source. Close to an isotropic source, it tries to align with the local gradient and becomes trapped in a limit cycle at a constant distance $R_{0}$ (indicated in green in fig. 2 (a) and (b)) while maintaining a constant angle $\Delta_{10}$ with the line joining its centre with the source. The angle $\Delta_{10}$ and radius $R_{0}$ are determined from conditions

[TABLE]

In an anisotropic chemical field, we find non-circular orbits of two distinct types - trajectories that resemble a figure of 8 and do not enclose the source and non-circular orbits that do (see fig. 5 (b) and brown regions of state diagram in fig. 2 (b)). These oscillations can occur in this inertia-less regime because of the dynamics of the extra degree of freedom, namely, the polar axis.

3.3 Effect of fluctuations on bound states

To check the stability of the bound states to thermal fluctuations, we add Gaussian white noise terms are added to the dynamics in eqs. 7, 8, 9 with the strength chosen as appropriate for a free swimmer of Peclet number $15$. The noisy trajectories so obtained and distributions of $R$, $\Delta_{1}$ and $\beta$ are shown in fig. 6. The values of $R\sigma^{-1}$ indicated on the axes are the mean values computed from the distribution and they match the deterministic values of $R_{0}$ shown in figs. 4 and 5.

3.4 Scattering

For large enough $h_{0}$, a chemotactic swimmer interacts briefly with the source before escaping it following paths that resemble scattering off an attractive centre, see fig. 7 (a). The scattering angle defined as $\Theta_{n1}=\lim_{t\to\infty}\theta_{n1}(t)$ is calculated by varying $h_{0}$ and the impact parameter $b_{0}$. $b_{0}$ is the initial lateral separation between the swimmer and the source. For a chemotactic swimmer at fixed $h_{0}$, there exists a threshold value of $b_{0}$ below which the swimmer is trapped in one of the two bound states. Above a maximum $h_{0}$, the swimmer is scattered for all $b_{0}$. Fig. 7 (b) shows $\Theta_{n1}$ for two different values of $h_{0}$. A swimmer that can form orbits revolves around the source several times before escaping which results in non-monotonic dependence of the scattering angle $\Theta_{n1}$ on $b_{0}$. An anti-chemotactic swimmer constantly turns away from the source and gets repelled, see fig. 7 (c), (d).

3.5 Fixed points in an isotropic source

To identify the two fixed points discussed so far, we will calculate the flow into these points (see fig. 8). We consider an isotropic source only, as the fixed points and phase portrait in this case can be exactly represented in 2D space spanned by $R$ and $\Delta_{1}$. Close to fixed point $(R_{0},\Delta_{10})$, the linearised equations can be written as $(\dot{\delta R},\dot{\delta\Delta_{1}})=\mathcal{M}\cdot(\delta R,\,\delta\Delta_{1})$. The signs of the real parts of the two eigenvalues of the dynamical matrix $\mathcal{M}$ determine system stability. At low $h_{0}$ the fixed point is stable node, see fig. 8 (a), characterised by two unequal negative real eigenvalues. With increasing $h_{0}$ the system undergoes a pitchfork bifurcation to a stable spiral, see fig. 8 (b), i.e. the two eigenvalues are complex conjugate of each other with negative real parts. On increasing $h_{0}$ further, the system undergoes a spiral-saddle bifurcation as the real part of the eigenvalues, with the eigenvector predominantly along $R$, turns positive causing $R$ to grow unbounded, see fig. 8 (h). At this value of $h_{0}$ transition for a bound state to a scattering. The other eigenvalue remains negative and $\Delta_{1}$ reaches the fixed value $\pi$ as seen in fig. 8 (i).

4 Two mobile swimmers

The motion of one swimmer in the body fixed frame of another, studied in section 3 and the features discussed can be used to study two mobile swimmers. We show that the final state is independent of the microscopic details and can be sorted into three broad classes: (1) Active dimers, whose relative distance stays fixed while the centre of mass translates uniformly, (2) Binary swimmers, whose polar axes become phase-locked, so that they revolve around a common centre in synchronised closed circular orbits while maintaining a fixed distance between their centres, like planetary orbits under gravity. (3) Scattering states, where the swimmers interact for a finite period of time and then take off at an angle to one another.

The interactions between the swimmers are not reciprocal. As a result they can be mutually attractive, mutually repulsive or show a mixed response; when one of the two is attracted to the other, while the other is repelled. The translational velocities of the swimmers 1 and 2 are $\mathbf{V}_{1}=V(\Delta_{1},\Delta_{2},R)\hat{\mathbf{R}}+\Omega(\Delta_{1},\Delta_{2},R)R\hat{\bm{\beta}},\,\,\mathbf{V}_{2}=-V(\Delta_{2},\Delta_{1},R)\hat{\mathbf{R}}+\Omega(\Delta_{2},\Delta_{1},R)R\hat{\bm{\beta}}$ and the angular velocities are $\dot{\theta}_{n1}=\omega(\Delta_{1},\Delta_{2},R),\,\,\dot{\theta}_{n2}=\omega(\Delta_{2},\Delta_{1},R)$. $R=|\mathbf{R}_{1}-\mathbf{R}_{2}|$ is the separation between the two swimmers and $\hat{\mathbf{R}}$ is the unit vector pointing from swimmer 2 to 1. The trajectories shown in the rest of this section are produced by solving the equations of motion with initial conditions $\theta_{n1}(0)-\theta_{n2}(0)=\pi$, $X_{1}(0)-X_{2}(0)=b_{0}$, where $b_{0}$ is the impact parameter, and $Z_{1}(0)-Z_{2}(0)=10^{3}\sigma$. We now discuss the three different types of dynamics in detail.

4.1 Active dimers

Two mutually attractive or mixed swimmers form a stable bound state where their relative separation freezes at a constant value $R_{0}$ and their polar axes align. We call this final state which can be non-moving or self-propelled, an active dimer. This state is seen when one of the two swimmers is designed so that it can be trapped by a source. Two mutually attractive swimmers form an active dimer when their polar axes point inwards along the line joining their centres and their relative velocity vanishes, so that $R_{0}$ satisfies: $V_{1}(0,0,R_{0})+V_{2}(0,0,R_{0})=0$. The dimer translates with the centre of mass velocity, defined in this inertia-less regime as the sum of the individual velocities as is given by $V_{1}(0,0,R_{0})-V_{2}(0,0,R_{0})$. A mixed pair also forms a dimer in a similar way, with the difference that one of them points towards the other while the other points away. The separation $R_{0}$ satisfies $V_{1}(0,0,R_{0})+V_{2}(0,\pi,R_{0})=0$ and the translation velocity is $V_{1}(0,0,R_{0})-V_{2}(0,\pi,R_{0})$, where 1 is attracted to 2, but 2 is repelled by 1. Bound states due to forces are formed when a particle sits in the minima of the effective potential produced by the other. In this system, long-ranged attractive and short-ranged repulsive drift velocities are generated by a balance between attraction due to chemotaxis and repulsion due to phoresis. This can be contrasted with other propelling clusters where the symmetry is broken by the shape of the cluster [9, 32, 23]. For two identical attractive swimmers, the condition to form a stationary dimer is $V(0,0,R_{0})=0$. The centre of mass of this system remains static always while the swimmers forms a stable dimer as ‘forces’ balance.

4.2 Binary swimmers

Two swimmers can go around a common point in synchronised circular orbits while maintaining fixed orientations $\Delta_{1}=\Delta_{10}$ and $\Delta_{2}=\Delta_{20}$, and a constant distance $R_{0}$ between their centres when at least one of them is designed such that it can orbit around a stationary source. We call this bound pair a binary swimmer for obvious reasons. The conditions of fixed relative distance and orientation are satisfied when

[TABLE]

The swimmers maintain fixed distances $R_{1,2}$ from their common centre. The total centre of mass velocity of the swimmer should be perpendicular to the line joining their centres to the common centre. This constraint and the equal orbital velocities of the swimmer provide three other conditions that are required to fully describe the bound state. We have the following relations for angles $\phi_{i0}$ between the sides of length $R_{0}$ and $R_{i}$

[TABLE]

Eqs. 14 and 15 can be solved in order to obtain the stationary angles $\Delta_{10}$, $\Delta_{20}$, $R_{0}$ and $R_{1,2}$. The swimmers orbit around a common centre with in circular orbits with angular frequency $\omega(\Delta_{10},\Delta_{20},R_{0})$. Fig. 10 shows cases where synchronous orbits are possible. For swimmer mobilities $m^{(1)}>m^{(2)}$, the frequency of orbital motion are different obtain loopy trajectories as seen in Figure 10 (c), (d).

For two identical swimmers, symmetry dictates that $\Delta_{10}=\Delta_{20}=\pi/2$. Thus identical binary swimmers are exactly like binary stars, as the component of their velocities parallel to the line joining the centre leads is zero and the non-zero perpendicular component leads to the rotation of the swimmers about a common point, see Figure 10 (b).

4.3 Effect of fluctuations on bound states

We end the discussion of bound states by testing the stability of these states to fluctuations. In order to do so additive brownian fluctuations are added to the equations of motion and distributions of $R$ and $\theta_{ni}$ are calculated in fig. 11. The distrbutions peak at vaues $R_{0}$ and $\Delta_{10}+\Delta_{20}$ as calulated for the noise-free dynamics.

4.4 Scattering

Two mutually attractive or mixed swimmers scatter apart if the combination of self-propulsion and interaction along their centres is enough to separate them, see 12 (a). Two anti-chemotactic swimmers continously anti-align and scatter away, see 12 (b). The active dimers and binary-swimmers discussed in the subsections 4.1 and 4.2 unbind for a large enough $b_{0}$, fig. 12 (d). The region of the state diagram in 2 where the swimmers scatter apart following trajectories typical of attractive or repulsive interactions are shown in white and blue respectively in fig. 2.

5 Conclusions

We have explored theoretically the varieties of dynamics exhibited by polar self-diffusiophoretic colloids, focusing on the case where the motion is planar. Gradients in the diffusing chemical species – generated by other colloids – affect the motion of the colloids, directly by translating them, and indirectly by rotating their polarity and thus their self-propelling velocities. These effects cooperate and compete, with remarkable dynamical consequences ranging from trapping, scattering and simple orbits around a fixed source to complex pairs of dances, which are the main results presented in this paper. The interactions are mediated by a diffusing field and are therefore naturally long-ranged. A distinguishing feature of our treatment is that the short-range interparticle repulsion is also dynamically generated by phoretic mechanisms. The orientation vector carried by each particle sets the direction of persistent motion. This endows the colloids with a kind of inertia. This is why, despite their vanishingly small Reynolds number, their behaviour bears some similarity to scattering in Newtonian mechanics. The analogy is incomplete as there is no equivalent of Newton’s III Law. It would appear therefore that there is no reduction of effective degrees of freedom in the 2-body problem; could it then show chaos? We do not see evidence for this in our numerical calculations. Possibly this is because the relaxational dynamics partially slaves some of the orientational degrees of freedom to each other and to the separation vectors; see eq. 7 and eq. 9. The case where the relaxation of chemical field occurs on timescales comparable to that of swimmer orientations, or or when the finite size of the swimmer creates orbits or active molecules [33, 34], would be an interesting extension, as would the case of non-planar, i.e., three-dimensional, motion. Meanwhile, our predictions for the planar case can be tested in suitable microfluidic setups [35].

A comparison to the Newtonian mechanics of particles interacting via central force fields is natural: the concentration of diffusing species plays a role similar to a potential, and the polar orientation resembles a momentum, endowing the dynamics with a character similar to inertia as it determines the direction of persistent motion (see [36] for another Stokesian driven system with an effective inertia). The absence of a conserved centre-of-mass momentum, however, renders the analogy incomplete, giving rise to the nontrivial joint motions of dimers discussed in section IV.

Hydrodynamics underlies the self-phoretic motility of our swimmers, but has not otherwise been included in our treatment. The absence of the hydrodynamic interaction in our treatment, even in the residual form of an incompressibility constraint, means the appropriate experimental test of our predictions should be under planar confinement between walls that are permeable to solvent and ions. In light of recent work [37] in which chemotaxis has been shown to compete against instabilities originating due to fluid flows, it is of importance to examine how the hydrodynamic interaction modifies the dynamics, especially of bound states, that we have described in this work. Whereas we have concentrated on the two-body problem, the implications of non-mutual pair interactions between dissimilar particles for the collective behaviour of active binary mixtures remain to be explored.

We expect our work to inspire efforts to fabricate particles with a range of catalytic and mobility coat patterns, using shape [38, 39] as an additional control parameter, This will allow a test of our predictions through exploration of our dynamical state diagrams.

SS thanks A Maitra, S Hablani, S Chatterjee, and K R Prathyusha for numerous discussions. SS thanks A Maitra particularly for discussions on fixed points and S Hablani for implementation of noisy dynamics. SS thanks R Seyboldt, A Maitra and T Adeleke-Larodo for detailed reading of the manuscript and offering many suggestions for improvement. SR was supported by a J C Bose Fellowship of the SERB, India, and a Homi Bhabha Chair Professorship of the Tata Education and Development Trust.

6 Appendix

The appendix is organised in three sections: in (A) we present details of the parameters chosen to construct the state diagrams in fig. 2, in (B) we outline the steps leading to the equations of motion in eqs. 7, 8 and 9 and in (C) the full equations of motion are written down.

6.1 Appendix A: Parameters for the state diagram

The state diagram is constructed by evaluating the trajectories for swimmers of different designs until a steady state is reached. Initial conditions are unchanged in evaluating a particular state diagram. Initial conditions for the single swimmer state diagrams in 2 (a), (b) and (c) are as follows: $\theta_{n1}(0)=\pi$ and $Z_{1}0)=10^{3}\sigma$, while $b_{0}=0.5\sigma$ for figs. 2 (a) and 2 (b) and $b_{0}=2\sigma$ for fig. 2 (c). Design parameters for single swimmer state diagrams: in fig. 2 (a) The source is isotropic with catalytic coat $\{1,0,0\}$. The swimmer catalytic coat design parameters are $\{1,-1,0\}$ and mobility $\mu_{0}^{(1)}=0.3$. $\mu_{1}^{(1)}$, $\mu_{2}^{(1)}$ are varied in the range $(-0.4,0.25)$ and $(-1,1)$ in steps of $0.025$. In fig. 2 (b) we use swimmer parameters same as in 2 (a), in presence of an anisotropic source with catalytic coat $\{1,1,0\}$. $\mu_{1}^{(1)}$, $\mu_{2}^{(1)}$ are varied in the range $(-0.4,0.25)$ and $(0.6,1.2)$ in steps of $0.01$. In fig. 2 (c) The swimmer and source design are the same as in fig. 2 (a); and $h_{0}=1$ is larger than the value $0.587$ for which we find bound states. The impact parameter is $b_{0}=2\sigma$. Note that only fig. 2 (c) is sensitive to choice of $b_{0}$. This is because given the choice of $h_{0}$, the capture radius for the swimmer is finite. If the impact parameter is increased larger than $30\sigma$, the bound state disappear completely.

The trajectories for two swimmers are evaluated for $\theta_{n1}-\theta_{n2}=2\pi$ and $Z_{1}-Z_{2}=2\times 10^{3}\sigma$. Parameters for two swimmer state diagrams: in fig. 2 (d), we consider two identical swimmers with $\mu^{(i)}_{0}=0.3$, $\mu^{(i)}_{1}=-0.5$ and $\mu^{(i)}_{2}$ is varied. The total catalytic activity $a^{(i)}$ is varied, the catalytic coat design is $\{1,-0.1,0\}$ and $b_{0}=2$. In 1(e) we consider two swimmers with mobility $\{0.3,-0.5,1\}$, catalytic coat $\{1,-1,0\}$, $b_{0}=4\sigma$ and surface activities $a^{(1)}$ and $a^{(2)}$ are varied. in fig. 2 (f) The total surface activity and $\mu^{(i)}_{2}$ is varied with other parameters being the same as in fig. 2 (e).

6.2 Appendix B: Evaluation of the chemical field

In this section we outline a method of reflections to calculate the chemical field $\rho$ produced by two swimmers of any orientation with a separation $R$. $\rho$ satisfies Neumann boundary condition on the swimmer surfaces, (eq. 2 in the main text).

[TABLE]

$\rho$ is a linear sum of contributions from both $A^{(i)}$, so it suffices to outline the calculation of $\rho$ for $a^{(1)}=0$. The solution can be split into two parts – $\rho={\rho}_{0}^{(2)}+\bar{\rho}$, where ${\rho}_{0}^{(2)}$ is the solution of the Poisson equation with the boundary condition eq. 17 on colloid 2 assuming colloid 1 to be absent. $\bar{\rho}$ is the correction due to presence of colloid 1 that we want to calculate. $\bar{\rho}$ is solved by a method of reflection, where at every order $n$, the boundary condition eq. 17 is enforced first on colloid 1 and then on 2. $\bar{\rho}$ is decomposed as follows

[TABLE]

The convexity of the spherical surface rules out multiple reflections at a given order $n$. ${\rho}^{(i)}_{n}$ satisfies the diffusion equation $\nabla^{2}{\rho}^{(i)}_{n}=0$ with boundary conditions

[TABLE]

where $n^{\prime}=n$ for $i=1$ and $n^{\prime}=n-1$ for $i=2$. Also note that $\rho_{0}^{(1)}=0$. $\nabla_{\perp}{\rho}^{(i)}_{n}(\mathbf{r}_{i})$ can be decomposed into surface harmonics as follows

[TABLE]

Substituting eq. 21 in eqs. 19 and 20, we derive the following recursion relation for coefficients $A^{(1)}_{n,\ell m}$.

[TABLE]

where $N^{12}_{j\ell m}$ is a geometrical factor that can be calculated as follows.

Recall from the main text, that the swimmers have radii $\sigma$, polar axes $\theta_{ni}$. In fig. 13, the frames of reference $\sum_{1}$ and $\sum_{2}$ are in the body fixed axes of the colloids and their $z$ axes are parallel to the line joining their centres. $\sum_{1}^{\prime}$ is the lab frame with its origin coinciding with that of $\sum_{1}$. We use spherical coordinates $(R_{1},\theta_{1},\phi_{1})$ and $(R_{2},\theta_{2},\phi_{2})$ in frames $\sum_{1}$ and $\sum_{2}$ respectively. The unit vectors representing polar axes are $(1,\theta_{n1},0)$ and $(1,\theta_{n2},0)$, see fig. 13. The expression for $\rho^{(1)}_{0}$ satisfying bounday conditions eq. 21 is

[TABLE]

The following transformations connect coordinates in $\sum_{1}$ and $\sum_{2}$

[TABLE]

Using eq. 24 in eq. 23, we obtain $\rho_{0}^{(1)}$ in $\sum_{1}$. Define

[TABLE]

A similar set of coefficients $N^{21}_{i\ell m}$ are defined by exchanging indices $1$ and $2$ in eq. 25. The recursion relations eq. 22 can be solved to obtain the field to any degree of precision. In order to obtain the equations of motions presented in the main text we truncate the expansion at $n=1$. For a uniformly coated colloid 2 with $\sigma^{(2)}_{0}=1$ we have the following result consistent with the potential due to a charge next to a sphere with Neumann boundary conditions (see [40])

[TABLE]

6.3 Appendix C: complete equations of motion

The angular momentum is

[TABLE]

The radial velocity is

[TABLE]

The orbital angular velocity is

[TABLE]