• Keine Ergebnisse gefunden

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

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 H2O2 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 v0 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 v0 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:

vD =−ζtr∇c. (4.1)

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

ωD =−ζrotei×∇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 ri and orientation ei 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 ˙ei = ωD ×ei. Then, we use Eqs.. (4.1) and (4.2) together with (ei×∇c)×ei = (1−ei⊗ei)∇c to write down the overdamped Langevin equations for each colloid:

i = v0ei−ζtr∇c(ri) +ξi, (4.3) e˙i = −ζrot(1−ei⊗ei)∇c(ri) +µi×ei. (4.4) Here, ξi is translational and µi rotational Gaussian white noise with zero mean and respective time correlation functions⟨ξi(t)⊗ξi(t)⟩= 2Dtr1δ(t−t) and⟨µi(t)⊗µi(t)⟩= 2Drot1δ(t−t), where we have introduced the translational (Dtr) and rotational (Drot)

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) = Dc2c−k

N

i=1

δ(r−ri), (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), c3d(r) = c0−(k/4πDc)N

i=11/|r−ri|, wherec0is 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:

c2d(r) = hc0− kh 4πDc

N

i=1

1

|r−ri|. (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 vD 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 lr = 

Dtr/Drot = 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] lr = 1.79a is measured. For historical reasons, we use here lr = 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 ζtrkh/(8πDtrDc) → ζtr, the rescaled rotational diffusiophoretic parameter ζrotkh/(8πDc

DtrDrot) →ζ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

∇c2d(ri) 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 rs = 2a(1 +ϵ), we replace the term 1/|r−ri| in Eq. (4.7) by

exp[−(r−ξ)/ξ]

|r−ri| , (4.8)

where we introduce the screening length ξ := rs. 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(ri) +

2dtDrotNrot, ri(t+dt) = ri(t) +dt

v0ei−ζtr∇c(ri) +

2dtDrotNtr

where we made use of equation (2.51). After updating the position ri 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 |ri − rj| < 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−6sec, 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(N2). 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(N2). For the ∇c(ri) 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(N2). However, we perform the update of c(ri) 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.