††thanks: Present Address: Department of Chemical and Materials Engineering,
University of Alberta, Edmonton, AB, T6G 1H9, Canada. Email:
<EMAIL_ADDRESS>
# Self-organization of active colloids mediated by chemical interactions

Zhiwei Peng  Raymond Kapral<EMAIL_ADDRESS>Chemical Physics Theory
Group, Department of Chemistry, University of Toronto, Toronto, ON, M5S 3H6,
Canada

###### Abstract

Self-propelled colloidal particles exhibit rich non-equilibrium phenomena and
have promising applications in fields such as drug delivery and self-assembled
active materials. Previous experimental and theoretical studies have shown
that chemically active colloids that consume or produce a chemical can self-
organize into clusters with diverse characteristics depending on the effective
phoretic interactions. In this paper, we investigate self-organization in
systems with multiple chemical species that undergo a network of reactions and
multiple colloidal species that participate in different reactions. Active
colloids propelled by complex chemical reactions with potentially nonlinear
kinetics can be realized using enzymatic reactions that occur on the surface
of enzyme-coated particles. To demonstrate how the self-organizing behavior
depends on the chemical reactions active colloids catalyze and their chemical
environment, we consider first a single type of colloid undergoing a simple
catalytic reaction, and compare this often-studied case with self-organization
in binary mixtures of colloids with sequential reactions, and binary mixtures
with nonlinear autocatalytic reactions. Our results show that in general
active colloids at low particle densities can form localized clusters in the
presence of bulk chemical reactions and phoretic attractions. The
characteristics of the clusters, however, depend on the reaction kinetics in
the bulk and on the particles and phoretic coefficients. With one or two
chemical species that only undergo surface reactions, the space for possible
self-organizations are limited. By considering the additional system
parameters that enter the chemical reaction network involving reactions on the
colloids and in the fluid, the design space of colloidal self-organization can
be enlarged, leading to a variety of non-equilibrium structures.

##  I Introduction

Self-organizing dynamics is often observed in active matter systems across a
wide range of length scales. Examples of self-organization in active matter
include aster formation of molecular motors and microtubules [1, 2, 3, 4],
swarming of bacteria [5, 6, 7], motility-induced phase separation of self-
propelled colloids [8, 9, 10, 11, 12, 13, 14, 15], and flocking behavior in
birds and fish [16, 17, 18, 19]. Through internal interactions, the
constitutive agents or particles of these systems are able to produce emergent
patterns, structures, or order. Among such interactions, chemical interaction
is often employed by both biological and synthetic active particles. Motile
microorganisms and cells can perform chemotaxis in response to chemical
concentration gradients [20, 21, 22]. In addition to responding to a
concentration gradient imposed in the environment, they may also produce or
consume the same chemical species, thus modifying the chemical concentration
in the environment. Synthetic diffusiophoretic colloids can move and reorient
in response to chemical gradients, mimicking the chemotactic behavior of
microorganisms.

Due to their ability to self-organize into dynamic structures, active colloids
are useful building blocks for self-assembled materials [23]. Previous
experimental and theoretical studies have shown that phoretic active colloids
can form clusters even at low particle densities [24, 25, 26, 27, 28, 29, 30,
31, 32]. Consider a self-diffusiophoretic particle that consumes a chemical in
a reaction on its surface and, in the process, creates a nonuniform chemical
concentration field. The chemical concentration gradients enable
diffusiophoresis in which particles may effectively attract or repel each
other depending on the details of the chemical interaction. If the particles
experience a diffusiophoretic force towards regions of lower concentration,
they tend to attract each other since they are consuming the chemical and act
as chemical sinks. Examples of the complex non-equilibrium dynamics systems of
this type exhibit, that include cluster and aster formation, as well as
plasma-like oscillations in response to perturbations, are given in Saha _et
al._ [33] who considered the chemotactic response self-diffusiophoretic
colloids to two chemical species $S$ and $P$ in the reaction $S\to P$.
Chemotactic aggregation or collapse has been observed in biological active
matter [34, 35, 36, 37]. When the translational and rotational phoretic motion
give competing attractive and repulsive interactions, dynamic clustering
states can be observed. Other variants of this model focusing on rotational
phoretic interaction also give rise to pattern formation and clustering [38,
39]; in this case, the active particles tend to align with or against the
gradient of the chemical field in addition to their self-propulsion. Continuum
theories have been developed that produce clustering and pattern formation
dynamics similar to the particle-based models [27, 33, 38, 39, 40].

So far we have outlined the dynamics of a single type of colloid interacting
via a single chemical species. That is, all colloids in the system respond to
the chemical field and produce or consume the chemical in the same fashion. A
natural extension of these previous studies is to consider the dynamics of
mixtures of colloids that interact with the single chemical species
differently. In a binary mixture of such colloids, the effective action-
reaction symmetry can be broken, which leads to non-reciprocal interactions
between the two types of colloids. Indeed, previous work has shown that such
mixtures exhibit aggregation or phase-separated states with distinct densities
depending on the interaction parameters [41].

In this paper, we generalize previous studies by considering a model for
active colloidal mixtures that interact chemically with multiple chemical
species that participate in a network of chemical reactions. Instead of
producing or consuming a single chemical, the colloids may participate in
multiple chemical reactions, thus potentially producing or consuming multiple
chemical species. By allowing the particles and fluid to participate in
multiple coupled reactions, the design space for colloidal self-organization
and self-assembly is greatly enlarged. Active colloids propelled by complex
chemical reactions with potentially nonlinear kinetics can be realized using
enzymatic reactions that occur on the surface of enzyme-coated particles [42,
43, 44, 45, 46, 47, 48, 49]. When an active particle participates in several
chemical reactions, its self-propulsion may be fueled by multiple chemical
species. Self-propelled active particles driven by multiple fuel species have
been realized in the laboratory [50, 51]. In contrast to active particles
powered by the decomposition of hydrogen peroxide, enzyme-powered particles or
motors are more desirable for biomedical applications due to their
nontoxicity, versatility, and biofuel availability [52, 53, 54, 55].

In Section II a general formulation of the chemical reaction network and the
Langevin equations of motion are presented in which both bulk reactions and
surface reactions on the colloids are considered. We model the surface
reaction as delta function sources or sinks but with general reaction
kinetics. Bulk reactions are also included, which act to maintain the chemical
system at a non-equilibrium state. We then present three case studies of
active colloids undergoing different reactions. In Section III, we show the
clustering dynamics of a single type of colloid undergoing a single surface
reaction. Because the bulk reaction is present, the chemical concentrations
are screened and we observe localized clusters with phoretic attraction. We
consider binary mixtures with different reaction kinetics in Sections IV and
V. In Section IV, a sequential reaction is considered in which the product of
the first reaction is the reactant of the second reaction in the sequence, and
the first reaction occurs on the first type of colloid while the second
reaction is on the second type of colloid. We show that this coupling induced
by the sequential reaction gives rise to diverse clustering behavior. In
Section V a binary mixture with a nonlinear autocatalytic reaction is
considered before concluding the paper in Section VI.

##  II  General formulation

Consider a suspension of chemically active colloidal particles dispersed in an
incompressible Newtonian fluid (solvent) of dynamic viscosity $\eta$. In
addition to the $N$ active colloids, the suspension contains $n_{s}$ reacting
solute species, $\bm{S}=\\{S_{j}\mid j=1,\dots,n_{s}\\}$. The local
concentrations of $\bm{S}$ at position $\bm{x}$ and time $t$ are denoted as
$\bm{c}(\bm{x},t)=\\{c_{j}(\bm{x},t)\mid j=1,\dots,n_{s}\\}$. The $n_{s}$
solute species form a chemical reaction network in which the solute $S_{j}$
can act as either reactants ($R_{j}$) or products ($P_{j}$) in potentially
multiple chemical reactions. Possibly reversible chemical reactions, which
were chosen from a set of $n_{r}$ such reaction schemes, occur on the surface
of the $N$ active colloids. For a particular surface reaction indexed by $K$
($K=1,\dots,n_{r}$), we may write the reaction scheme compactly as
$\bm{R}^{K}\rightleftharpoons\bm{P}^{K}$. Here, the set of reactants is a
subset of $\bm{S}$, $\bm{R}^{K}=\\{S_{j}\lvert j\in I^{K}\\}\subset\bm{S}$,
where $I^{K}$ is the index set that contains the chemical species acting as
reactants for the reaction indexed by $K$. Similarly, the set of products of
reaction $K$ is $\bm{P}^{K}=\\{S_{j}\lvert j\in O^{K}\\}\subset\bm{S}$, where
$O^{K}$ is the index set that contains the chemical species acting as products
for the reaction indexed by $K$. The local concentrations corresponding to
these sets of species are similarly defined as $\bm{c}_{R}^{K}(\bm{x},t)$ and
$\bm{c}_{P}^{K}(\bm{x},t)$.

The local chemical concentration fields $\bm{c}(\bm{x},t)$ evolve by the
reaction-diffusion equations that account for reactions on the surfaces of the
active particles as well as reactions in the fluid phase that both participate
in the reaction network and serve to maintain the system in a nonequilibrium
state so that active motion is possible. The reaction-diffusion equation is
given by

$\displaystyle\partial_{t}c_{j}(\bm{x},t)=$ $\displaystyle
D_{j}\nabla^{2}c_{j}-\sum_{K=1}^{n_{r}}\mathcal{R}^{b}_{j,K}(\bm{c}^{K})$
$\displaystyle+\sum_{K=1}^{n_{r}}\sum_{i_{K}=1}^{N_{K}}\mathcal{R}^{s}_{j,K}(\bm{c}^{K})\delta(\bm{x}-\bm{X}_{i_{K}}),$
(1)

