© The Author(s) 2024. Published by Science Press and EDP Sciences.

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-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-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.

MODEL

Single particle dynamics

Any overdamped colloidal particle, labeled here by the index i, executes a regular 2D Brownian motion with Langevin equation(1)where r_i=(x_i, y_i) are the coordinates of its center of mass and is a stationary source of Gaussian noise with and 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₀, or (ii) short-range pair-repulsions with truncated Lennard-Jones potential [19],where r_ij is the distance between particles i and j, r_m=2^1/6σlocates the potential minimum, and σr=2r₀ 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 [10, 8]. Without entering the details of the specific sensing mechanisms, we can define the sensing function of particle i as [14](2)where 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_ij≤d_c (see Figure 1A).

Figure 1 Coexisting phase patterns in a suspension of N=470 hard discs of radius r0=1 in a square box of side L=100. (A) Schematics of the QS protocol for the hot-cold transition with DM=0.1 (red dots) and Dm=0.01 (blue dots); (B)–(H) suspension snapshots at t=105 for Lp=10 and increasing QS distances, respectively, dc=5, 8, 10, 12, 15, 17, and 19.

The particle diffusivity is then governed by the following simple quorum sensing protocol:(3)where D_M>D_m and the transition threshold is defined as P_th=ρ₀L_p, with ρ₀=N/L² 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 . 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, , , and with , shows that the free dynamical parameters for our suspension are L_p/L, d_c/L, DM/D_m, and σ/L. Accordingly, the quantities defining the QS protocol scale like , , and . 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. (2), where it plays the role of a natural cut-off when . On the other hand, for sensing distances, d_c, larger than the mean interparticle distance, , 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, and, therefore, the packing fraction , 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₀ 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 and , where the QS protocol becomes respectively to discreteness (finite interparticle separation) and periodicity effects (sensing area encroaching the box boundaries).

Figure 2 (A) Phase diagram in the space parameter (dc, Lp) for the hard-discs suspension of Figure 1. Seven distinct phase topologies are distinguishable, say, for fixed Lp and increasing dc: 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 Lp, i.e., of the QS threshold in Eq. (2) (see text). (B) Fraction of cold particles, Nm/N, vs. dc for Lp=10. All remaining parameters are as in Figure 1. Vertical colored bands in (B) locate the different phase topologies from (A).

Under stationary conditions, that is for sufficiently long simulation runs (typically t>10⁵), the number of cold particles fluctuates around a stable average value, N_m, also displayed in Figure 2Bvs. 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 dc~L_p. This result comes as no surprise in view of the definition of the QS threshold in Eq. (3), P_th=ρ₀L_p, and recalling that the particles’ average sensing function for a uniform distribution would be .

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 , 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.

Figure 3 (A) Radial distribution function g(r) for the hard-disc suspension of Figure 2B with different values of dc (see legend in (C)). (B) Maximum amplitude of the g(r) curves in (A), gmaxvs. dc. 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=105. The local densities, ρ in (C) have been computed over square bins of side δL=10.

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⁵. We then averaged the bin density, ρ=δN/δL² 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 ρ₀=N/L² (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 dc/L_p fixed. Vice versa, this scaling property fails for exceedingly small () or large L_p values (), where discreteness and finite box size effects come into play.

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

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 : 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_maxvs. 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.

Figure 5 Dependence of (A) Nm/N and (B) gmax on the diffusivity parameters, DM and Dm. Other suspension parameters: N=382, Lp=10, and particles interacting via the potential of Eq. (2) with σ=2 and .

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, , 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(4)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 at d_c=0 down to 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 and are smaller than the relevant free-particle diffusivity, D_M and D_m, but their ratio remains the same, . 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 [10, 8, 24], with , which holds for φ≪1.

Figure 6 Effective particle diffusivity, Deff, vs. dc in a two-phase suspension with N=382, Lp=10, DM=0.1, Dm=0.01, and pair potential of Eq. (2) with σ=2 and . Single-particle Deff data (crosses) are compared with better statistics ensemble averaged data (dots) and the fitting law of Eq. (5) (solid curve), obtained from the corresponding Nm/N vs. dc curve in Figure 5. Inset: MSD vs. t for different phase topologies.

In the intermediate regime, d_c~L_p, for a stationary two-phase suspension a simple two-state argument suggests that with N_M=N-N_m; hence(5)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-at-distance 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(6)

In this case we can expect an anti-associative behavior with all individuals trying to distance themselves from one another so as to avoid local gatherings. This variation of the present model is the subject of ongoing investigations.

Funding

Y.L. was supported by the National Natural Science Foundation of China (12375037 and 11935010).

Author contributions

Y.Z. collected the data; Y.L. and F. M. conceived and analyzed the work; F.M. wrote the original draft. All authors have read and are agreed to publish the manuscript.

Conflict of interest

The authors declare no conflict of interest.

References

All Figures

	Figure 1 Coexisting phase patterns in a suspension of N=470 hard discs of radius r0=1 in a square box of side L=100. (A) Schematics of the QS protocol for the hot-cold transition with DM=0.1 (red dots) and Dm=0.01 (blue dots); (B)–(H) suspension snapshots at t=105 for Lp=10 and increasing QS distances, respectively, dc=5, 8, 10, 12, 15, 17, and 19.
In the text

Figure 2 (A) Phase diagram in the space parameter (dc, Lp) for the hard-discs suspension of Figure 1. Seven distinct phase topologies are distinguishable, say, for fixed Lp and increasing dc: 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 Lp, i.e., of the QS threshold in Eq. (2) (see text). (B) Fraction of cold particles, Nm/N, vs. dc for Lp=10. All remaining parameters are as in Figure 1. Vertical colored bands in (B) locate the different phase topologies from (A).

In the text

Figure 3 (A) Radial distribution function g(r) for the hard-disc suspension of Figure 2B with different values of dc (see legend in (C)). (B) Maximum amplitude of the g(r) curves in (A), gmaxvs. dc. 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=105. The local densities, ρ in (C) have been computed over square bins of side δL=10.

In the text

	Figure 4 Cold particle fraction, Nm/N, for the hard-disc suspension of Figure 1, but different Lp (solid symbols). Crosses represent the results obtained for Lp=10, i.e., same parameters as in Figure 2B, except here particles interact via the potential of Eq. (2) with σ=2 and .
In the text

	Figure 5 Dependence of (A) Nm/N and (B) gmax on the diffusivity parameters, DM and Dm. Other suspension parameters: N=382, Lp=10, and particles interacting via the potential of Eq. (2) with σ=2 and .
In the text

Figure 6 Effective particle diffusivity, Deff, vs. dc in a two-phase suspension with N=382, Lp=10, DM=0.1, Dm=0.01, and pair potential of Eq. (2) with σ=2 and . Single-particle Deff data (crosses) are compared with better statistics ensemble averaged data (dots) and the fitting law of Eq. (5) (solid curve), obtained from the corresponding Nm/N vs. dc curve in Figure 5. Inset: MSD vs. t for different phase topologies.

In the text