One of the most intriguing phenomena in active matter has been the gas–liquid-like motility-induced phase separation (MIPS) observed in repulsive active particles. However, experimentally, no particle can be a perfect sphere, and the asymmetric shape, mass distribution, or catalysis coating can induce an active torque on the particle, which makes it a chiral active particle. Here, using computer simulations and dynamic mean-field theory, we demonstrate that the large enough torque of circle active Brownian particles in two dimensions generates a dynamical clustering state interrupting the conventional MIPS. Multiple clusters arise from the combination of the conventional MIPS cohesion, and the circulating current caused disintegration. The nonvanishing current in non-equilibrium steady states microscopically originates from the motility “relieved” by automatic rotation, which breaks the detailed balance at the continuum level. This suggests that no equilibrium-like phase separation theory can be constructed for chiral active colloids even with tiny active torque, in which no visible collective motion exists. This mechanism also sheds light on the understanding of dynamic clusters observed in a variety of active matter systems.
INTRODUCTION
Active matter can spontaneously form structures not restricted by equilibrium thermodynamics as they keep dissipating energy and break the time-reversal symmetry locally. For example, while the gas–liquid-like phase separation in equilibrium requires cohesive interactions,1 active colloids can undergo motility-induced phase separation (MIPS) resulting from the combination of self-propulsion and steric repulsion. Essential physics of MIPS has been well captured by linear swimmer models such as active Brownian particles (ABPs),2,3 run-and-tumble particles (RTPs),4,5 and active Ornstein–Uhlenbeck particles (AOUPs).6 The coarse graining of those microscopic linear models usually restores the detailed balance at the continuum level, and MIPS was mostly understood based on equilibrium-like phase separations.7–9
Despite the success of linear swimmer models in capturing the essential physics of MIPS, they fail in describing the commonly observed chiral swimmers. Many biological micro-organisms, including bacteria,12,13 sperm cells, and some algae,14 swim in circular and helical trajectories near the surface or in chemical gradients.10,11 Those chiral swimming patterns play an important role in surface selection, attachment, and forming microcolonies of micro-organism.15 It is also convenient to design synthetic circle swimmers by an asymmetric shape,16–18 mass distribution,19 and catalysis coating20 to control the radius and frequency of circular trajectories. Individually, these chiral swimmers exhibit intriguing phenomena such as gravitaxis and21 abnormal transportation.22–24
The emergent behavior of circle swimmers with explicit alignments was studied by generalizing the Vicsek model, where slow rotations enhance the polarization in macroflocks, while fast rotations induce secondary instabilities, leading to phase synchronized microflocks.25 However, these flocks are integrated by phase-lock mechanism and suffer from strong fluctuations once including excluded volume interactions, which occur naturally among circle swimmers.26 In simulations of purely repulsive circle swimmers, a significant suppression of MIPS, counterclockwise rotation of macrodroplet at some optimal slow rotation frequencies,27 formation of the vortex array structure,28 and collective oscillation of the density under a strong chirality29 were observed. Recently, a field theory quantitatively accounting for the suppression of MIPS was proposed.30 However, as they consider the effect of self-propelling torque using an effective rotational diffusion coefficient, it is only valid in the slow rotation region. In contrast, non-equilibrium hyperuniform fluids were found in a deterministic chiral active particle model,31,32 in which the active rotation rate can be seen as extremely fast as the rotational noise is zero. Therefore, the unified understanding of the collective assembly of chiral active swimmers remains unknown. To this end, we formulate a hydrodynamic description from the microscopic dynamics of circle active Brownian particles (cABPs), and it is valid at both slow and fast rotations. We find a short wavelength instability leading to a dynamical clustering state at fast rotation, and the theoretically predicted phase boundary agrees well with large-scale Brownian dynamics simulations. By discussing the hydrodynamic matrix, we propose a novel instability mechanism originating from the combination of conventional MIPS cohesion and the circulating current that caused disintegration. We also verify that the system spanning circulating current originates from the motility “relieved” by fast rotation from temporal fluctuations. Our results may help understand the general mechanism of self-limiting size cluster formation without introducing explicit alignment to synchronize phases25 or chemical signaling.33,34 More generally, nonvanishing current in steady states is a distinct feature of the non-equilibrium systems, as not restricted by the thermodynamic balance. Our results suggest that no equilibrium-like theory, based on the detailed balance, can be constructed to understand chiral active colloids.
MODEL OF CIRCLE ACTIVE BROWNIAN PARTICLES
We consider N self-propelled particles in two dimensions interacting via the repulsive, pairwise additive, Weeks–Chandler–Andersen potential
with the cutoff at rc = 21/6σ, beyond which V = 0. Here, σ is the nominal particle diameter, ϵ determines the interaction strength, and r is the center-to-center distance between two particles. Particle i with positions ri and orientations ei = (cos θi, sin θi) evolves in response to a systematic driving force and rotational torque, according to the following overdamped Langevin equations:35
Here, μ is the translational mobility and v0 is the magnitude of the self-propulsion velocity. The Gaussian white noise ξi models the interaction with the solvent, but as the translational diffusion effect is much smaller than that of the self-propulsion, we set Dt = 0. To model the circular motion of active particles, the orientation θi displays a drift of the constant angular velocity ω0 alongside the rotational diffusion coefficient Dr, with νi being the unit Gaussian white noise ⟨νi(t)νj(t′)⟩ = δijδ(t − t′). In the dilute and deterministic limit, an individual particle moves in a circular trajectory with revolution radius R = v0/ω0. Then, there are three characteristic time scales: rotational diffusion τr = 1/Dr, circular motion τω = 2π/ω0, and ballistic motion . To construct the full phenomenology of this model requires scanning a four-parameter phase diagram with space and time units σ and 1/Dr, respectively, parameterized by the Pclet number Per = v0/(σDr), the dimensionless angular velocity Γ = ω0/Dr, the potential stiffness κ = μϵ/(v0σ), and the area fraction ϕ = πσ2N/(4L2) = ρ0πσ2/4 with L and ρ0 being the side lengths of the simulation box and density, respectively. Here, we set σ = 1 and μ = 1 and fix 1/κ = v0/ϵ ≡ 24 so that two head-to-head colliding cABPs reach the dynamic balance between their self-propulsion and pairwise interaction. In our Brownian dynamics simulations, periodic boundary conditions are applied in both directions, and the time step Δt satisfies v0Δt < 10−3σ. To ensure that the system reaches a steady state, we calculate the total potential energy of the system, which is ensured to be stable for at least 1000τb before the measurement of any quatity in the steady state.
RESULTS
Hydrodynamic description and finite wavelength instability
To understand the emergence of inhomogeneous states, with the mean-field approximation,36–38 we obtain the dynamic equations for the density and polarization fields ρ(r) = ∑iδ(r − ri) and W(r) = ∑iδ(r − ri)ei, respectively, as
Here, ve(ρ) = v0 − ζρ is the mean-field effective velocity with ζ reflecting how strongly the motility is slowed down by its neighbors, and it can be written as ve(ϕ) = v0(1 − ϕ/ϕ*) with ϕ* being the packing fraction damping ve to zero. is the unit vector indicating the direction of torque. Although this linear relationship was originally proposed for linear ABPs, we verify that it is still valid in homogeneous cABPs, and ϕ* does not depend on ω0 (see the supplementary material). D describes the effective diffusion caused by particle collisions.
Apparently, the isotropic homogeneous steady state H(ρ, W) = (ρ0, 0) is a solution for Eq. (3). The linearized dynamics of fluctuations are
where the Fourier transform is defined as ; is the polarization fluctuation longitudinal/transversal to q; the damping rates are Γρ = Dq2 and Γr = Dq2 + Dr; and the density and longitudinal polarization evolutions are coupled by the average effective velocity and the effective compressional modulus . Intuitively, the positive η drives the current to avoid the denser zone, while the negative η results in the current pointing to the denser zone. MIPS instability in the linear ABP models is usually intrigued by ∂ve/∂ρ < 0 and η < 0, suggesting that particles mobility slows down by the steric repulsion from its neighbors, and particles tend to accumulate where they move slower.9 Such slowdown-accumulation feedback loop provides an effective cohesion for phase separation, and we verify here that this mechanism is also preserved in chiral swimmers. However, the transversal component cannot be decoupled from the other two modes at the linear order as in normal fluids due to the convection induced by automatic rotation ω0.
Looking at the wave mode exp(λt + iq · r), the dispersion relation λ(q = |q|) unveils two types of instability of the homogeneous state depending on the relative strength of the angular velocity and rotational diffusion Γ = ω0/Dr. For slow rotation Γ < 1, an increase in ϕ and Per leads to the MIPS-like instability [Fig. 1(d)] when , where accounts for the suppression of MIPS by automatic rotation similar to Ref. 30, while for fast rotation Γ > 1, we find a qualitatively different picture, where an instability starts at q > q* > 0 [Fig. 1(e)] corresponding to the short wavelength fluctuation. This type II instability arises when . As in this regime, the MIPS instability is interrupted [with pseudo-phase boundary shown by the cyan plane in Fig. 1(a) and green dashed curves in Fig. 1(b)].
Theoretically predicted phase diagram of circle active Brownian particles. (a) Predicted phase boundaries for the control parameter [ϕ, Per, Γ] at Dr = 0.1 with the phenomenological mean-field parameter D = 151σ2Dr, ϕ* = 0.63. The blue plane separates the stable homogeneous state with type I instability, the red plane separates the stable homogeneous state with type II instability, and the cyan plane is the pseudo-separation plane for type I and type II instability. The green/black curves are the cross section of type I/II phase plane with the fixed Per plane. The yellow curve marks Γ = 1. (b) Cross section at Per = 240. (c)–(e) Dispersion relationship λ(q) for the stable homogeneous state (■), type I instability (▾), and type II instability (•).
Theoretically predicted phase diagram of circle active Brownian particles. (a) Predicted phase boundaries for the control parameter [ϕ, Per, Γ] at Dr = 0.1 with the phenomenological mean-field parameter D = 151σ2Dr, ϕ* = 0.63. The blue plane separates the stable homogeneous state with type I instability, the red plane separates the stable homogeneous state with type II instability, and the cyan plane is the pseudo-separation plane for type I and type II instability. The green/black curves are the cross section of type I/II phase plane with the fixed Per plane. The yellow curve marks Γ = 1. (b) Cross section at Per = 240. (c)–(e) Dispersion relationship λ(q) for the stable homogeneous state (■), type I instability (▾), and type II instability (•).
Dynamical clustering state
To specify different inhomogeneous states from two types of instability, we simulate the collective behavior of N = 40 000 identical cABPs with fixed Per and increasing ϕ, and the typical snapshots are summarized in Fig. 2(a). For slow rotations Γ ≲ 1, we reproduce the MIPS as in ABPs: the system remains homogeneous (state H) below the spinodal (), whereas induces a bulk phase separation with a single dense liquid droplet coexisting with a dilute gas phase (state P). At with fast rotation Γ ≳ 1, we observe the emergence of multiple dynamic clusters (state C), and they continuously merge, split, and decay but never merge into a stable dense bulk phase39,40 (see the supplementary material video). Although there are small temporary clusters emitted out from the bulk dense phase in MIPS, they are negligible compared with the system size. The finite size scaling is studied in Fig. 4 to characterize the states. The H–C phase boundary from computer simulation agrees with the linear instability predicted with the fitting parameter ϕ* = 0.58 and D = 461σ2Dr as shown by the red dashed curve in Fig. 2(b). The green dot marks Γ = 1, below which the predicted shown by cyan line deviates from the H–P phase boundary in simulation, and this may result from the large hysteresis loop in MIPS.
Simulated phase diagram of circle active Brownian particles. (a) Typical snapshots of three states, i.e., homogeneous state (H), phase separated state (P), and dynamical clustering state (C). Color of the particles indicates the self-propulsion orientation shown in the top-left inset. (b) Simulated phase diagram of cABPs, in which blue dots, black squares, and red triangles represent the homogeneous state (H), the MIPS state (P), and the dynamic clustering state (C), respectively. The red/cyan dashed curves are the theoretically predicted phase boundary between the stable homogeneous state and type II/I instability with the fitting parameter ϕ* = 0.58 and D = 461σ2Dr. Here, Dr = 0.1, and Per = 360. The yellow line is the H–C phase boundary obtained in simulations for Per = 1200 and Dr = 0.1.
Simulated phase diagram of circle active Brownian particles. (a) Typical snapshots of three states, i.e., homogeneous state (H), phase separated state (P), and dynamical clustering state (C). Color of the particles indicates the self-propulsion orientation shown in the top-left inset. (b) Simulated phase diagram of cABPs, in which blue dots, black squares, and red triangles represent the homogeneous state (H), the MIPS state (P), and the dynamic clustering state (C), respectively. The red/cyan dashed curves are the theoretically predicted phase boundary between the stable homogeneous state and type II/I instability with the fitting parameter ϕ* = 0.58 and D = 461σ2Dr. Here, Dr = 0.1, and Per = 360. The yellow line is the H–C phase boundary obtained in simulations for Per = 1200 and Dr = 0.1.
To further examine the phase boundary predicted in the hydrodynamic theory, we derive the slope of the H–C phase boundary in the parameter space [ϕ, Γ],
which suggests that the slope of type II instability increases with . Thus, we increase v0 by three times in simulation and do observe the gentle increase in the boundary (the red dashed curve) changing to nearly vertical, as shown by the yellow line in Fig. 2(b). The agreement between simulations and the instability analysis verify that the dynamical clustering state is induced by the short wavelength instability at fast rotation.
To eliminate the possibility that the dynamical clustering can be a metastable state or near-critical fluctuation of the bulk phase separation, we further simulate the system deep in the dynamical clustering state of [ϕ, Per] for Γ > 1, as shown in Fig. 3. Although clusters start emerging at ϕ ≃ 0.32, Dr = 0.1, and Per ≃ 240 with further increase in ϕ to 0.40 and Per by a factor of 5, these clusters remain highly dynamical instead of merging into a stable giant cluster as in MIPS. To ensure that the systems have reached the dynamical clustering steady states, we also perform simulations starting from MIPS configurations, in which the single giant droplet splits into multiple smaller clusters after switching on fast rotation. This suggests that the dynamical clustering is a well-defined non-equilibrium steady state.
As active matter systems generally suffer from strong finite size effects,41 we further examine the effect of the system size in two inhomogeneous states, as shown in Fig. 4. n(s) is the cluster size distribution (CSD) for the cluster containing s particles. Particles with distance less than rc are considered to be in the same cluster. The homogeneous state exhibits a CSD decaying in the form n(s) ∝ s−α exp(−s/sξ) with α ≃ 1.9 (see the supplementary material). In MIPS, a separated dot represents the single giant cluster [Fig. 4(a)], and the time averaged largest cluster size Smax increases linearly with the system size confirming the bulk phase separation scenario. For the dynamical clustering state, a pronounced bump emerges abruptly at an intermediate size s ∼ 103–104 after crossing . Remarkably, increasing the system size by eight times does not produce any visible change in the CSD of the dynamical clustering state. In addition, Smax increases much slower than ∝N for the dynamical clustering state. The gradual slope may result from the emergence of a temporal extremely large cluster in fluctuation as the system size increases. These confirm that the characteristic length scale is independent of the system size N, while determined by its self-propulsion parameters, consistent with the predicted short wavelength instability.
Finite size analysis for the cluster size distribution. (a) n(s) at Γ = 0.24 for different number of particles N in the system indicated by colors. (b) n(s) at Γ = 2.4. (c) Maximum cluster size Smax as a function of system size N, and the colors indicate Γ. The red dashed line indicates ∝N, and the black dashed line represents a constant. Here, Dr = 0.1, Per = 240, and ϕ = 0.40.
Finite size analysis for the cluster size distribution. (a) n(s) at Γ = 0.24 for different number of particles N in the system indicated by colors. (b) n(s) at Γ = 2.4. (c) Maximum cluster size Smax as a function of system size N, and the colors indicate Γ. The red dashed line indicates ∝N, and the black dashed line represents a constant. Here, Dr = 0.1, Per = 240, and ϕ = 0.40.
Circulating current caused disintegration
As the conventional MIPS cohesion effect is preserved in the linear order hydrodynamics of chiral swimmers, to form dynamic clusters, there should be a disintegration mechanism. Thus, we explore the steady state current, which vanishes in equilibrium systems as a result of the balance of chemical potential and pressure everywhere.42 For slowly varying fields, we neglect the temporal derivative and the viscosity term in Eq. (3) to obtain the adiabatic solution of the polarization field,
where j(ρ) = ρve(ρ) represents the scalar strength of the effective current. As Jad = veWad, we find a non-curl-free current in steady states of chiral swimmers, . The circulating current originates from the mixture of spatial gradients induced by automatic rotation. Intuitively, such system spanning circulating current can continuously disassemble the dynamic clusters. Different from the equilibrium-like MIPS in ABPs, the steady state current in cABPs is a distinct feature of the broken detailed balance.
To understand the microscopic origin of the circulating current, we measure the self-intermediate scattering function (sISF), which reveals dynamic patterns at different spatiotemporal regimes,
Here, ⟨·⟩ represents the ensemble average, and we define nq = qL/(2π). The sISF of an individual circle swimmer was studied in Refs. 43 and 44. First, in the large time and length scale, a diffusion-like behavior is expected, as the persistent swimming is smeared out by frequent particle collisions (in the parameter region used in simulation, τb/τr ≃ 4 × 10−3). Thus, the correlation follows a near-exponential decay at t ≫ τb and small nq [the dashed curve in Fig. 5(a)]. While in MIPS at slow rotation Γ ≲ 1, Fs(q, t) relaxes in a much slower manner than the diffusive motion, reflecting the particles frozen inside the stable dense droplet, confirming that the particle exchange is an interfacial process in bulk phase separation. Second, to observe the circular motion, one should look at a large enough wavelength Λ = L/nq ≳ R, and thus, we choose nq = 1. Intriguingly, for fast rotation Γ ≳ 1, Fs(q, t) oscillations are around a finite plateau at the time scale around τw, suggesting that the chiral swimming pattern of an individual particle is partially preserved in the dynamical clustering state, even when dense liquid droplets are formed. Finally, at a large wave number nq = 5, the oscillations of Fs(q, t) around zero within τr indicate a persistent motion |Δr(t)| = vefft at a short wavelength; until at longer times, the rotational diffusion washes out the oscillations. The increasing oscillation amplitude with Γ indicates a larger veff, thus suggests less velocity damping effect at faster rotation. And the increase of veff finally saturates when Γ reaches the homogeneous state. We notice that the instantaneous effective swimming velocity ve is independent of ω0, while the time averaged persistence velocity veff increases with ω0, as shown in the sISF. This inconsistency suggests that fast enough rotation can “relieve” trapped particles from temporal fluctuations before being smeared out by the rotational diffusion, thus preserving the chiral swimming pattern and giving rise to the system spanning circulating current.
Self intermediate scattering function of circle active Brownian particles. (a) nq = 1, the dots indicate the time scale τω/τr. (b) nq = 5, the color indicates Γ, and the inside bracket is τω/τr. Here, Dr = 0.1, ϕ = 0.40, and Per = 360.
Self intermediate scattering function of circle active Brownian particles. (a) nq = 1, the dots indicate the time scale τω/τr. (b) nq = 5, the color indicates Γ, and the inside bracket is τω/τr. Here, Dr = 0.1, ϕ = 0.40, and Per = 360.
CONCLUSION AND DISCUSSION
In conclusion, a dynamical clustering state emerges in cABPs with sufficiently fast rotation and interrupts the conventional MIPS. The underlying short wavelength instability originates from the combination of the conventional MIPS cohesion, and the circulating current caused disintegration. Different from the most linear active particles, a detailed balance cannot be restored at the continuum level for chiral active particles, which allows for the nonzero curl current in steady states. This is surprising, as there is no explicit alignment interaction between chiral active colloids to induce any collective spontaneous flow in the system.25 More generally, the instability introduced by simply switching on monochromatic rotation might help to design active colloids that can achieve complex spatiotemporal pattern formation. Finally, our key finding of the well-defined non-equilibrium steady state may also provide a platform to examine the generalization of effective thermodynamic concepts9,45 previously developed in the context of MIPS.
SUPPLEMENTARY MATERIAL
See the supplementary material for the derivation of the continuum theory and the linear stability analysis and a video for the dynamical clustering phase.
ACKNOWLEDGMENTS
This work was supported by the Singapore Ministry of Education through the Academic Research Fund under Grant No. MOE2019-T2-2-010. We thank NSCC for granting computational resources.
AUTHOR DECLARATIONS
Conflict of Interest
The authors declare no conflict of interest.
DATA AVAILABILITY
The data that support the findings of this study are available within the article and its supplementary material and from the corresponding author upon reasonable request.