where $\bm{c}^{K}$ includes the concentrations of both the reactants
$\bm{c}_{R}^{K}(\bm{x},t)$ and the products $\bm{c}_{P}^{K}(\bm{x},t)$. The
second term on the right side accounts for reactions in the bulk fluid phase
that may help to maintain the nonequilibrium state of the system. The
nonequilibrium state could also be established through chemostats at the
boundaries. The third term accounts for surface reactions on the active
particles and in writing this term we adopted a monopole approximation for the
surface reactions, following previous studies [27, 28, 38, 56, 39]. The
positions of the particle centers are denoted by $\bm{X}_{i_{K}}$. Note that
the reaction indexed by $K$ occurs on the particle labeled $i_{K}$, and there
are in total $N_{K}$ such particles. As a result, in the third term on the
right side of (II), we first sum over all particles ($N_{K}$) that carry the
reaction indexed by $K$ and then sum over all $n_{r}$ reactions. In general,
both the bulk and surface rates, $\mathcal{R}^{b}_{j,K}$ and
$\mathcal{R}^{s}_{j,K}$, are nonlinear functions of the chemical
concentrations, and both forward and reverse rates are included in the case of
reversible reactions. In (II), if $c_{j}$ does not participate in the reaction
indexed by $K$, the rate coefficients $\mathcal{R}^{b}_{j,K}$ and
$\mathcal{R}^{b}_{j,K}$ are understood to be zero. For the remainder of this
paper, we assume that all chemical species have the same molecular diffusivity
$D$.

As a consequence of microscopic reversibility, the surface reaction rates also
depend on the forces the colloids experience [57, 58, 59], including external
forces as well as colloid-colloid interactions arising from chemical gradients
and direct interactions. The magnitude and sign of the effect depends strongly
on the orientation of the Janus colloid relative to the force. In our
simulations neither the forces a motor experiences nor its orientation vector
are controlled; hence, in what follows, we may neglect this effect on the
reaction rate (see ESI). We note that an external force applied to a
collection of Janus colloids with controlled orientation vectors has been
shown to change the net reaction rate of the collection. [40]

We consider monodisperse and spherical active particles of radius $a$. The
linear equations of motion for the active particles may be written as

$\displaystyle\frac{\text{d}}{\text{d}t}\bm{X}_{i_{K}}=V_{\rm
sd}\big{(}\bm{c}^{K}\big{)}\bm{u}_{i_{K}}+\sum_{j=1}^{n_{s}}\mu_{i_{K}j}\bm{\nabla}c_{j}+\frac{\bm{F}_{i_{K}}^{P}}{\zeta_{t}}+\bm{V}^{B}_{i_{K}}.$
(2)

The first term on the right is the self-diffusiophoretic velocity directed
along $\bm{u}_{i_{K}}$ that depends on the concentrations in the $K$ reaction,
while the second term is the diffusiophoretic velocity due to the
concentration gradients in the surroundings. The coefficient $\mu_{i_{K}j}$
represents the response of particle $i_{K}$ to the chemical species $j$;
$\mu_{i_{K}j}$ may be zero for some species and take either sign ($+$ or $-$)
depending on the physical context. In the third term $\bm{F}_{i_{K}}^{P}$ is
the impulsive force on particle $i_{K}$ due to hard interactions with the
other active particles, and $\zeta_{t}$ is the translational friction
coefficient of the particles. The last term is the Gaussian, fluctuating
Brownian velocity, which satisfies

$\displaystyle\langle\bm{V}^{B}_{i_{K}}(t)\rangle=\bm{0}\quad{\rm
and}\quad\langle\bm{V}^{B}_{i_{K}}(t)\bm{V}^{B}_{i_{K}}(t^{\prime})\rangle=2D_{t}\delta(t-t^{\prime})\bm{I},$
(3)

where $D_{t}$ is the translational diffusivity, and $\bm{I}$ is the identity
tensor. In writing (2), hydrodynamic interactions among particles are
neglected and the translational friction coefficient is taken to be
$\zeta_{t}=6\pi\eta a$. The translational diffusivity is related to the
thermal energy $k_{B}T$ via the relation $k_{B}T=\zeta_{t}D_{t}$. We note that
previous studies suggest that chemotactic forces may dominate the clustering
dynamics while hydrodynamic interactions play a less important role [60].

The angular equations of motion for the active particles are given by

$\displaystyle\frac{\text{d}}{\text{d}t}\bm{u}_{i_{K}}=\sum_{j=1}^{n_{s}}\gamma_{i_{K}j}\left(\bm{u}_{i_{K}}\times\bm{\nabla}c_{j}\right)\times\bm{u}_{i_{K}}+\bm{\Omega}^{B}_{i_{K}}\times\bm{u}_{i_{K}},$
(4)

where the signed constant $\gamma_{i_{K}j}$ is the angular diffusiophoretic
coefficient for species $j$, and $\bm{\Omega}^{B}_{i_{K}}$ is the fluctuating
Brownian angular velocity given by

$\displaystyle\langle\bm{\Omega}^{B}_{i_{K}}(t)\rangle=\bm{0}\quad{\rm
and}\quad\langle\bm{\Omega}^{B}_{i_{K}}(t)\bm{\Omega}^{B}_{i_{K}}(t^{\prime})\rangle=2D_{r}\delta(t-t^{\prime})\bm{I},$
(5)

where $D_{r}$ is the rotational diffusivity that satisfies
$\zeta_{r}D_{r}=k_{B}T$, and $\zeta_{r}=8\pi\eta a^{3}$ is the rotational
friction coefficient. In (4), notice that
$\left(\bm{u}_{i_{K}}\times\bm{\nabla}c_{j}\right)\times\bm{u}_{i_{K}}=\left(\bm{I}-\bm{u}_{i_{K}}\bm{u}_{i_{K}}\right)\cdot\bm{\nabla}c_{j}$.

Equations (2) and (4) are a generalization of the equations considered by
previous studies that include one chemical species [27, 28]. Modeling the
chemical as a diffusing species in the presence of delta function sinks
located at the particle centers, Pohl and Stark [27, 28] showed that the
colloids can establish dynamic clustering states or collapse into a single
cluster, depending on the phoretic parameters. In this model, the chemical
concentration at a distance of $r$ away from a sink scales as $1/r$ (in 3D).
As a result, if the phoretic interaction is purely attractive, this long-
ranged interaction (the chemical gradient scales as $1/r^{2}$) ultimately
leads to the aggregation of colloids into a single large cluster. In the
current work, our aim is to demonstrate that complex kinetics involving
multiple chemical species can lead to more diverse clustering dynamics, which
may be useful for the design of self-assembling active materials.

##  III Active colloids propelled by a single surface reaction

We now consider the simplest case in which all particles participate in the
same surface reaction, given by $S_{1}\to S_{2}$. In the bulk, $S_{1}$ is
replenished using the effective reaction $S_{2}\to S_{1}$ so as to maintain a
non-equilibrium state. We remark that the bulk reaction should not be treated
as the reverse of the surface reaction but as a separate reaction with other
participating chemical species. These additional chemical species are held at
fixed concentrations by chemical reservoirs [61]. As a result, the
concentration of the reservoir species do not appear in the effective bulk
reaction. This simple reaction scheme serves as an introduction to the general
mechanisms of phoretic interactions and to the effects of chemical screening
due to bulk reactions.

For this reaction scheme, the reaction-diffusion equation (II) becomes

$\displaystyle\frac{\partial c_{1}}{\partial t}$
$\displaystyle=D\nabla^{2}c_{1}-\kappa\sum_{i=1}^{N}c_{1}\delta(\bm{x}-\bm{X}_{i})+k_{b}c_{2},$
(6a) $\displaystyle\frac{\partial c_{2}}{\partial t}$
$\displaystyle=D\nabla^{2}c_{2}+\kappa\sum_{i=1}^{N}c_{1}\delta(\bm{x}-\bm{X}_{i})-k_{b}c_{2},$
(6b)

where $\kappa$ is an effective reaction rate for the surface reaction and
$k_{b}$ is the bulk reaction rate. We note that in 3D, $\kappa$ is the surface
rate coefficient (which has units of length/time) multiplied by the surface
area; in 2D, it is multiplied by the circumference. Because both species have
the same diffusivity, we obtain $\partial(c_{1}+c_{2})/\partial t\equiv 0$.
Suppose that initially $c_{1}+c_{2}=c_{0}=const.$, then the total
concentration $c_{0}$ remains constant in time and homogeneous in space. Upon
substituting $c_{2}=c_{0}-c_{1}$ into (6), one only needs to solve (6a). The
equations of motion (2) and (4) are now written as

$\displaystyle\frac{\text{d}}{\text{d}t}\bm{X}_{i}$
$\displaystyle=\alpha\,c_{1}\,\bm{u}_{i}+\mu\,\bm{\nabla}c_{1}+\frac{\bm{F}_{i}^{P}}{\zeta_{t}}+\bm{V}_{i}^{B},$
(7a) $\displaystyle\frac{\text{d}}{\text{d}t}\bm{u}_{i}$
$\displaystyle=\gamma\left(\bm{u}_{i}\times\bm{\nabla}c_{1}\right)\times\bm{u}_{i}+\bm{\Omega}^{B}_{i}\times\bm{u}_{i},$
(7b)

where the self-diffusiophoretic speed is taken to be linearly proportional to
the local concentration, $V_{\rm sd}=\alpha c_{1}$. This form of a linear
self-diffusiophoretic speed was considered in Ref. [33]. We have also assumed
that all particles respond to the chemical gradient $\bm{\nabla}c_{1}$ in the
same fashion and the subscripts of $\mu$ and $\gamma$ in (2) and (4) are
dropped.

In the presence of bulk “refueling” (i.e., $S_{2}\to S_{1}$), the system (6)
is chemically screened; the inverse screening length is

$\displaystyle\lambda=\sqrt{\frac{k_{b}}{D}}.$ (8)

In the remainder of the paper, we focus on 2D systems with monodisperse disks
of radius $a$. Equations (6) and (7) are governed by the following
dimensionless groups:

$\displaystyle\Lambda=\lambda a,\quad Da=\frac{\kappa}{4\pi
aD},\quad\hat{\alpha}=\frac{\alpha c_{0}a}{D},$ (9a)
$\displaystyle\quad\hat{\mu}=\frac{\mu
c_{0}}{D},\quad\hat{\gamma}=\frac{\gamma c_{0}a}{D},\quad\frac{D_{t}}{D}.$
(9b)

