Clustering of quorum sensing colloidal particles

: We propose a simple model of colloidal suspension, whereby individual particles change their di ﬀ usivity from high (hot) to low (cold), as the local concentration of their closest peers grows larger than a certain threshold. Such a non-reciprocal interactive mechanism is known in biology as quorum sensing. Upon tuning the parameters of the adopted quorum sensing protocol, the suspension is numerically shown to go through a variety of two-phase (hot and cold) conﬁgurations. This is an archetypal model with potential applications in robotics and social studies.


INTRODUCTION
Our understanding of how the microscopic details of the single-particle dynamics lead to different collective behaviors is presently far from satisfactory.For this reason, of late practitioners often resort to non-reciprocal interactions as a tool to model the autonomous clustering of active matter [1][2][3].By non-reciprocity we refer to the situation when one particle perceives the presence of a second particle, which, therefore, influences the dynamics of the first particle without modifying its own.Being typically long-ranged, such non-Newtonian interactions involve neither fields of force [4][5][6], nor (reciprocal) pair-potentials [7], nor steric effects [8,9], as assumed, for instance, in the early literature on motility-induced phase separation (MIPS) [10].Indeed, finite-size self-propelled particles may form aggregates by adjusting their velocity according to the direction and the local density of their peers, a mechanism known in biology as quorum sensing (QS) [11,12].
In its simplest form, MIPS has been shown to be the analogous of a gas-liquid phase separation and can be observed both in biological and synthetic active matter.However, there is still no consensus on how to induce these phases and numerical studies report different behaviors depending on model details [13].
Elementary models of QS active suspensions have been numerically demonstrated to exhibit intriguing clustering capabilities depending on the encoded sensing protocol.Bechinger and coworkers [14] showed that disordered swarms can form when the self-propulsion of individual particles switches off either in regions of high peer concentration, independently of their relative orientation (active-to-passive transitions), or, on reverse, when their sensing cone points away from high density regions (passive-to-active transitions) [15].Chiral QS protocols can be invoked to generate rotating swarms, or swirls [16].Moreover, QS particles switching between two active states have been reported to exhibit either social or antisocial behaviors, depending on whether their self-propulsion speed grows smaller or larger for densities above a given concentration threshold (active-to-active transitions) [17].
Inspired by the recent literature on QS active matter, we propose here an even simpler QS protocol, which applies to regular colloids, that is to passive particle suspensions.We assume that the diffusion constant of the overdamped suspension particles can switch from high to low, as the peer density they perceive in their closest surroundings raises above a set threshold.Under certain conditions, such a QS protocol can trigger a clustering process, which eventually leads to a dynamical phase separation.
The content of this paper is organized as follows.In Section "Model", we detail our QS protocol, characterized by a tunable sensing distance and transition threshold.In Section "Clustering via quorum sensing" , we analyze quantitatively our numerical results for the clustering of a two-dimensional (2D) suspension of hard discs, by introducing appropriate quantifiers.In Section "Model parameter dependence", we discuss the robustness of the proposed QS clustering mechanism against changes of the model parameters.In Section "Effective diffusivity in a two-phase suspension", we investigate the effective diffusivity of the suspension particles to stress the dynamical nature of the coexisting phases.Finally, in Section "conclusions", we briefly mention variations of the present model and possible applications beyond colloids.

Single particle dynamics
Any overdamped colloidal particle, labeled here by the index i, executes a regular 2D Brownian motion with Langevin equation where r i = (x i , y i ) are the coordinates of its center of mass and ξ i (t) = [ξ ix (t), ξ iy (t)] is a stationary source of Gaussian noise with ⟨ξ iq (t)⟩ = 0 and ⟨ξ iq (t)ξ ip (0)⟩ = 2D i δ qp δ(t) for q, p = x, y, modeling the equilibrium thermal fluctuations in the suspension fluid.The particle diffusivity, D i , is allowed to switch between a higher and lower value through the QS protocol detailed below.The stochastic differential Eq. ( 1) was numerically integrated by means of a standard Euler-Maruyama scheme [18].To ensure numerical stability, the numerical integrations have been performed using an appropriately short time step, δt = 10 −3 (See the figure captions for more numerical details).

Colloidal suspension dynamics
We numerically simulated a colloidal suspension by placing N identical, independent colloidal particles of Eq. ( 1) in a square box of side L and imposing periodic boundary conditions.In a suspension, effects due to the particle finite size cannot be neglected.We modeled steric interactions as either: (i) hard-disk collisions, whereby the particles were modeled as hard discs of radius r 0 , or (ii) short-range pair-repulsions with truncated Lennard-Jones potential [19], where r i j is the distance between particles i and j, r m = 2 1/6 σ locates the potential minimum, and σ = 2r 0 is the particle effective diameter.Further reciprocal interactions and hydrodynamical interactions [20,21] have been neglected.

