Swarm Hunting and Cluster Ejections in Chemically Communicating Active Mixtures

A large variety of microorganisms produce molecules to communicate via complex signaling mechanisms such as quorum sensing and chemotaxis. The biological diversity is enormous, but synthetic inanimate colloidal microswimmers mimic microbiological communication (synthetic chemotaxis) and may be used to explore collective behaviour beyond the one-species limit in simpler setups. In this work we combine particle based and continuum simulations as well as linear stability analyses, and study a physical minimal model of two chemotactic species. We observed a rich phase diagram comprising a “hunting swarm phase”, where both species self-segregate and form swarms, pursuing, or hunting each other, and a “core-shell-cluster phase”, where one species forms a dense cluster, which is surrounded by a (fluctuating) corona of particles from the other species. Once formed, these clusters can dynamically eject their core such that the clusters almost turn inside out. These results exemplify a physical route to collective behaviours in microorganisms and active colloids, which are so-far known to occur only for comparatively large and complex animals like insects or crustaceans.

x =x x u , t =t t u (1) Here, all the quantities with a tilde are dimensionless while those with a subscript " u " contain the units. We also introduce the dimensionless chemical field,c p,h = c p,h x d u , where d denotes the spatial dimension (d = 2 in our case). Eq. (2) in the main text then reads in dimensionless form as where we have used that δ(αr) = 1 α d δ(r). Next, we choose x 0 = R, where R denotes the (soft) radius of a particle and t u = 1 k0 yielding: where the normalised coefficientsD c and µ are then given by:D c = Dc k0R 2 and µ = k d k0 .

B. Particle ODEs
The same normalization is applied to Eq. (1) in the main text which leads to: where η s i (t) represents unit-variance Gaussian white noise with zero mean. With the additional definitions ofα kl = α kl γk0R d+2 andD = D k0R 2 we obtain: We model V using the Weeks-Chandler-Anderson potential, withṼ = 1 2 i,j =iṼ ij , where the sums run over all particles and whereṼ Here r ij denotes the distance between particle i and j,r ij = r ij /x u ,σ = σ/x u and˜ = γk0R 2 .
In the following, we omit tildes. The specific value of hardly affects our results, so the following independent and relevant dimensionless parameters remain: 1) effective diffusion constants D c and D, 2) reduced coupling constants α kl , 3) ratio of decay and production of the chemical µ 4) density of species p,h (which is conserved throughout the dynamics)

II. DESCRIPTION OF NUMERICAL METHODS
A. Numerical methods for solving the equations of the particle-based chemotaxis model We solve our model Eqs. (3,6,7) by discretizing space and time. We perform simulations in two spatial dimensions (2D) with periodic boundary conditions. Interactions between the particles are described using the Weeks-Chandler-Anderson potential. The PDEs (3) for the dynamics of the chemical fields are solved in time with a forward Euler method. This does not impose any restriction on the timestep dt due to the Courant-Friedrich-Lewy (CFL) condition [1] of the heat equation (D c dt dx 2 + dt dy 2 < 1 4 , where dt and dx, dy denote the temporal and spatial discretisations), since the necessary timestep dt for the particle dynamics is small enough to not violate the above CFL condition. The Laplace operator of the heat term is approximated with central differences. The particle dynamics of the colloids (Eqs. 6,7) is solved in time with the Euler-Maruyama scheme [2] for incorporation the noise term. The interactions of the colloids is effectively treated by applying a cell-list summation [3].
In order to speed up the calculations, we make use of cell lists [3], with a minimal image convention. This is motivated by the fact that the repulsive force caused by the potential V is short ranged and hence we can introduce a cutoff radius r c = 2 1/6 σ. More explicitly, we consider only interactions within a cutoff radius r c such that the pair interaction V ij = 0 for r ij > r c . Different simulations are performed by varying the diffusion constants D and D c , the coupling constants α kl and the ratio of decay and production of the chemical µ. In addition to the particle-based simulations, we have also numerically solved the system of the four coupled dimensionless PDEs of the continuum model, Eqs. (5,6) of the main text as before in two spatial dimensions with periodic boundary conditions. We use a forward Euler method to propagate the densities ρ s and the chemical fields c s in time. The Laplace operator is approximated with central differences.

Repulsion between cells
We must be careful when implementing the chemotactic continuum model, which unlike our particle based equations, neglect the steric repulsions among the particles, preventing sharp density peaks, destabilizing the simulation. In view of this effect, we add an additional nonlinear term to the continuum model to effectively describe local repulsions among the particles [4]: Here, D s rep with s ∈ {p, h} are dimensionless repulsion coefficients and we assume that the strength of the repulsion is the same between the two different species, D p rep = D h rep .

Regularization of the density
As a second adjustment of the continuum model, avoiding a blow up [5] we use: where κ is a small regularization parameter and κ → 0 leads to the original system. In [5] it was shown that the solutions of equations of the kind of Eq. (5) of the main text typically have a spiky structure and Eq. (10) can be used to model the aggregation phenomena.