Here, $Da$ is a Damköhler number that compares the surface reaction rate to
the rate of diffusion. We note that $a^{2}D_{r}/D_{t}=3/4$ for isolated
spheres of radius $a$. Since the chemical species often diffuse faster than
the colloids, in the remainder of the paper, we fix $D_{t}/D$ to be $0.1$ and
focus on the variation of the other dimensionless parameters. Finally, the
area fraction of the particles is important and denoted as $\phi$.

Figure 1: Schematic of effective attractions between two active particles
mediated by chemical interactions. The red arrows surrounding the blue
particle indicate the concentration gradient $\bm{\nabla}c_{1}$. (I): For
$\mu<0$, the gray particle moves towards the blue particle due to the
diffusiophoretic velocity in the presence of the gradient,
$\mu\bm{\nabla}c_{1}$. (II): For $\gamma<0$ and $\alpha>0$, the gray particles
are first rotated ($\gamma\bm{u}\times\bm{\nabla}c_{1}$) and then can swim
towards the blue particle with their self-diffusiophoretic velocity, $\alpha
c_{1}\bm{u}$.

The single surface reaction model specified by (6) and (7) can be treated as a
generalization of previous models [27, 28, 39]. More specifically, the swim
speed (the first term on the right side of (7a)) in these previous models is
taken to be a constant and does not explicitly depend on the local chemical
concentration $c_{1}$. Furthermore, in Pohl and Stark [27, 28], the bulk
refueling is absent; in Liebchen _et al._ [39], $\mu\equiv 0$, and the active
particles do not exhibit translational diffusiophoretic motion as a result of
the gradient $\bm{\nabla}c_{1}$ induced by other particles.

Figure 2: Representative snapshots of active particles at area fraction
$\phi=0.1$. The red background is the concentration field of the product
species $c_{2}$. (a): Snapshot of active Brownian particles,
$Da,\hat{\mu},\hat{\gamma}=0$ and $\hat{\alpha}=10$. In this case, the
concentration remains a constant in space and time, $c_{1}\equiv c_{0}$ and
$c_{2}\equiv 0$, because initially $c_{1}=c_{0}$. (b): Clustering of active
particles due to phoretic attraction (mechanism I in Fig. 1),
$\hat{\alpha},\hat{\gamma}=0;\hat{\mu}=-10$. (c): Homogeneous state in the
presence of phoretic repulsion,$\hat{\alpha},\hat{\gamma}=0;\hat{\mu}=10$.
(d): Clustering of active particles due to phoretic attraction (mechanism II
in Fig. 1), $\hat{\alpha}=10,\hat{\mu}=0,\hat{\gamma}=-10$. (e): Homogeneous
state in the presence of phoretic repulsion,
$\hat{\alpha}=10,\hat{\mu}=0,\hat{\gamma}=10$. (f): Clustering of active
particles in the presence of both phoretic attraction mechanisms,
$\hat{\alpha}=10,\hat{\mu}=-10,\hat{\gamma}=-10$. For all panels,
$\Lambda=0.25$. For (b)-(f), $Da=1/(2\pi)\approx 0.16$.

In our colloid model, which is similar to that of Saha, et al. [33], there are
two phoretic mechanisms that lead to effective attraction between particles
(see Fig. 1). The first phoretic attraction mechanism originates from the
phoretic translation of a particle (gray) in the presence of the chemical
gradient induced by another active particle (blue). As shown in Fig. 1, the
surface reaction on the blue particle induces a radial concentration gradient
in the vicinity. For $\mu<0$, the concentration gradient will induce a
translational velocity, $\mu\bm{\nabla}c_{1}$, that directs the second
particle nearby towards the blue particle. We note that for $\mu>0$, the
effective interaction becomes repulsive. The second mechanism for attraction
results from the combined effect of self-diffusiophoretic swimming, $\alpha
c_{1}\bm{u}$, and phoretic alignment, $\gamma\bm{u}\times\bm{\nabla}c_{1}$. In
the radial concentration gradient, the gray particles experience a
diffusiophoretic torque that tends to align the swimming direction $\bm{u}$
against the gradient $\bm{\nabla}c_{1}$ provided that $\gamma<0$. Once
aligned, the gray particle can swim towards the blue particle via self-
diffusiophoretic translational motion if $\alpha>0$. Again, if $\gamma>0$, the
effective interaction becomes repulsive. Because of the conservation,
$c_{0}=c_{1}+c_{2}$, we have $\bm{\nabla}c_{1}=-\bm{\nabla}c_{2}$; one can
then reformulate (7) in terms of $c_{2}$, in which case the sign of $\mu$ and
$\gamma$ needs to be reversed in order to have attractive interactions.

Typical snapshots of the system at long times with an area fraction $\phi=0.1$
are shown in Fig. 2. In the absence of surface reaction ($Da=0$), the chemical
concentration remains a constant in space and time if initially
$c_{1}=c_{0}=const.$ As a result, the system reduces to the bare dynamics of
active Brownian particles with hard-sphere interactions; the configuration of
the system is shown in Fig. 2(a). At this low density, motility-induced phase
separation [8, 9, 10, 11, 12, 13, 14] is absent and the system remains in a
homogeneous state with visible density fluctuations. In Fig. 2(b), the self-
diffusiophoretic swimming is turned off ($\alpha=0$), and particles form
localized clusters due to phoretic attraction induced by the chemical gradient
(mechanism I in Fig. 1). If one reverses the sign of $\mu$ ($-$ to $+$), the
effective interaction becomes repulsive, and particles remain in a homogeneous
state (see Fig. 2(c)).

In Fig. 2(d), $\mu=0$ and $\gamma<0$, particles form clusters due to the
second mechanism of phoretic attraction (see Fig. 1). In contrast to Fig.
2(a), we see that the rotation and alignment from phoretic interactions allow
one to turn the random distribution of active Brownian particles into
localized clusters. If one reverses the sign of $\gamma$ from negative to
positive, particles experience effective repulsion and do not form clusters
(see Fig. 2(e)). Lastly, in Fig. 2(f), both mechanisms of attraction are
present and one observes cluster formation. Because the chemical system is
screened, the clusters are localized and one does not observe the aggregation
of all particles into a single cluster (cf. Pohl and Stark [27]). We note that
with screening large-scale clusters can be obtained for higher area fractions.

##  IV Binary mixture with sequential reactions

Figure 3: Schematic of the sequential reaction scheme. The surface reaction
$S_{1}\to S_{2}$ occurs on the first type of particles while $S_{2}\to S_{3}$
occurs on the second type.

We now turn our attention to binary mixtures. The system of active particles
is classified into two types. The surface reaction $S_{1}\to S_{2}$ occurs on
the surface of type 1 particles. For type 2 particles, we have the surface
reaction $S_{2}\to S_{3}$. Similar to the previous section, bulk reactions
$S_{3}\to S_{2}$ and $S_{2}\to S_{1}$ (see Fig. 3 for a schematic) are
incorporated to maintain the system at a nonequilibrium state (see Fig. 3).
Similar to the reactions considered in the previous section, the bulk
reactions are effective reactions that are maintained out of equilibrium using
chemical reservoirs of other species. The reaction-diffusion system governing
the concentration $c_{i}$, $i=1,2,3$, is given by

$\displaystyle\frac{\partial c_{1}}{\partial t}=$ $\displaystyle
D\nabla^{2}c_{1}-\kappa_{1}\sum_{i\in\sigma_{1}}c_{1}\delta\left(\bm{x}-\bm{X}_{i}\right)+k_{b,1}c_{2},$
(10a) $\displaystyle\frac{\partial c_{2}}{\partial t}=$ $\displaystyle
D\nabla^{2}c_{2}+\kappa_{1}\sum_{i\in\sigma_{1}}c_{1}\delta\left(\bm{x}-\bm{X}_{i}\right)-k_{b,1}c_{2}$
$\displaystyle-\kappa_{2}\sum_{j\in\sigma_{2}}c_{2}\delta\left(\bm{x}-\bm{X}_{j}\right)+k_{b,2}c_{3},$
(10b) $\displaystyle\frac{\partial c_{3}}{\partial t}=$ $\displaystyle
D\nabla^{2}c_{3}+\kappa_{2}\sum_{j\in\sigma_{2}}c_{2}\delta\left(\bm{x}-\bm{X}_{j}\right)-k_{b,2}c_{3}.$
(10c)

In the above, $\sigma_{1}=\\{i\,\lvert\,i\in\mathbb{Z},1\leq i\leq N_{1}\\}$
denotes the index set of $N_{1}$ type 1 particles. There are $N-N_{1}$ type 2
particles denoted by the set
$\sigma_{2}=\\{i\,\lvert\,i\in\mathbb{Z},N_{1}<i\leq N\\}$. The rate
coefficients, $\kappa_{1}$ and $\kappa_{2}$, are the surface rate of the first
and second reaction, respectively; $k_{b,1}$ and $k_{b,2}$ are the bulk
reaction rate coefficients of the reactions $S_{2}\to S_{1}$ and $S_{3}\to
S_{2}$, respectively. Adding (10a)–(10c), we see that the chemical species
satisfy the conservation condition, $c_{1}+c_{2}+c_{3}=c_{0}=const.$ The
active velocity for the type $1$ particle with index $i$ is given by
$\alpha_{1}c_{1}\bm{u}_{i}+\mu_{11}\bm{\nabla}c_{1}+\mu_{12}\bm{\nabla}c_{2}$;
for the type 2 particle with index $j$, we have the active velocity given by
$\alpha_{2}c_{2}\bm{u}_{j}+\mu_{21}\bm{\nabla}c_{1}+\mu_{22}\bm{\nabla}c_{2}$.
Because $\bm{\nabla}c_{3}=-\left(\bm{\nabla}c_{1}+\bm{\nabla}c_{2}\right)$,
the term $\mu_{13}\bm{\nabla}c_{3}$ is not needed and has been absorbed into
the coefficients $\mu_{11}$ and $\mu_{12}$. The active velocity for type 2
particles is treated similarly. The active angular velocities for type 1 and 2
are, respectively,
$\bm{u}_{i}\times\left(\gamma_{11}\bm{\nabla}c_{1}+\gamma_{12}\bm{\nabla}c_{2}\right)$
and
$\bm{u}_{j}\times\left(\gamma_{21}\bm{\nabla}c_{1}+\gamma_{22}\bm{\nabla}c_{2}\right)$.
The fluctuating and hard-sphere terms are the same as those given in (7).