Quorum sensing protocol
We then assumed that the diffusivity of each particle depends on the spatial distribution of its neighbors.In biological systems this process is mediated by some form of inter-particle communication (mostly chemical in bacteria colonies [11,12]).On the other hand, the diffusivity of colloids is suppressed with increasing density [8,10].Without entering the details of the specific sensing mechanisms, we can define the sensing function of particle i as [14] where V d c i denotes its perception area of radius d c .This means that each suspension particle senses the presence of its peers in all directions within a restricted horizon, r i j d c (see Figure 1A).
The particle diffusivity is then governed by the following simple quorum sensing protocol: where D M > D m and the transition threshold is defined as P th = ρ 0 L p , with ρ 0 = N/L 2 denoting the suspension density in the case of uniform spatial distribution.Note that for a uniform suspension the particles' sensing function of Eq. ( 2) is approximatively P(d c ) = ρ 0 d c .Clearly, this form of particle interaction is non-reciprocal, since j may be perceived by i without being itself affected by the presence of i.

Tunable model parameters
A simple rescaling of the space and time variables, x → x/L, y → y/L, and t → t/τ with τ = L 2 /D m , shows that the free dynamical parameters for our suspension are  2), where it plays the role of a natural cut-off when r i j → 0. On the other hand, for sensing distances, d c , larger than the mean interparticle distance, l N = √ L 2 /N, QS is expected to control the suspension diffusivity irrespective of any short-ranged reciprocal interactions.For this reason, in our simulations we fix the suspension parameters, N and L, the particle size, σ = 2r 0 and, therefore, the packing fraction ϕ = πN(σ/L) 2 , and vary the remaining tunable parameters L p , d c , and D M /D m .

CLUSTERING VIA QUORUM SENSING
We start presenting our numerical results for a suspension of N hard discs of radius r 0 confined to a square box of side L. As illustrated in Figure 1, on increasing the ratio d c /L p the suspension separates into two phases, a hot and a cold one respectively with D i = D M and D i = D m .For d c /L p ≪ 1 (d c /L p ≫ 1), all particles turn hot (cold).As d c /L p approaches 1, we observe the appearance of several small cold clusters, which eventually merge into one large cold cluster.Analogously, upon lowering d c /L p toward 1, small cavities open up in the uniform cold phase, which host sparse hot particles.As d c /L p approaches 1, these cavities also merge into one large cavity.One can regard the cavity formation in the cold phase as a sort of symmetric counterpart of the cluster formation in the hot phase.It is worth noticing that for d c /L p ∼ 1 the transition between cluster and cavity phase separation is marked by the emergence of band structures that cut across the entire simulation box.This is a symmetry breaking effect induced by the periodic conditions at the box boundaries.The orientation of the emerging bands can be either vertical or horizontal, depending on the system initial conditions.The coexistence of extended compact hot and cold phases would require an infinitely large simulation box.
The corresponding phase diagram in the (d c , L p ) plane is shown in Figure 2A for a fixed ratio D M /D m .The pattern sequence of Figure 1 repeats itself for any value of L p , except L p l N and L p L/2, where the QS protocol becomes respectively to discreteness (finite interparticle separation) and periodicity effects (sensing area encroaching the box boundaries).
Under stationary conditions, that is for sufficiently long simulation runs (typically t > 10 5 ), the number of cold particles fluctuates around a stable average value, N m , also displayed in Figure 2B vs. d c at fixed L p .The corresponding topology of the coexisting phases is marked by the same color code as in Figure 2A for reader's convenience.As already apparent in Figure 2A, the transition from a predominantly hot to a predominantly cold suspension occurs for d c ∼ L p .This result comes as no surprise in view of the definition  of the QS threshold in Eq. ( 3), P th = ρ 0 L p , and recalling that the particles' average sensing function for a uniform distribution would be P(d c ) = ρ 0 d c .
To further characterize the QS induced hot-to-cold transition, we considered the suspension-averaged radial distribution function g(r) [22].A few g(r) curves for the hard disc suspension of Figures 1 and 2B are displayed in Figure 3A for different values of the QS distance, d c .The strength of the QS effect is quantified by the height, g max , of the first g(r) peak around r σ, reported in Figure 3B as a function of d c .The clustering effect turns out to be the strongest as the one cold cluster grows into a periodic band, namely, for the two-phase topology that separates the predominantly hot from the predominantly cold configurations.Recalling that for a uniform particle distribution g(∞) = 1, we are not surprised to observe that in the presence of clustering, d c < L p , the tails of g(r) approach asymptotic values g(∞) < 1. Vice versa for d c > L p , the formation of hot cavities increases the density of the predominant cold phase, so that g(∞) > 1.
Finally, we also considered the distribution of the local suspension density ρ.We divided the simulation box into square bins of size δL × δL and counted the number of particles, δN, in each bin under stationary conditions, i.e., for t > 10 5 .We then averaged the bin density, ρ = δN/δL 2 over time for better statistics.Some ρ distributions, P(ρ), are plotted in Figure 3C.The coexistence of two phases is signaled by bimodal ρ distributions, with peaks to the left (right) of ρ 0 = N/L 2 (uniform particle distribution) respectively for the hot (cold) phase.The two peaks equilibrate in the intermediate band phase topology.For a more sophisticated description of the clustered suspension we implemented Voronoi's tessellation to consistently partition the simulation box into non-overlapping local polyhedra without introducing the arbitrary mesh length δL [23].
The area distribution of the Voronoi cells (not shown) proved consistent with the curves of Figure 3C.

