**4. Dynamic clustering and chemotactic collapse in a system of active colloids 49**

**4.3. Model for active colloids with interactions transmitted by a chemical field 51**

We model a system ofN active colloids that effectively interact via diffusiophoretic forces induced by a chemical field. Our model system resembles the experimental setup of Ref.

[14] consisting of a collection of Janus-particles half covered with gold on one side and platinum on the other.

### 4.3.1. Derivation of the model equations

Typically, the colloids are heavier than the suspending solvent. They settle to the bottom of the experimental cell, where they form a colloidal monolayer. To induce activity, Ref.

[14] employs gold particles half covered with platinum, which catalyzes the separation
of H_{2}O_{2} towards water and oxygen. Through a combination of self-diffusio- and

elec-trophoresis the particles self-propel (see i.e., [123–125]) with a swimming speed v_{0} along
direction e, which we assume to be fixed in the particle. During the catalytic reactions
at the surfaces of the colloids fuel is converted into several chemical products. This has
two consequences. First of all, in the surrounding of the particles fuel is diminished. If
the swimming speed is very sensible to the fuel concentration, particles will move slower
in regions with large colloidal density. As we have seen in Sec. 2.2.3 a sufficiently strong
slow-down at regions with higher density leads to phase separation. However, in the
experiments of Ref. [14] a 20-fold increase of the fuel was needed to increase v_{0} by a
factor of about 2.5. Thus small inhomogeneities of the chemical field do not strongly
alter the colloidal velocity. Applying the nomenclature of Sec. 3.2.4, we can thus neglect
chemokinesis. On the other hand, as we shall see, chemotaxis does play a crucial role.

A second consequence of the particles consuming fuel is that spatial gradients of re-actants and products form around them. They induce diffusiophoretic motion of neigh-boring particles. Eq. (3.41) states the translational swimming velocity resulting from diffusiophoresis. We neglect the quadrupolar term for our half-coated Janus particles.

Defining ζ_{tr} :=−⟨ζ⟩ with the slip velocity coefficient given by Eq. (3.34), we obtain:

v_{D} =−ζ_{tr}∇c. (4.1)

Defining the rotational diffusiophoretic parameter ζ_{rot} via (9/4a)⟨ζn⟩= −ζ_{rot}e_{i} with e_{i}
being the swimming direction of the i-th particle, Eq. (3.42) takes the form:

ω_{D} =−ζ_{rot}e_{i}×∇c . (4.2)

The slip velocity coefficient ζ depends on the surface potential, with which the solute
interacts with the colloid surface. So, by choosing appropriate materials for the Janus
colloids and their caps, the phoretic parameters ζ_{rot} and ζ_{tr} might each be tuned to have
a positive or negative value. The phoretic parameters also depend on the geometry and
number of caps covering the colloidal surface [40].

We now formulate the overdamped Langevin equations for position r_{i} and orientation
e_{i} of the i-th colloid. Adding up active and phoretic contributions, the deterministic
translational velocity becomes ˙ri =v0ei+vD, and the time variation of the orientation
vector is described by the kinematic equation ˙e_{i} = ω_{D} ×e_{i}. Then, we use Eqs.. (4.1)
and (4.2) together with (e_{i}×∇c)×e_{i} = (1−e_{i}⊗e_{i})∇c to write down the overdamped
Langevin equations for each colloid:

r˙i = v0ei−ζtr∇c(ri) +ξi, (4.3)
e˙_{i} = −ζ_{rot}(1−e_{i}⊗e_{i})∇c(r_{i}) +µ_{i}×e_{i}. (4.4)
Here, ξ_{i} is translational and µ_{i} rotational Gaussian white noise with zero mean and
respective time correlation functions⟨ξ_{i}(t)⊗ξ_{i}(t^{′})⟩= 2D_{tr}1δ(t−t^{′}) and⟨µ_{i}(t)⊗µ_{i}(t^{′})⟩=
2D_{rot}1δ(t−t^{′}), where we have introduced the translational (D_{tr}) and rotational (D_{rot})

Other than, for example Ref. [40], we consider only the chemical concentration field c which is diminished by the fuel consumption of the colloidal particles. It evolves according to the diffusion equation:

˙

c(r) = D_{c}∇^{2}c−k

N

i=1

δ(r−r_{i}), (4.6)

where the colloids enter as sinks diminishing the field with ratek. Typically, the chemical diffuses much faster along a particle radius than the swimming colloids need for the same distance. So, we can neglect the time derivative ofcin Eq. (4.6) meaning that each colloid instantly establishes a static chemical sink around itself, which moves with the colloid.

The chemical diffuses in three dimensions and we use the static solution of Eq. (4.6),
c_{3d}(r) = c_{0}−(k/4πD_{c})N

i=11/|r−r_{i}|, wherec_{0}is the background chemical concentration.

Strictly speaking, we would need to implement a non-flux boundary condition at the boundary of the infinite half-space, but this will not change the basic 1/r dependence of the chemical field. To introduce a two-dimensional concentration field, in which the colloidal monolayer moves, we integrate over a thin layer of thickness h= 2a and obtain c2d(r) = h

0 c3d(r)dz ≈hc3d. Finally, the chemical concentration field reads:

c_{2d}(r) = hc_{0}− kh
4πD_{c}

N

i=1

1

|r−r_{i}|. (4.7)

So, self-phoretic colloids induce long-range chemical concentration gradients and via Eq.

(4.3) and (4.4) long-range interactions between the colloids, which have indeed been experimentally measured by Ref. [35].