Figure 4: Representative snapshots of binary mixtures of active particles
undergoing sequential surface reactions at total area fraction $\phi=0.1$.
Type 1 particles are shown in blue and type 2 are orange. The background is
the concentration field of the intermediate species $c_{2}$. For all
simulations, the particles are initialized with random positions and
orientations, the initial concentrations are $c_{1}=c_{0}$, $c_{2}=c_{3}=0$,
$\Lambda_{1}=\Lambda_{2}=0.25$, $Da_{1}=Da_{2}=1/(2\pi)\approx 0.16$,
$\hat{\mu}_{12}=\hat{\mu}_{21}=0$, $\hat{\gamma}_{12}=\hat{\gamma}_{21}=0$,
and the number of type 1 particles are the same as 2. For
(a)–(c),$\hat{\alpha}_{1}=\hat{\alpha}_{2}=\hat{\gamma}_{11}=\hat{\gamma}_{22}=0$;
for (d)–(f), $\hat{\mu}_{11}=\hat{\mu}_{22}=0$, and
$\hat{\alpha}_{1}=\hat{\alpha}_{2}=10$. (a): Type 1 particles form clusters
while type 2 are dispersed, $\hat{\mu}_{11}=\hat{\mu}_{22}=-20$. (b):
Heterogeneous clusters are observed that contain both types of particles,
$\hat{\mu}_{11}=-20,\hat{\mu}_{22}=20$. (c): Both 1 and 2 are dispersed,
$\hat{\mu}_{11}=\hat{\mu}_{22}=20$. (d): Type 1 particles form clusters while
type 2 are dispersed, $\hat{\gamma}_{11}=\hat{\gamma}_{22}=-20$. (e):
Heterogeneous clusters are observed that contain both types of particles,
$\hat{\gamma}_{11}=-20,\hat{\gamma}_{22}=20$. (f): Both 1 and 2 are dispersed,
$\hat{\gamma}_{11}=\hat{\gamma}_{22}=20$.

In this binary model, the chemical coupling between the dynamics of the two
types of particles is through the intermediate species $S_{2}$. The equations
are non-dimensionalized following the scheme given in (9). In contrast to the
single reaction case considered in the previous section, we now have an even
larger set of dimensionless parameters. For simplicity, we consider the case
in which the cross terms $\mu_{12},\mu_{21},\gamma_{12}$ and $\gamma_{21}$ are
zero. We focus on the effects of $\mu_{11},\mu_{22},\gamma_{11}$ and
$\gamma_{22}$ on the clustering behavior.

In Fig. 4, we show typical simulation snapshots of the system at long times
with a total (all particles) area fraction $\phi=0.1$. The blue particles in
Fig. 4 denote type 1 while the orange are type 2. When $\mu_{11}<0$, type 1
particles tend to aggregate because they consume $S_{1}$ (sinks of $S_{1}$),
which induces a phoretic attraction. This aggregation of type 1 particles is
shown in Fig. 4(a) and (b). In Fig. 4(a), $\mu_{22}<0$, type 2 particles are
repelled from sources of $S_{2}$, i.e., type 1 particles. Since $\mu_{22}<0$,
type 2 particles are also attracted to sinks of $S_{2}$, or type 2 particles.
Because $S_{2}$ is only produced by type 1 particles (and in the fluid from
$S_{3}$) and the reaction is screened, the concentration of $S_{2}$ ($c_{2}$)
is very small at locations of type 2 particles (see the background contour in
Fig. 4(a)). In addition, the bulk reaction $S_{2}\to S_{1}$ diminishes the
concentration of $S_{2}$ near type 2 particles. This, in turn, reduces the
production of $S_{3}$ (on type 2 particles) and hence $S_{2}$ by the bulk
reaction $S_{3}\to S_{2}$. As a result, the attraction among type 2 particles
is weak and one observes that they remain loosely dispersed away from islands
of type 1 particles. Effectively, in Fig. 4(a), the binary system exhibits a
phoretic attraction among type 1 particles and type 2 particles are repelled
from type 1. This phoretic repulsion constitutes a non-reciprocal chemical
interaction because type 2 are repelled from type 1, but not vice versa. For
weaker screening, or larger screening length, the size of type 1 particle
clusters is larger. In this case, we observe that type 2 particles are pushed
away from type 1 clusters and appear to be more closely aggregated (see ESI,
Fig. 1). As multiple type 1 clusters move towards each other due to longer
range attraction, they may squeeze a collection of type 2 particles together
and enclose them as part of the merged cluster (see Fig. 5).

Figure 5: Illustration of the 1-1 attraction induced trapping of type 2
particles. Left: the blue clusters are attracted towards each other, their
motion are indicated by the black arrows. Right: as the blue clusters
aggregate, some gold particles escape from the two open channels shown on the
left while the rest are enclosed by the large cluster. The schematic is based
on the simulation shown in Fig. 9 of the ESI.

When the sign of $\mu_{22}$ is switched from negative to positive (Fig. 4(a)
$\to$ Fig. 4(b)), type 1 particles become attractive to type 2 and we observe
heterogeneous clusters consisting of both types of particles. The clustering
dynamics in Fig. 4(b) is subtle and warrants a more detailed consideration. In
particular, we notice that type 2 particles are often on the periphery of the
clusters while type 1 particles occupy the center. Because type 1 particles
consumes $S_{1}$ and $\mu_{11}<0$, they tend to move towards particles of the
same type and form clusters of type 1. Once a type 1 cluster of several
particles is formed, a noticeable density gradient of $S_{2}$ is established
locally. Since $\mu_{22}>0$, type 2 particles tend to move towards sources of
$S_{2}$, which are the localized clusters of type 1 particles. As a result,
type 2 particles are attracted to type 1 from the term
$\mu_{22}\bm{\nabla}c_{2}$ but type 1 particles are not attracted to type 2.
This non-reciprocal attraction ultimately leads to clusters in which type 2
particles are attached to the periphery. We also note that type 2 particles do
not form clusters on their own because they consume $S_{2}$ but are attracted
to sources of $S_{2}$.

In Fig. 4(c), $\mu_{11}>0$, which implies that type 1 particles are repelled
from sinks of $S_{1}$. In other words, they exhibit effective repulsion and do
not form clusters or establish a noticeable concentration gradient of $S_{2}$.
In the absence of $\bm{\nabla}c_{2}$, type 2 particles will also become
dispersed since there are no type 1 clusters to which they are attracted. As a
result, we see that both types of particles are dispersed.

For Figs. 4(a)–(c), the phoretic interactions result from the type I mechanism
outlined in Fig. (1), though with more than one chemical species. Just like
the single reaction case considered in the previous section, one can achieve
mixture clustering dynamics similar to those shown in Figs. 4(a)–(c) using the
type II mechanism. These results are shown in Figs. 4(d)–(f), where the type I
mechanism is turned off. That is, $\mu_{ij}=0$, $i,j=1,2$, and
$\alpha_{1},\alpha_{2},\gamma_{11}$ and $\gamma_{22}$ are nonzero. In this
case, particles self-propel with their self-diffusiophoretic velocities and
are rotated to align with or against density gradients depending on the sign
of $\gamma_{ij}$. In Fig. 4(d), $\gamma_{11},\gamma_{22}<0$, and we observe
clustering of type 1 particles as a result of their phoretic attraction;
similar to Fig. 4(a), type 2 particles are dispersed in the bulk. In Fig.
4(e), $\gamma_{11}<0$ and $\gamma_{22}>0$, type 2 particles are attracted to
type 1 clusters. When $\gamma_{11},\gamma_{22}>0$, both type 1 and type 2 are
dispersed, as shown in 4(f).

In Fig. 4(b) and (e), we see that there are type 2 particles trapped in the
middle of the clusters. This occurs because two small type 1 clusters with
type 2 particles attached to their periphery can merge into a single cluster
as a result of their phoretic attraction, which leaves some type 2 particles
trapped in the middle of the new cluster. To show this merging process, in
Fig. 7 we plot the system snapshots at different times. In contrast to the
simulation of Fig. (4)(b), the phoretic coefficient $\hat{\mu}_{22}$ in the
simulation of Fig. 7 has a delayed activation. That is, initially
$\hat{\mu}_{22}=0$, once the system establishes type 1 clusters (Fig. 7(b)),
$\hat{\mu}_{22}$ is activated ($\hat{\mu}_{22}=-20$). Once $\hat{\mu}_{22}$ is
turned on, type 2 particles move towards the already established type 1
clusters. As the system evolves, clusters will tend to merge if they are close
enough to “sense” the chemical gradients from each other. In Fig. 7(c), the
numbered circles 1–3 highlight the merging process of two clusters into one.
In circle 4, the two small clusters have already merged and we observe a line
of trapped type 2 particles. (Lastly, we note that an earlier limited study
using particle-based simulations of a binary mixture of sphere-dimer motors
has shown that sequential kinetics also leads to distinctive cluster formation
in this system. [62])

Figure 6: Particle orientations in typical clusters. The particle orientation
vector is denoted by the white arrow. The left and right panels of the top row
correspond to Fig. 4(a) and (b), respectively. For the top row,
$\gamma_{\alpha\beta}=0$. The left and right panels of the bottom row
correspond to Fig. 4(d) and (e), respectively. For the bottom row,
$\mu_{\alpha\beta}=0$.