MODEL PARAMETER DEPENDENCE
To establish the robustness of the QS induced phase transition we discuss here the dependence of the quantifiers introduced above on the model parameters.In Figure 4 we focus on the dependence of the average cold disc fraction, N m /N, on the QS threshold, P th , i.e., on the parameter L p .As anticipated in Section "Clustering via quorum sensing", the simulation output is insensitive to the actual value of L p as long as we keep the ratio d c /L p fixed.Vice versa, this scaling property fails for exceedingly small (L p l N ) or large L p values (L p L/2), where discreteness and finite box size effects come into play.
We then replaced the disc hard-core collision with particle pair interactions mediated by the truncated Lennard-Jones potential of Eq. (2) of equal (effective) diameter σ = 2r 0 : no appreciable changes in the two phase fractions were detected.In the following we will keep reporting data for suspensions of particles interacting via the potential of Figure 2 to further substantiate this conclusion.
As anticipated in Section "Model", we expect that under stationary conditions QS induced diffusivity phase transitions only depend on the ratio D M /D m , one diffusivity constant, say D m , entering the time scaling factor τ. Our expectations are confirmed by the simulation data for N m /N and g max vs. d c plotted in Figure 5. Curves for the same ratios D M /D m overlap, no matter what the choice of D M and D m .It is worth pointing out that on increasing D M /D m the hot-to-cold transition is slightly anticipated (Figure 5A), and the cluster packing enhanced (Figure 5B).This is an effect of the diffusivity difference between coexisting phases, as also detected in hard-disc suspensions.Finally, we looked at the dependence of the curves N m /N versus d c on the suspension parameters N and σ (not shown).Due to local density fluctuations, the onset of the QS transition happens for sensing distances, d c , appreciably lower than L p (see Figure 4).However, we noticed that the transition onset shifts to higher d c values, d c → L p −, with increasing N.Moreover, the transition onset distance also decreases (almost proportionally) with σ.This effect is due to the fact that, as anticipated in Section "Model", the particle size enters the QS protocol as a natural cut-off of the sensing function of Eq. ( 2).These remarks concern the transition onset alone, the overall hot-to-cold transition being localized by the condition d c ∼ L p , as discussed above.

EFFECTIVE DIFFUSIVITY IN A TWO-PHASE SUSPENSION
We conclude our analysis by investigating the effective diffusivity of a particle in a two-phase suspension.Of course, hot and cold particles have instantaneous diffusivity D M and D m , respectively.The question then rises of how they diffuse within the suspension.Our numerical simulations show that all particles obey the same normal diffusion law where ⟨. . .⟩ denotes a time average taken over the trajectory of a single particle.For a better statistics and without loss of generality, in Figure 6 the single-particle effective diffusivity has been ensemble averaged over all N suspension particles.When plotted versus d c , the single-particle D eff drops from D * M at d c = 0 down to D * m for d c >> L p .These two limiting situations refer to the hot and cold uniform one-phase suspensions, respectively.Subject to QS transitions, a single particle thus switches from hot to cold, and vice versa, without belonging permanently to either phase.Moreover, one notices by inspection that D * M and D * m are smaller than the relevant free-particle diffusivity, D M and D m , but their ratio remains the same, D * M /D * m = D M /D m .Indeed, we checked that the effective diffusivity in a uniform one-phase suspension decays with the packing fraction, ϕ, according to the well known power law [8,10,24], D * M(m) = D M(m) (1 − λϕ) 2 with λ ≃ 0.9, which holds for ϕ ≪ 1.
In the intermediate regime, d c ∼ L p , for a stationary two-phase suspension a simple two-state argument suggests that This equation relates the d c dependence of D eff and N m /N, in good agreement with the simulation data of Figure 6 (even for a single particle).