III. LINEAR STABILITY ANALYSIS
To understand the origin of the observed pattern formation, we perform a linear stability analysis of Eqs.
where the corresponding PDEs for the phoretic field c s = c s (r, t), s ∈ {p, h} are given by: For the linear stability analysis, we consider small perturbations from the homogeneous steady state of the particle density ρ s = ρ s 0 + δρ s and the phoretic field c s = c s 0 + δc s . The resulting linearised system of four equations describes the time-evolution of the deviations of the particle density and the chemical field from the homogeneous state ρ s = ρ s 0 and c s = ρ s 0 /µ.
Since the linearised equations Eqs. (14-16) do not depend explicitly on time, we perform a Fourier transformation in space and make a separation Ansatzρ s (q, t) = e λtρs (q),ĉ s (q, t) = e λtĉs (q), leading immediately to the following eigenvalue problem (we set α hh = 0 as in the main text): where the stability of the system is determined by the eigenvalues λ of the matrix, as follows: (a) The steady state is stable, if the eigenvalues of the matrix all have real parts strictly less than zero.
(b) The steady state is unstable, if at least one of the eigenvalues of the matrix has a positive real part.
(c) Otherwise in the marginal case higher order terms determine the stability of the problem.
The four eigenvalues λ of the linearised system are explicitly given by (α hh = 0, For the parameter range under consideration, eigenvalues with non-vanishing imaginary part exist (Im(λ) = 0), if the condition −4α ph α hp > α 2 pp is fulfilled. In order to develop an instability for long wavelengths (|q| −→ 0), a Taylor series extension of (18) and an examination where Re(λ) > 0 provides us with the criterion If the eigenvalues are real (here, Im(λ) = 0 ⇐⇒ −4α ph α hp < α 2 pp ), then we obtain the following criterion for instability from calculating Re(λ) > 0: Altogether, an instability is given, when the following criterion is fulfilled: and this instability is also oscillatory if If the local repulsions among the particles is taken into account and the linear stability analysis is performed with Eq. (9) instead of Eq.(5) of the main text, the uniform phase is unstable, when the following criterion is fulfilled: which is then also oscillatory if Since the term µD rep is very small for the considered parameters in Fig. 2 (of the main text), the instability criterion is barely influenced by this additional D rep -term approximately leading to Eqs. (21,22) again.
Note that the instability criterion, Eq. (21), does not depend on the chemical diffusion constant D c . Therefore, a faster alternative to obtain Eq. (21), would be to assume that the dynamics of the chemical is fully enslaved by the motion of the colloids and relaxes quasi-instantaneously to its steady state, i.e. to assume ∂ t c s = 0. Note however, that only the instability criterion itself is independently of D c and not the instability band. In addition, once the instability has emerged, the fact that ∂ t c s = 0 is of crucial importance for the phenomenology of the resulting patterns, which can be seen best for the clusters ejecting their inner particles, which hinge on chemical delay (or memory) effects.
We have numerically tested the instability criterion in the parameter regime where it is oscillatory (i.e. non-vanishing imaginary part) to see if it is shifted due to perturbation convection, see [6]. Snapshots of the additional performed simulations are shown below in Fig. 1. The snapshots are arranged in a row with increasing decay-rate of the chemical µ. Above a certain value, this rate is too large, such that the concentration of the chemical substances is insufficient and the homogeneous phase is stable. The dashed line indicates for which µ the instability criterion (Eq. (7)) is fulfilled. As can be seen, this transition line fits well to the simulations and we did not find any shift, suggesting that the advective and absolute instability are very close to each other in the present case. In order to distinguish between the different phases, in particular between the clusters and the hunting swarm phase, we calculate the mean particle velocity, which is defined as where we average over all particles of species s ∈ {p, h} and · denotes the ensemble average.

V. DIFFUSION EQUATION IN 3D
To verify that our results do not change qualitatively when solving the diffusion equation in 3D, we have performed additional simulations for verification. Figure 2 shows an exemplaric simulation snapshot of the particle based model in three dimensions for a cluster in the blue domain (α hp = 0.01, α ph = −10).
FIG. 2: Simulation snapshot of the particle based model in 3D. Left panel shows the diffusion of c p by looking at the cross section of the chemical field in a plane perpendicular to the z-axis. Right panel shows the prey and the hunters, colored in black and red, respectively. Parameters as in Fig. 2(b) of the main text.

VI. DESCRIPTION OF THE MOVIES
The movies 1-8 show the time-evolution of the different patterns, hunting swarms, mixed clusters and core-shell clusters, for the particle-based and the continuum model for the same parameters as used in Fig. 2 of the main text (panels (a-h) corresponding to the files mov1-mov8). Movie 9 shows a dynamical inside out reversal for an appropriate initial condition, starting with a dense and comparatively large cluster of prey particles, which is then invaded by the hunter particles. Parameters used for movie 9: α pp = 1, α hh = 0, α hp = 10, α ph = −0.01, µ = 0.001, D c = 0.1, D = 0.003, = 1, 2N = 300 particles and L box = 80. Movie 10 shows clusters ejecting their inner particles that occur starting from a uniform initial state. Parameters used for movie 10: α pp = 100, α hh = 0, α hp = 10, α ph = −1000, µ p = 0.1, µ h = 0.01, D p c = 0.5, D h c = 10, D = 0.001, = 10, L box = 200 with 500 prey-particles and 2000 hunter-particles.