While the heterogeneous clusters formed due to type I (Fig. 4(b)) and type II
(Fig. 4(e)) mechanisms are similar, they can be distinguished by the
orientational order within clusters. In Fig. 4(b), particles aggregate due to
attraction, i.e., they move towards each other with induced phoretic
translational velocities. When they form a cluster, the particles are held
together by attraction. The particle orientation vector is not affected by the
attraction and remains uniformly randomly distributed in clusters (see the top
row of Fig. 6). For type II mechanism shown in Fig. 4(e), particles reorient
due to phoretic alignment. In order to form a cluster, particles have to swim
towards each other. In order to maintain a cluster, particles on the outer
shells have to continue pushing towards the interior of the cluster (see the
bottom row of Fig. 6). A particle at the outer shell with an orientation
vector pointing away from the cluster simply swims away.

Figure 7: Snapshots of the binary mixture undergoing sequential surface
reactions with total area fraction $\phi=0.1$ at different times: (a)
$t/\tau_{D}=0$, (b) $t/\tau_{D}=250$, and (c) $t/\tau_{D}=500$, where
$\tau_{D}=a^{2}/D$ is the diffusive time scale. In terms of the bulk rate, we
have for (b) $k_{b,1}t=15.625$, and for (c) $k_{b,1}t=31.25$. Initially,
$\hat{\mu}_{22}=0$; after (b), $\hat{\mu}_{22}$ is activated,
$\hat{\mu}_{22}=-20$. All other parameters are the same as Fig. 4(b). Numbered
circles 1–3 show the merging process of two clusters into one while circle no.
4 shows a merged cluster.  Figure 8: Maximum cluster size, defined as the
number of particles in a cluster, in the system as a function of time. The
black line, $\mu_{22}<0$, corresponds to the simulation presented in Fig.
4(a). The blue line, $\mu_{22}>0$, corresponds to the simulation presented in
Fig. 4(b). In counting the number of particles in a cluster, both types of
particles are included.

In Fig. 8, we plot the maximum cluster size in the system as a function of
time for simulations corresponding to Figs. 4(a) and 4(b). The maximum cluster
size is defined as the total number of particles, regardless of their type, in
a single cluster. The system first establishes small clusters and then these
clusters may merge due to attraction (see also Fig. 7 and ESI). Such merging
events are signified by the jumps in the cluster size curves shown in Fig. 8.
The maximum cluster size for $\mu_{22}<0$ (black) is larger than that for
$\mu_{22}>0$ (blue) because type 2 particles are attracted to type 1 particles
or clusters when $\mu_{22}<0$. The fluctuations in the curves indicate that
particles can randomly attach to or detach from a cluster as a result of
Brownian motion. The jumps in the cluster size curve are more significant if
two large clusters merge into a single one; this can occur at higher particle
densities or for weaker chemical screening (see Figs. 1 and 6 in the ESI).

Figure 9: The self and cross radial distribution functions (RDFs) of the
binary mixture that forms heterogeneous clusters. The RDFs are computed using
a single frame of the system snapshot at $t/\tau_{D}=2497.5$. The system
parameters are the same as those in Fig. 4(b). The subscripts
$\alpha,\beta=1,2$ denote the particle type. The RDFs are computed with 150
bins for $r/(2a)\in[0,6]$, where $r$ is the radial distance.

For the system parameters used in Figs. 4(b) and 7, the binary mixture forms
heterogeneous clusters with a layering structure. In Fig. 9 we plot the self
and cross radial distribution functions (RDFs) for a snapshot taken at
$t/\tau_{D}=2497.5$. As can be seen in Fig. 4(b), the clusters tend to
approach a hexagonal packing with more loosely attached particles in the
exterior. Because of the layering structure, the 1-1 RDF ($g_{11}$) has the
largest peaks, which is followed by the 1-2 RDF ($g_{12}$). The 2-2 RDF has
the lowest peaks because there are no attraction among type 2 particles; the
peaks in $g_{22}$ largely results from the trapped particles (the layers)
between type 1 particles.

We note that various system parameters may be used to tune the dynamics and
self-organization of the binary mixture. In Fig. 4(b) and 7, the phoretic
coefficients $\mu_{11}$ and $\mu_{22}$ have the same magnitude but opposite
sign, $\hat{\mu}_{22}=-\hat{\mu}_{11}$. By varying the relative magnitude
$\mu_{22}/|\mu_{11}|$ while keeping the signs fixed, the system can be placed
in different states (see ESI, Fig. 3). If $\mu_{22}$ is small, type 2
particles experience a weak attraction to type 1 clusters; as a result, the
clusters are mostly of type 1 particles with many type 2 particles dispersed
in bulk. For intermediate values of $\mu_{22}$, type 2 particles are attached
to the clusters. If $\mu_{22}$ is large, type 1 clusters can be destroyed by
type 2 particles as a result of non-reciprocal attraction. Whenever a blue
cluster forms, it attracts type 2 particles, which then move towards the
cluster and can push through it; as a result, the clusters are destroyed.
Another parameter that may be accessible is the screening length, which can be
modified by tuning the bulk reaction rate (see ESI, Fig. 1). If the screening
length is much smaller than the particle size, no clusters can form since
other particles do not sense the chemical gradient. For large screening length
(or weak screening), the system can form larger clusters even at low particle
densities.

Figure 10: Clustering phase diagram for type 1 particles in the binary mixture
undergoing the sequential reaction. Filled circles denote the clustering state
while open circles imply that type 1 particles do not form clusters. The
system parameters are $\Lambda_{1}=\Lambda_{2}=0.25$,
$Da_{1}=Da_{2}=1/(2\pi)\approx 0.16$,
$\hat{\mu}_{12}=\hat{\mu}_{21}=\hat{\mu}_{22}=0$,
$\hat{\gamma}_{12}=\hat{\gamma}_{21}=\hat{\gamma}_{22}=0$,
$\hat{\alpha}_{1}=\hat{\alpha}_{2}=2$, and the total area fraction is
$\phi=0.1$.

Instead of varying the phoretic coefficients for different colloidal types,
one can also consider different combinations of translational and rotational
phoretic coefficients for the same particle type. We present a phase diagram
in Fig. 10 for the clustering of type 1 particles by varying $\mu_{11}$ and
$\gamma_{11}$ while keeping all other parameters fixed. In the lower left
corner of the diagram, both $\mu_{11}$ and $\gamma_{11}$ are negative, which
gives rise to attractive interactions from both type I and II mechanisms. As a
result, in this region type 1 particles form clusters (filled circles) . When
both $\mu_{11}$ and $\gamma_{11}$ are positive (upper right corner), no
cluster is formed (open circles). Near the phase boundary, for example,
$\hat{\gamma}_{11}=-20$ and $\hat{\mu}_{11}=4$, the translational and
rotational phoretic effects are in competition. For this particular case, the
system can establish clusters; if $\hat{\mu}_{11}$ is increased, the clusters
vanish and the system moves into a dispersed state.

In this section, we have illustrated the dynamics of a binary mixture with
sequential reactions using three chemical species $S_{1}$, $S_{2}$ and
$S_{3}$. The three chemical species form a closed chemical cycle (or loop) in
the sense that all species are produced and consumed within the system, albeit
at different “sites”. As sketched in Fig. 3, $S_{1}$ is produced in the bulk
and consumed on type 1 particles; $S_{2}$ is produced and consumed in the bulk
and on particles; $S_{3}$ is produced on type 2 particles and consumed in the
bulk. (The cycle can be “short-circuited” since $S_{2}$ produced in the bulk
can react with type 2 particles.) In Refs. [63, 64], colloidal mixtures that
participate in reactions that form a closed cycle are considered where all
chemical reactions occur on particle surfaces and chemical species freely
diffuse in bulk. In contrast, the chemical cycle in our system is formed only
when the particles interact with the surrounding chemical environment (The
chemical reactions on the particles do not form a cycle.). Notably, the bulk
reactions in the environment introduces chemical screening, which
fundamentally changes the chemical interaction from long-ranged (pure
diffusion) to short-ranged (screening). By tuning the chemical environment,
i.e., the bulk reaction rate, one can control the system dynamics via the
screening lengths. Because the system in general has multiple screening
lengths, they can be designed to be different, thus may give rise to more
complex dynamics.

##  V Binary mixture with a nonlinear autocatalytic reaction

We now consider a binary mixture with nonlinear reaction kinetics based on the
Selkov model [65]. This model has its origin as simple description of a step
in glycolsis involving the phosphofructokinase enzyme that converts
fructose-6-phosphate and adenosine triphosphate (ATP) to
fructose-1,6-bisphosphate and adenosine diphosphate (ADP). The reversible
version of this model [66] comprises the following reactions: on enzyme-coated
type 1 particles we have
$S_{1}+2S_{2}\xrightleftharpoons[k_{2}^{-}]{k_{2}^{+}}3S_{2}$. In the bulk
fluid, we have $F\xrightleftharpoons[k_{1}^{-}]{k_{1}^{+}}S_{1}$ and
$S_{2}\xrightleftharpoons[k_{3}^{-}]{k_{3}^{+}}G$, where $F$ and $G$ are held
fixed by reservoirs and their concentrations do not enter the reaction
kinetics directly. On type 2 particles, we have $S_{2}\xrightarrow{k_{4}}G$
(see Fig. 11 for a schematic). The destruction of $S_{2}$ on type 2 is not
part of the Selkov model but is introduced to consume excess $S_{2}$ produced
by the autocatalytic reaction on type 1. In the Selkov model $S_{1}$ and
$S_{2}$ signify ATP and ADP, respectively, but here we regard it as a
representative model for nonlinear autocatalytic kinetics.

The concentration of $S_{1}$, denoted by $c_{1}$, is governed by

$\displaystyle\frac{\partial c_{1}}{\partial t}=$ $\displaystyle
D\nabla^{2}c_{1}+k_{1}^{+}-k_{1}^{-}c_{1}$
$\displaystyle-\sum_{i\in\sigma_{1}}\left(k_{2}^{+}c_{1}c_{2}^{2}-k_{2}^{-}c_{2}^{3}\right)\delta(\bm{x}-\bm{X}_{i}).$
(11)

The concentration of $S_{2}$, or $c_{2}$, satisfies the equation