CONCLUSIONS
We have analyzed a simple model of QS induced phase transition, whereby the particles of a colloidal suspension switch from high to low diffusivity as the concentration of their peers in their immediate vicinity grows above a certain threshold.We have numerically proven that such mechanism suffices to trigger the formation of cold clusters (hot cavities) in a hot (cold) suspension.The coexistence of two phases under stationary conditions has been characterized by means of several quantifiers, like particle phase fractions, radial distributions, an effective diffusivity.In the current literature, like in the present paper, QS synthetic matter has been investigated in 2D, only, because no actual 3D experiments have been performed so far.However, experiments of living swimmers suggest that QS clustering also occurs in 3D [11,12].While in our presentation we refer for simplicity to the simulated system as an ideal colloidal suspension, the proposed model can be extended to describe assemblies of less elementary units, where sensing-atdistance is much easier to implement.For instance, it is everyday's experience that individuals in a social setting tend to slow down the pace as they perceive, for instance visually, the presence of an even casual gathering, thus triggering the formation of large, persisting crowds.Stated otherwise, the proposed clustering mechanism relies on non-reciprocal interactions, which can operate at much longer distances than any pair interaction, steric, electrostatic and hydrodynamic, alike.
By the same token, one can think of reversing the QS mechanism of Eq. ( 3), as

Figure 1
Figure 1 Coexisting phase patterns in a suspension of N = 470 hard discs of radius r 0 = 1 in a square box of side L = 100.(A) Schematics of the QS protocol for the hot-cold transition with D M = 0.1 (red dots) and D m = 0.01 (blue dots); (B)-(H) suspension snapshots at t = 10 5 for L p = 10 and increasing QS distances, respectively, d c = 5, 8, 10, 12, 15, 17, and 19.
L p /L, d c /L, D M /D m , and σ/L.Accordingly, the quantities defining the QS protocol scale like ρ 0 → L 2 ρ 0 = N, P th → LP th = N(L p /L), and P(d c ) → L P(d c ) = N(d c /L).The choice of the diffusion constants, D m and D M , affects the stationary properties of the suspension only through the ratio D M /D m .The particle diameter, σ, enters the simulation output mostly through the definition of the actual sensing function P i (d c ) in Eq. (

Figure 2
Figure 2 (A) Phase diagram in the space parameter (d c , L p ) for the hard-discs suspension of Figure 1.Seven distinct phase topologies are distinguishable, say, for fixed L p and increasing d c : all hot, small clusters, one cluster, bands, one cavity, small cavities, all cold.Configurations with many small clusters or cavities are hard to detect at too low or too large values of L p , i.e., of the QS threshold in Eq. (2) (see text).(B) Fraction of cold particles, N m /N, vs. d c for L p = 10.All remaining parameters are as in Figure 1.Vertical colored bands in (B) locate the different phase topologies from (A).

Figure 3
Figure 3 (A) Radial distribution function g(r) for the hard-disc suspension of Figure 2B with different values of d c (see legend in (C)).(B) Maximum amplitude of the g(r) curves in (A), g max vs. d c .The color bands denote the different phase topologies as in Figure 2B.(C) Distribution of the local suspension density for the same model parameters as in (A).All curves reported in (A)-(C) have been averaged over 1000 stationary configuration taken 100 time units apart, starting from t = 10 5 .The local densities, ρ in (C) have been computed over square bins of side δL = 10.

Figure 4
Figure 4  Cold particle fraction, N m /N, for the hard-disc suspension of Figure1, but different L p (solid symbols).Crosses represent the results obtained for L p = 10, i.e., same parameters as in Figure2B, except here particles interact via the potential of Eq. (2) with σ = 2 and ϵ = 1/6.

Figure 5
Figure 5 Dependence of (A) N m /N and (B) g max on the diffusivity parameters, D M and D m .Other suspension parameters: N = 382, L p = 10, and particles interacting via the potential of Eq. (2) with σ = 2 and ϵ = 1/6.

Figure 6
Figure 6 Effective particle diffusivity, D eff , vs. d c in a two-phase suspension with N = 382, L p = 10, D M = 0.1, D m = 0.01, and pair potential of Eq. (2) with σ = 2 and ϵ = 1/6.Single-particle D eff data (crosses) are compared with better statistics ensemble averaged data (dots) and the fitting law of Eq. (5) (solid curve), obtained from the corresponding N m /N vs. d c curve in Figure5.Inset: MSD vs. t for different phase topologies.