We modeled each particle as a chemical sink. Hence, if the translational
diffusio-phoretic parameter ζ_{tr} is positive, the diffusiophoretic velocity v_{D} from Eq. (4.1) is
di-rected towards the neighboring particles and it acts as an effective attraction, while
ζ_{tr}<0 gives rise to an effective repulsion. Similarly, a positive rotational parameter ζ_{rot}
in Eq. (4.2) rotates the swimming direction of an active colloid towards a neighboring
chemical sink and the colloid moves towards the neighbor. Hence, rotational phoresis

also acts like an attractive colloidal interaction while it becomes repulsive for ζ_{rot} < 0.

Tuning the signs of the phoretic parameters, one can realize different types and combi-nations of effective interactions and, as we shall show in the following chapters, thereby induce a variety of collective dynamics.

### 4.3.2. Rescaling

To reduce the parameters of our system, we rescale time by tr = 1/(2Drot) and length
by l_{r} =

D_{tr}/D_{rot} = 2.33a such that the rescaled diffusion constants are set to one
in the Langevin equations (4.3) and (4.4). Using the thermal values from the
Stokes-Einstein relations for the diffusion coefficients, one has

Dtr/Drot = 1.15a. Since the
colloids move close to the bottom wall, this value is changed. Indeed, in the experiments
of Ref. [14] l_{r} = 1.79a is measured. For historical reasons, we use here l_{r} = 2.33a,
which does not change the qualitative behavior of our system. The actual value of lr

is needed for implementing screening of the phoretic interaction inside a colloidal clus-ter, as explained above. After rescaling we are left with four relevant system parame-ters: the P´eclet number Pe =v0/(2√

DtrDrot), the rescaled translational diffusiophoretic
parameter ζ_{tr}kh/(8πD_{tr}D_{c}) → ζ_{tr}, the rescaled rotational diffusiophoretic parameter
ζ_{rot}kh/(8πD_{c}√

D_{tr}D_{rot}) →ζ_{rot}, and the area fraction σ defined as the projected area of
all particles divided by the area of the simulation box. Note that the factor kh/(4πDc)
from Eq. (4.7) is subsumed into the rescaled diffusiophoretic parameters, when using

∇c_{2d}(r_{i}) in the Langevin equations (4.3) and (4.4).

### 4.3.3. Screening, repulsive particle-particle and particle-wall interactions

The particles with radius a interact via a hard-core repulsion. Whenever they overlap during the simulations, we separate them along the line connecting their centers. We have also tested a hard-core potential of the form U =

i

j(2a/|ri−rj|)^{36} revealing
essentially the same results.

When colloids form compact clusters, the chemical substance cannot diffuse freely be-tween the particles. Therefore, we implement a screened chemical field as follows.

Whenever a colloid is surrounded by six closely packed neighbors with distances below
r_{s} = 2a(1 +ϵ), we replace the term 1/|r−r_{i}| in Eq. (4.7) by

exp[−(r−ξ)/ξ]

|r−r_{i}| , (4.8)

where we introduce the screening length ξ := r_{s}. We set ϵ = 0.3 and checked that our
results do not change, if ϵ varies by 50%. In Ref. [40] screening occurs if the reaction
rate k depends on concentration c. However, as we have discussed before, this does not
apply to the experimental situation of Ref. [14], which we model.

case cluster aggregation emerges at the walls. This, however, has not been observed in experiments, for which reason we chose the random boundary condition. A more detailed description of the wall-particle interaction could be an interesting extension of our work.

### 4.3.4. Numerical implementation

We implemented a numerical solution of the Langevin equations in a two-dimensional square simulation box. To study different area fractions σ, we always use 800 particles and adapt the box size. In two dimensions, the moving direction is controlled by one angle, which we call φ. With e = [cos(φ),sin(φ)] the equations (4.3) and (4.4) are integrated by an Euler scheme:

φ_{i}(t+dt) = φ_{i}(t)−dtζ_{rot}

−sin(φ),cos(φ)

·∇c(r_{i}) +

2dtD_{rot}N_{rot},
r_{i}(t+dt) = r_{i}(t) +dt

v_{0}e_{i}−ζ_{tr}∇c(r_{i})
+

2dtD_{rot}N_{tr}

where we made use of equation (2.51). After updating the position r_{i} for each particle
i, we check for overlaps with other particles or the wall. The hardcore interaction that
we apply acts only on particles with distance |r_{i} − r_{j}| < 2a, i.e., it is strictly
short-ranged. The systems we study contain typically large clusters, which implies many
particle-particle interactions. To resolve them accurately, we apply a small time step of
dt = 5·10^{−6}sec, which makes the hard-core interactions the numerically most expensive
part of the program. Naively coding, for each time step the hard-core interactions have
numerical complexity O(N^{2}). To reduce it, we calculate so-called neighbor lists. To this
end, we subdivide the box in smaller boxes with width 2a and assign each particle to
one of these boxes [with complexity O(N)]. Then, neighbors for a particle i are in one
of the 8 neighboring boxes, each of which maximally contains one particle (given that
there are no overlaps). Hence, the total numerical complexity to search for possibly
overlapping particles is O(N + 8N) =O(N)< O(N^{2}). For the ∇c(r_{i}) term, we need to
consider the contribution of each particle j (see Eq. 4.7) due to the long-ranged nature
of the interactions. Hence, the numerical complexity is necessarily O(N^{2}). However, we
perform the update of c(r_{i}) only every 50th time step saving computational time again.

Figure 4.2.: Full state diagram ζ_{tr} versus ζ_{rot} at Pe = 19 and density σ = 0.05. More
details are given in the main text.