$\displaystyle\frac{\partial c_{2}}{\partial t}$ $\displaystyle=$
$\displaystyle D\nabla^{2}c_{2}+k_{3}^{-}-k_{3}^{+}c_{2}$ (12)
$\displaystyle+\sum_{i\in\sigma_{1}}\left(k_{2}^{+}c_{1}c_{2}^{2}-k_{2}^{-}c_{2}^{3}\right)\delta(\bm{x}-\bm{X}_{i})$
$\displaystyle-\sum_{j\in\sigma_{2}}k_{4}c_{2}\delta(\bm{x}-\bm{X}_{j}).$

Because of the linear reactions in the bulk, the chemical system is screened;
the nonlinear reactions occur only on the surface of the particles. The
inverse screening lengths can be defined as $\lambda_{1}=\sqrt{k_{1}^{-}/D}$
and $\lambda_{2}=\sqrt{k_{3}^{+}/D}$. Scaling the screening lengths by the
radius $a$, we have the non-dimensional screening lengths given by
$\Lambda_{1}=\lambda_{1}a$ and $\Lambda_{2}=\lambda_{2}a$. For simplicity, we
only consider the phoretic motion induced by concentration gradients and set
the self-diffusiophoretic coefficients and the rotational diffusiophoretic
coefficients to be zero. In the simulations, the number of type 1 particles is
the same as type 2, $N_{1}=N_{2}=500$; the area fraction counting both types
of particles is $\phi=0.1$.

In the single and sequential reactions considered in previous sections, the
total concentration of all species is conserved, which provides a reference
concentration for non-dimensionalization. In the reversible Selkov model, such
a simple conservation statement is absent. To non-dimensionalize the system,
we define the reference concentration $\tilde{c}=k_{1}^{+}a^{2}/D$, which is
obtained by balancing the diffusive term with the bulk production of $S_{1}$
and taking the radius $a$ of the particles as the length scale. With this, we
introduce the following non-dimensional reaction parameters:

$\hat{k}_{2}^{\pm}=\frac{k_{2}^{\pm}\tilde{c}^{2}}{D}=\frac{k_{2}^{\pm}\left(k_{1}^{+}\right)^{2}a^{4}}{D^{3}},\quad\hat{k}_{3}^{-}=\frac{k_{3}^{-}}{k_{1}^{+}},\quad\text{and}\quad\hat{k}_{4}=\frac{k_{4}}{D}.$
(13)

The non-dimensional phoretic coefficients are given by
$\hat{\mu}_{ij}=\mu_{ij}\tilde{c}/D$ for $i,j=1,2$.

Figure 11: Schematic of the Selkov reaction scheme. The autocatalytic reaction
occurs on the first type of particle while type 2 particles consume $S_{2}$
with a first-order reaction. Reversible reactions that produce or consume
$S_{1}$ and $S_{2}$ are present in the bulk and the species $F$ and $G$ are
maintained by reservoirs and their concentrations do not enter the reaction
kinetics.  Figure 12: Heterogeneous clusters of the binary Selkov model. (a):
The spatially-averaged chemical concentrations vs. the scaled time
$t/\tau_{D}$. (b): The system snapshot at $t/\tau_{D}=200$ and the background
contour plot denotes the concentration $c_{1}$. (c): The system snapshot at
$t/\tau_{D}=200$ (same as (b)) and the background contour plot denotes the
concentration $c_{2}$. In this simulation, the self-diffusiophoretic
coefficients and the angular diffusiophoretic coefficients are zero. The
active velocities of the particles result from the translational
diffusiophoretic motion: $\hat{\mu}_{11}=\hat{\mu}_{12}=\hat{\mu}_{21}=-2$ and
$\hat{\mu}_{22}=2$. Other physical parameters used are as follows:
$\Lambda_{1}=0.1,\Lambda_{2}=1,\hat{k}_{2}^{\pm}=0.04,\hat{k}_{3}^{-}=0.45$
and $\hat{k}_{4}=4$. Notice that the dimensional screening length of $c_{1}$
is $10a$, or five particle diameters, while the screening length of $c_{2}$ is
equal to the radius.

In Fig. 12, we show the clustering behavior of the reversible Selkov model
when $\hat{\mu}_{11}=\hat{\mu}_{12}=\hat{\mu}_{21}=-10$ and
$\hat{\mu}_{22}=10$. With these phoretic coefficients, the binary mixture
forms localized clusters that contain both types of particles. Since
$\hat{\mu}_{11}<0$ and $\hat{\mu}_{21}<0$, particles are attracted to the the
sinks of $S_{1}$, and at long times the clusters are localized in regions of
low concentration of $S_{1}$. The concentration of $S_{1}$, or $c_{1}$, is
shown in the contour plot in Fig. 12(b). The same clusters are shown in Fig.
12(c) with the contour representing the concentration of $S_{2}$ ($c_{2}$). We
notice that the concentration gradient of $S_{2}$ is not visible in Fig. 12(c)
because the dimensional screening length of $c_{2}$ is equal to the radius of
the particle. As a result, gradients of $c_{2}$ are highly localized at each
particle and large scale concentration gradients of $c_{2}$ can not be
established. Within a cluster, type 1 particles produce $S_{2}$ and
subsequently $S_{2}$ is consumed by type 2 particles. Because the screening
length of $c_{1}$ is five particle diameters, the gradients of $c_{1}$ are
visible. In the results shown in Fig. 12, the dominant attraction mechanism is
from $\bm{\nabla}c_{1}$, which dictates the clustering dynamics.

In the simulation of Fig. 12, the initial concentrations for both $S_{1}$ and
$S_{2}$ are zero. The spatially averaged concentrations are plotted in Fig.
12(a) as a function of the scaled time. Initially, the bulk production of
$S_{1}$ and $S_{2}$ dominates and nonzero concentrations of $S_{1}$ and
$S_{2}$ are gradually established. Once type 1 particles ‘detect’ a finite
concentration of $S_{1}$ and $S_{2}$, the autocatalytic reaction will be
activated and we observe a sharp decline in $\overline{c}_{1}$ and an increase
in $\overline{c}_{2}$. During this time, the system is also starting to form
small clusters. At long times, we observe mostly steady average concentrations
for both $S_{1}$ and $S_{2}$ with fluctuations. The small clusters are
attracted to each other and at long times form larger clusters as shown in
Figs. 12(b) and 12(c).

Figure 13: Schematic of two colloidal molecules composed of three particles of
type 1 and four particles of type 2. The molecule in (a) has rotational
symmetry; the molecule in (b) has a line of symmetry as denoted by the dashed
line.

The binary mixture with Selkov kinetics can also form self-assembled colloidal
molecules. In the absence of Brownian motion, colloidal molecules that are
stationary or exhibit net motion can be obtained. Using a simple model in
which the particles act as sources or sinks for a single otherwise diffusing
chemical field, previous studies have shown that stable colloidal molecules of
various compositions can be formed [67, 68]. In contrast, in our current model
we have two chemical species and nonlinear reaction kinetics. For simplicity,
we consider in isolation a collection of seven particles in which three are of
type 1 (blue in Fig. 11) and the other four are of type 2 (gold in Fig. 11).
In Fig. 13, we show examples of two colloidal molecules that can be formed by
taking the same parameters as those in Fig. 12. The molecule in Fig. 13(a) has
rotational symmetry; in the absence of Brownian fluctuations, it can maintain
a stationary state. The molecule in Fig. 13(b) has a line of symmetry and
exhibits net translational motion along this line. The blue particles consume
$S_{1}$ and create a concentration gradient with lower concentrations of
$S_{1}$ near the blue particles and higher concentrations near the gold ones
(see Fig. 14). Because the gold particles are attracted to regions of low
concentrations of $S_{1}$, the molecule translates in a direction with the
blue-particle side at its head (see ESI, Movie 1).

Figure 14: Snapshot of the translating molecule and the concentration field
$c_{1}$. The direction of motion is denoted by the arrow.

By treating the colloidal particles as a spatially homogeneous system, we note
that the chemical system is at a stable fixed point. Denoting the spatially
constant densities of type 1 and 2 particles as $n_{1}$ and $n_{2}$,
respectively, one may write
$\partial_{t}c_{1}=D\nabla^{2}c_{1}+k_{1}^{+}-k_{1}^{-}c_{1}-(k_{2}^{+}c_{1}c_{2}^{2}-k_{2}^{-}c_{2}^{3})n_{1}$,
and
$\partial_{t}c_{2}=D\nabla^{2}c_{2}+k_{3}^{-}-k_{3}^{+}c_{2}+(k_{2}^{+}c_{1}c_{2}^{2}-k_{2}^{-}c_{2}^{3})n_{1}-k_{4}c_{2}n_{2}$.
At steady state for system parameters of Fig. 12, these equations give
$c_{1}\approx 19.764$ and $c_{2}\approx 0.087$. A linear stability analysis
about this fixed point reveals that the chemical system is stable to
inhomogeneous concentration perturbations for all wavenumbers, which suggests
that the instability (i.e., emergence of a non-homogeneous density state)
shown in Fig. 12 results from the clustering instabilities induced by phoretic
interactions. In a full description that includes the equation of motion for
the inhomogeneous density fields, for appropriate system parameters, the
chemical system itself may be unstable and the interaction of the chemical
instability with the clustering instability may give rise to more complex
bifurcations and dynamics.

The above results provide some simple examples of the phenomena that can be
seen in systems where the chemical mechanism involves nonlinear kinetics on
specific colloidal particles, along with chemical reactions on other colloids
and in the surrounding fluid. The possible types of behavior is far richer
than that explored here. For example, the formation of clusters with specific
compositions and surface reaction kinetics constitute a compartmentalized
reaction-diffusion system where portions of the reaction network reside in
spatially localized compartments that communicate by diffusion. Depending on
the diffusion length relative to the distance between compartments, the system
dynamics can change its bifurcation structure. [69] A laboratory example is
the Belousov-Zhabotinsky reaction carried out in miroemulsions. [70] In this
system different parts of the mechanism are confined to the different phases
in the emulsion, and the dynamics depends strongly on the composition of the
microemulsion. In biological systems reactions are often compartmentalized in
specialized organelles and communicate with reactions in the surrounding
cytoplasm through diffusion [71]. Dynamic clustering in active colloidal
systems, where the reaction network distributed over different system
constituents, is expected to display a similarly rich bifurcation structure
that now depends on the dynamics of the clustering process. While such
phenomena have not been explored in this paper, they deserve further study.

##  VI Conclusion

In this paper, we have considered the collective motion and self-organization
of chemically active colloids propelled by complex chemical reactions
involving multiple chemical species and potentially nonlinear kinetics. The
multiple chemical species are coupled by a network of chemical reactions. In
general, the results have shown that active colloids can form clusters in the
presence of phoretic attractions. The characteristics of the clusters,
however, depend on the reaction kinetics and the resulting chemical fields and
phoretic coefficients. In binary mixtures, we have shown that diverse
clustering behaviors can be obtained including the clustering of one type of
colloid while the other type is dispersed, heterogeneous clusters containing
both types of particles, and heterogeneous clusters with one type of colloid
on the periphery. We have also shown that chemical screening typically leads
to localized clusters in the dilute regime. With a single chemical species,
the tunable parameters for self-organization are limited. By allowing multiple
chemical reactions both on the colloids and in the surrounding fluid, our
study shows that the design space of colloidal self-organization can be
broadened by including chemical reaction kinetics as an additional parameter.

The complex chemical reaction kinetics considered in our work can be
implemented using enzyme-coated colloidal particles or natural cells. In the
context of designing micro- and nano-motors for biomedical applications,
enzymatic reactions that make use of readily available biofuel is desirable.
For example, an endogenous enzyme-powered Janus platelet micro-motor is shown
to exhibit phoretic motion due to the uneven decomposition of urea in
biofluids[49]. Micro- and nano-motors powered by other enzymatic reactions
have also been realized[50, 51, 42, 43, 44, 45, 46, 47, 48, 49]. The enzymatic
reactions that can be used for self-propulsion and self-organization may
depend on the biological environment and whether the required biofuel species
are present in such an environment. As a result, it is beneficial to develop a
variety of enzyme-powered motors that make use of different reaction schemes.

While we focused on demonstrating the diverse clustering behaviors that can be
achieved, our more general model outlined in section II includes nonlinear
kinetics and bulk reactions. In both the sequential and nonlinear models
considered in this paper, the bulk reaction kinetics are linear. We note that
generic nonlinear reaction kinetics in the bulk can be employed to establish
external gradients that influence the clustering dynamics of the particles. In
the case of linear bulk kinetics, the resulting multiple screening lengths can
also be used to control the system dynamics. In addition to phenomena related
to compartmentalized reaction-diffusion systems discussed above, the research
reported here can also form the basis for of chemical computation in active
colloidal systems. The surface reactions (especially those that are nonlinear)
can be constructed to be logic gates [72], and clusters of colloids supporting
different gates can lead to the self assembly of logic circuits. In this way
it may be possible for the system to carry out chemical computations to
influence its dynamics to perform specific tasks.

To study the self-organization of active colloids mediated by complex chemical
reactions, we presented a minimal computational model. In constructing our
model, for simplicity we used a monopole approximation for the chemical fields
following previous studies and we have neglected hydrodynamic interactions. To
account for the near-field chemical interactions, one needs to solve the
reaction-diffusion equations subject to surface reaction boundary conditions.
For constant or linear surface reaction kinetics in the absence of bulk
reactions, the reaction-diffusion equation of a single chemical species in the
presence of one or two active particles have been considered[73, 74, 75, 76,
77, 78]. Mesoscopic particle-based [79] and other numerical methods [80, 81]
can be used to take into account full complex chemical dynamics as well as
hydrodynamic interactions.

## Author Contributions

Z. P. and R. K. contributed at all stages of this work. Z. P. developed the
simulation software and performed simulations.

## Conflicts of interest

There are no conflicts to declare.

## Acknowledgements

This work was supported in part by the Natural Sciences and Engineering
Research Council (NSERC) of Canada. Computations were performed on SciNet HPC
Consortium computers. SciNet is funded by the Canada Foundation for
Innovation, the Government of Ontario, the Ontario Research Excellence Fund,
and the University of Toronto.

##  Appendix A Numerical and Brownian dynamics simulations

Consider a collection of $N$ spherical particles in a square simulation box
$[-L,L]\times[-L,L]$ with periodic boundary conditions in both directions. To
evolve the reaction-diffusion equations numerically, we discretize space using
a uniform grid of spacing $h=2L/M$, where $M$ is the number of intervals in
each direction. The reaction-diffusion equations are discretized using second-
order centered finite differences. Taking the second equation in the
sequential reaction scheme as an example, the discretization gives

$\displaystyle\frac{\text{d}c_{2}[i,j]}{\text{d}t}=D\frac{c_{2}[i+1,j]-2c_{2}[i,j]+c_{2}[i-1,j]}{h^{2}}$
$\displaystyle+D\frac{c_{2}[i,j+1]-2c_{2}[i,j]+c_{2}[i,j-1]}{h^{2}}$
$\displaystyle-k_{b,1}c_{2}[i,j]+k_{b,2}c_{3}[i,j]$
$\displaystyle+\kappa_{1}\sum_{k=1}^{N_{1}}\frac{c_{1}[i_{k},j_{k}]}{h^{2}}-\kappa_{2}\sum_{k=N_{1}+1}^{N}\frac{c_{2}[i_{k},j_{k}]}{h^{2}},$
(14)

where $c_{1}$ is now the values at grid points and $c_{1}[i,j]$ denotes the
node value at indices ($i,j$), and $i_{k}$ and $j_{k}$ denotes the index of
the grid cell in which a particle resides. Because we need to know the
particle locations in order to evolve the reaction-diffusion equation in time,
an occupancy matrix (bit array) is maintained for each particle type. We then
discretize time using the first-order explicit Euler scheme. At each time
step, the reaction-diffusion equations are first evolved, from which we
calculate the active linear and angular velocities. With these velocities, we
then perform Brownian dynamics simulations that will update the particle
positions and orientations.

We simulate the Langevin equations of motion using Brownian dynamics, where
the fluctuating velocities are treated using the Euler-Maruyama scheme. The
hard-sphere interactions among particles are treated using a potential-free
geometric optimization approach and all possible collisions during a discrete
time step are resolved simultaneously. This scheme requires solving a linear
complementarity problem, where a projected gradient descent method is used
(see Refs.[82, 83] for details). The simulations are performed using an in-
house program that runs on CUDA-enabled NVIDIA GPUs.

## References

  * Ndlec _et al._ [1997] F. Ndlec, T. Surrey, A. C. Maggs and S. Leibler, _Nature_ , 1997, 389, 305–308.
  * Surrey _et al._ [2001] T. Surrey, F. Nédélec, S. Leibler and E. Karsenti, _Science_ , 2001, 292, 1167–1171.
  * Lee and Kardar [2001] H. Y. Lee and M. Kardar, _Phys. Rev. E_ , 2001, 64, 056113.
  * Ross _et al._ [2019] T. D. Ross, H. J. Lee, Z. Qu, R. A. Banks, R. Phillips and M. Thomson, _Nature_ , 2019, 572, 224–229.
  * Kearns [2010] D. B. Kearns, _Nat. Rev. Microbiol_ , 2010, 8, 634–644.
  * Ariel _et al._ [2015] G. Ariel, A. Rabani, S. Benisty, J. D. Partridge, R. M. Harshey and A. Be’Er, _Nat. Commun._ , 2015, 6, 8396.
  * Be’er and Ariel [2019] A. Be’er and G. Ariel, _Mov. Ecol._ , 2019, 7, 1–17.
  * Fily and Marchetti [2012] Y. Fily and M. C. Marchetti, _Phys. Rev. Lett._ , 2012, 108, 235702\.
  * Redner _et al._ [2013] G. S. Redner, M. F. Hagan and A. Baskaran, _Phys. Rev. Lett._ , 2013, 110, 055701.
  * Cates and Tailleur [2015] M. E. Cates and J. Tailleur, _Annu. Rev. Condens. Matter Phys._ , 2015, 6, 219–244.
  * Takatori and Brady [2015] S. C. Takatori and J. F. Brady, _Phys. Rev. E_ , 2015, 91, 032117\.
  * Solon _et al._ [2018] A. P. Solon, J. Stenhammar, M. E. Cates, Y. Kafri and J. Tailleur, _New J. Phys._ , 2018, 20, 075001.
  * Hermann _et al._ [2019] S. Hermann, P. Krinninger, D. de las Heras and M. Schmidt, _Phys. Rev. E_ , 2019, 100, 052604.
  * Omar _et al._ [2023] A. K. Omar, H. Row, S. A. Mallory and J. F. Brady, _Proc. Natl. Acad. Sci. U.S.A._ , 2023, 120, e2219900120.
  * Zhao _et al._ [2023] H. Zhao, A. Košmrlj and S. S. Datta, _Phys. Rev. Lett._ , 2023, 131, 118301.
  * Toner and Tu [1998] J. Toner and Y. Tu, _Phys. Rev. E_ , 1998, 58, 4828.
  * Cavagna and Giardina [2014] A. Cavagna and I. Giardina, _Annu. Rev. Condens. Matter Phys._ , 2014, 5, 183–207.
  * Bialek _et al._ [2012] W. Bialek, A. Cavagna, I. Giardina, T. Mora, E. Silvestri, M. Viale and A. M. Walczak, _Proc. Natl. Acad. Sci. U.S.A._ , 2012, 109, 4786–4791.
  * Weihs [1973] D. Weihs, _Nature_ , 1973, 241, 290–291.
  * Berg [1975] H. C. Berg, _Annu. Rev. Biophys. Bioeng._ , 1975, 4, 119–136.
  * Adler [1975] J. Adler, _Annu. Rev. Biochem._ , 1975, 44, 341–356.
  * Keller and Segel [1971] E. F. Keller and L. A. Segel, _J. Theor. Biol._ , 1971, 30, 225–234.
  * Mallory _et al._ [2018] S. A. Mallory, C. Valeriani and A. Cacciuto, _Annu. Rev. Phys. Chem._ , 2018, 69, 59–79.
  * Theurkauff _et al._ [2012] I. Theurkauff, C. Cottin-Bizonne, J. Palacci, C. Ybert and L. Bocquet, _Phys. Rev. Lett._ , 2012, 108, 268303.
  * Thakur and Kapral [2012] S. Thakur and R. Kapral, _Phys. Rev. E_ , 2012, 85, 026121.
  * Palacci _et al._ [2013] J. Palacci, S. Sacanna, A. P. Steinberg, D. J. Pine and P. M. Chaikin, _Science_ , 2013, 339, 936–940.
  * Pohl and Stark [2014] O. Pohl and H. Stark, _Phys. Rev. Lett._ , 2014, 112, 238303.
  * Pohl and Stark [2015] O. Pohl and H. Stark, _Eur. Phys. J. E_ , 2015, 38, 1–11.
  * Huang _et al._ [2019] M. J. Huang, J. Schofield, P. Gaspard and R. Kapral, _J. Chem. Phys._ , 2019, 150, 124110.
  * Varma _et al._ [2018] A. Varma, T. D. Montenegro-Johnson and S. Michelin, _Soft Matter_ , 2018, 14, 7155–7173.
  * Schmidt _et al._ [2019] F. Schmidt, B. Liebchen, H. Löwen and G. Volpe, _J. Chem. Phys._ , 2019, 150, 094905.
  * Meredith _et al._ [2020] C. H. Meredith, P. G. Moerman, J. Groenewold, Y.-J. Chiu, W. K. Kegel, A. van Blaaderen and L. D. Zarzar, _Nat. Chem._ , 2020, 12, 1136–1142.
  * Saha _et al._ [2014] S. Saha, R. Golestanian and S. Ramaswamy, _Phys. Rev. E_ , 2014, 89, 062316.
  * Keller and Segel [1970] E. F. Keller and L. A. Segel, _J. Theor. Biol._ , 1970, 26, 399–415.
  * Brenner _et al._ [1998] M. P. Brenner, L. S. Levitov and E. O. Budrene, _Biophys. J._ , 1998, 74, 1677–1693.
  * Mittal _et al._ [2003] N. Mittal, E. O. Budrene, M. P. Brenner and A. Van Oudenaarden, _Proc. Natl. Acad. Sci. U.S.A._ , 2003, 100, 13259–13263.
  * Peruani _et al._ [2012] F. Peruani, J. Starruß, V. Jakovljevic, L. Søgaard-Andersen, A. Deutsch and M. Bär, _Phys. Rev. Lett._ , 2012, 108, 098102.
  * Liebchen _et al._ [2015] B. Liebchen, D. Marenduzzo, I. Pagonabarraga and M. E. Cates, _Phys. Rev. Lett._ , 2015, 115, 258301.
  * Liebchen _et al._ [2017] B. Liebchen, D. Marenduzzo and M. E. Cates, _Phys. Rev. Lett._ , 2017, 118, 268001.
  * Gaspard and Kapral [2020] P. Gaspard and R. Kapral, _Research_ , 2020, 2020, 9739231.
  * Agudo-Canalejo and Golestanian [2019] J. Agudo-Canalejo and R. Golestanian, _Phys. Rev. Lett._ , 2019, 123, 018101.
  * Dey _et al._ [2015] K. K. Dey, X. Zhao, B. M. Tansi, W. J. Méndez-Ortiz, U. M. Córdova-Figueroa, R. Golestanian and A. Sen, _Nano Lett._ , 2015, 15, 8311–8315.
  * Schattling _et al._ [2015] P. Schattling, B. Thingholm and B. Städler, _Chem. Mater._ , 2015, 27, 7412–7418.
  * Ma _et al._ [2016] X. Ma, A. C. Hortelão, T. Patiño and S. Sánchez, _ACS Nano_ , 2016, 10, 9111–9122.
  * Ma _et al._ [2016] X. Ma, X. Wang, K. Hahn and S. Sánchez, _ACS Nano_ , 2016, 10, 3597–3605.
  * Zhao _et al._ [2018] X. Zhao, K. Gentile, F. Mohajerani and A. Sen, _Acc. Chem. Res._ , 2018, 51, 2373–2381.
  * Toebes _et al._ [2019] B. J. Toebes, F. Cao and D. A. Wilson, _Nat. Commun._ , 2019, 10, 5308\.
  * Arqué _et al._ [2019] X. Arqué, A. Romero-Rivera, F. Feixas, T. Patiño, S. Osuna and S. Sánchez, _Nat. Commun._ , 2019, 10, 2826.
  * Tang _et al._ [2020] S. Tang, F. Zhang, H. Gong, F. Wei, J. Zhuang, E. Karshalev, B. Esteban-Fernández de Ávila, C. Huang, Z. Zhou, Z. Li _et al._ , _Sci. Robot._ , 2020, 5, eaba6137.
  * Gao _et al._ [2013] W. Gao, M. D’Agostino, V. Garcia-Gradilla, J. Orozco and J. Wang, _Small_ , 2013, 9, 467–471.
  * Schattling _et al._ [2017] P. S. Schattling, M. A. Ramos-Docampo, V. Salgueiriño and B. Städler, _ACS Nano_ , 2017, 11, 3973–3983.
  * Sun _et al._ [2019] J. Sun, M. Mathesh, W. Li and D. A. Wilson, _ACS Nano_ , 2019, 13, 10191–10200.
  * Chen _et al._ [2019] Z. Chen, T. Xia, Z. Zhang, S. Xie, T. Wang and X. Li, _Chem. Eng. J._ , 2019, 375, 122109.
  * Yuan _et al._ [2021] H. Yuan, X. Liu, L. Wang and X. Ma, _Bioact. Mater._ , 2021, 6, 1727–1749.
  * Mathesh _et al._ [2020] M. Mathesh, J. Sun and D. A. Wilson, _J. Mater. Chem. B_ , 2020, 8, 7319–7334.
  * Kondrat _et al._ [2016] S. Kondrat, O. Zimmermann, W. Wiechert and E. von Lieres, _Eur. Phys. J. E_ , 2016, 39, 1–10.
  * Gaspard and Kapral [2017] P. Gaspard and R. Kapral, _J. Chem. Phys._ , 2017, 147, 211101\.
  * Gaspard and Kapral [2019] P. Gaspard and R. Kapral, _Adv. Phys. X_ , 2019, 4, 1602480.
  * Corato1 and Pagonabarraga [2022] M. D. Corato1 and I. Pagonabarraga, _J. Chem. Phys._ , 2022, 157, 084901\.
  * Huang _et al._ [2017] M.-J. Huang, J. Schofield and R. Kapral, _New J. Phys._ , 2017, 19, 125003.
  * Huang _et al._ [2018] M.-J. Huang, J. Schofield, P. Gaspard and R. Kapral, _J. Chem. Phys._ , 2018, 149, 024904.
  * Colberg _et al._ [2014] P. H. Colberg, S. Y. Reigh, B. Robertson and R. Kapral, _Acc. Chem. Res._ , 2014, 47, 3504–3511.
  * Ouazan-Reboul _et al._ [2023] V. Ouazan-Reboul, R. Golestanian and J. Agudo-Canalejo, _Phys. Rev. Lett._ , 2023, 131, 128301.
  * Ouazan-Reboul _et al._ [2023] V. Ouazan-Reboul, J. Agudo-Canalejo and R. Golestanian, _Nat. Commun._ , 2023, 14, 4496.
  * Sel’Kov [1968] E. Sel’Kov, _Eur. J. Biochem._ , 1968, 4, 79–86.
  * Richter _et al._ [1981] P. H. Richter, P. Rehmus and J. Ross, _Prog. Theor. Phys._ , 1981, 66, 385–405.
  * Soto and Golestanian [2014] R. Soto and R. Golestanian, _Phys. Rev. Lett._ , 2014, 112, 068301\.
  * Soto and Golestanian [2015] R. Soto and R. Golestanian, _Phys. Rev. E_ , 2015, 91, 052304\.
  * Chavez and Kapral [2000] F. Chavez and R. Kapral, _Phys. Rev. E_ , 2000, 63, 016211.
  * Vanag and Epstein [2001] V. Vanag and I. Epstein, _Science_ , 2001, 284, 835.
  * Glass and Kauffman [1973] L. Glass and S. Kauffman, _J. Theor. Biol._ , 1973, 39, 103.
  * Chen _et al._ [2023] J.-X. Chen, J.-Q. Hu and R. Kapral, _arXiv:2302.07670_ , 2023.
  * Michelin _et al._ [2013] S. Michelin, E. Lauga and D. Bartolo, _Phys. Fluids_ , 2013, 25, 061701\.
  * Michelin and Lauga [2014] S. Michelin and E. Lauga, _J. Fluid Mech._ , 2014, 747, 572–604.
  * Michelin and Lauga [2015] S. Michelin and E. Lauga, _Eur. Phys. J. E_ , 2015, 38, 1–16.
  * Michelin and Lauga [2017] S. Michelin and E. Lauga, _Sci. Rep._ , 2017, 7, 42264.
  * Nasouri and Golestanian [2020] B. Nasouri and R. Golestanian, _Phys. Rev. Lett._ , 2020, 124, 168003\.
  * Nasouri and Golestanian [2020] B. Nasouri and R. Golestanian, _J. Fluid Mech._ , 2020, 905, A13.
  * Kapral [2008] R. Kapral, _Adv. Chem. Phys._ , 2008, 140, 89–146.
  * Varma and Michelin [2019] A. Varma and S. Michelin, _Phys. Rev. Fluids_ , 2019, 4, 124204\.
  * Rojas-Pérez _et al._ [2021] F. Rojas-Pérez, B. Delmotte and S. Michelin, _J. Fluid Mech._ , 2021, 919, A22.
  * Yan _et al._ [2019] W. Yan, H. Zhang and M. J. Shelley, _J. Chem. Phys._ , 2019, 150, 064109\.
  * Yan _et al._ [2020] W. Yan, E. Corona, D. Malhotra, S. Veerapaneni and M. Shelley, _J. Comput. Phys._ , 2020, 416, 109524.