Nonequilibrium strongly hyperuniform fluids of circle active particles with large local density fluctuations

See allHide authors and affiliations

Science Advances  25 Jan 2019:
Vol. 5, no. 1, eaau7423
DOI: 10.1126/sciadv.aau7423


Disordered hyperuniform structures are an exotic state of matter having vanishing long-wavelength density fluctuations similar to perfect crystals but without long-range order. Although its importance in materials science has been brought to the fore in past decades, the rational design of experimentally realizable disordered strongly hyperuniform microstructures remains challenging. Here we find a new type of nonequilibrium fluid with strong hyperuniformity in two-dimensional systems of chiral active particles, where particles perform independent circular motions of the radius R with the same handedness. This new hyperuniform fluid features a special length scale, i.e., the diameter of the circular trajectory of particles, below which large density fluctuations are observed. By developing a dynamic mean-field theory, we show that the large local density fluctuations can be explained as a motility-induced microphase separation, while the Fickian diffusion at large length scales and local center-of-mass-conserved noises are responsible for the global hyperuniformity.


Perfectly ordered structures, such as crystals or quasi-crystals at zero temperature, are usually associated with some discrete symmetries and exhibit long-range correlations, leading to the structure factor of the system S(q → 0) = 0 (1). Similarly, the local density variance in these structures 〈δρ2〉 scales with the window size of observation L as 〈δρ2〉 ~ L−λ with λ = d + 1, where d is the dimensionality of the system. In contrast, in normal disordered structures, e.g., conventional gases, liquids and amorphous solids, and even thermalized crytals, the long-wavelength density fluctuation makes S(q → 0) = const. > 0 and λ = d. Recently, the concept of hyperuniformity was introduced to study the state of matter (2). A structure is hyperuniform when it has vanishing long-wavelength density fluctuations, i.e., S(q → 0) ~ qα → 0 with α > 0 and 〈δρ2〉 ~ L−λ with d < λ ≤ d + 1 (2). It has been found in the past two decades that, besides the ordered hyperuniform structures, i.e., perfect crystals and quasi-crystals, a number of disordered structures are also hyperuniform, including the maximally random jammed packing (3), avian photoreceptor patterns (4), and some nonequilibrium systems (510).

Disordered hyperuniform structures have received an increasing amount of scientific attention, as some strongly hyperuniform disordered structures with λ = d + 1 exhibit similar properties as crystals with even better performance, e.g., large isotropic photonic bandgaps insensitive to defects (11) can be opened at low dielectric contrast (12). Although ideal hyperuniform structures similar to perfect crystals are unavoidably affected by thermal excitation, it still shows promise in designing robust disordered materials with novel functionalities (1, 1316). By far, various protocols were developed to design particle interactions to form disordered hyperuniform ground states in classical many-particle systems. However, the resulting interactions normally have delicate long-range or multibody terms (1), making their experimental realization highly challenging. An alternative approach is to use driven systems, e.g., self-organized colloidal suspensions under periodic shearing (17), to form nonequilibrium dynamic hyperuniform states, which can effectively avoid the dynamic trapping, and, in principle, have a higher resistance to thermal perturbations (7). Now, experimentalists only succeeded in generating weakly hyperuniform structures with λ ≃ d + 0.5 (8) using this method. Nevertheless, this is certainly a direction that deserves further investigation (10), as some self-driven systems, or active matter systems, have produced a number of unexpected emergent phenomena never found in corresponding equilibrium systems (1822). However, in conventional active matter systems, i.e., active nematic systems and active Brownian particles systems, giant number fluctuations characterized by λ < d (19) and motility-induced phase separation (MIPS) with λ ≃ 0 (2023) are usually observed. These large density fluctuations seemingly prohibit the formation of hyperuniform structures in active matter systems.

Recently, chiral active matter whose motion is chiral-symmetry broken, e.g., active particles/swimmers with circular motion in two-dimensional (2D) (2428) or active spinner fluids (29), has attracted considerable attention. Both experiments and simulations have shown many interesting collective phenomena in these systems (25, 3033). In this work, using computer simulations combined with analytic theories, we study a 2D system of circle active particles, which perform independent circular motion with the same handedness and random circling phases. We show that in the limit of strong driving or zero thermal noise, with increasing density of particles or radius of circular motion R, the system undergoes an absorbing-active transition forming a nonequilibrium strongly hyperuniform fluid phase with density variance 〈δρ2〉 ~ L−3 (L → ∞) the same as in perfect crystals. Further increasing the density or R triggers the formation of dynamic clusters, which results in large local density fluctuations. These fluctuations are “confined” within the length scale of R, while the strong hyperuniformity persists at large length scales. This surprising coexistence of large local density fluctuations and the global hyperuniformity is explained by dynamic mean-field theories at different length scales.



As illustrated in Fig. 1A, we consider a 2D suspension of N active colloidal particles with diameter σ. Each particle experiences an in-plane force Fp with random initial orientation and a constant torque Ω perpendicular to the plane, which drive the particles to perform circular motion with the same handiness (30, 31). The dynamics of particle i at finite temperature T is governed by the over-damped Langevin equations (26, 30)r.i(t)=γt1[iU(t)+Fpei(t)]+2kBT/γtξit(t)(1)e.i(t)=[γr1Ω+2kBT/γrξir(t)]×ei(t)(2)where ri and ei are the position of particle i and its self-propulsion orientation, respectively, and kB is the Boltzmann constant. γt/r is the translational/rotational friction coefficient. For simplicity, we set γt = γr2. ξit(t) and ξir(t) represent Gaussian noises with zero mean and unit variance. We use Weeks-Chandler-Andersen (WCA) potential to mimic the excluded volume interaction between colloidal particles U(t) (see Methods). The packing fraction of the system is defined as ϕ = ρσ2π/4, where ρ is the particle density. The self-propulsion speed of particle is v0=γt1Fp. The reduced noise strength in the system is defined as TR = kBT/(Fpσ) that measures the strength of thermal noise compared with the self-propulsion. In the zero noise limit, i.e., TR = 0, isolated active particles perform circular motions with fixed radius R = Fpσ2/Ω and period Γ = 2πγr/Ω. This athermal noise-free situation is the major focus of this work, and the effect of thermal noise is discussed later.

Fig. 1 Absorbing-active transition.

(A) Schematic of the 2D system of circle active particles. (B) Dynamic phase diagram in the representation of packing fraction ϕ and circle radius R. (C and D) Typical snapshots of the active and adsorbing states near the critical point with ϕ = 0.20 with R = 1.75σ, where the color indicates the self-propulsion orientation of each particle. These two states are marked as magenta and cyan solid symbols, respectively, in (E). (E) MSD as functions of time for system with R = 1.75σ started from random configurations. Red line, active state (ϕ = 0.22); blue line, absorbing state (ϕ = 0.01); green line, system near the critical point (ϕ = 0.20) in which the system ultimately falls into the absorbing state after a long simulation time. (F and G) Diffusion constant as functions of ϕ near the critical point ϕc for systems with different R. The dashed lines are the fitting of power law (ϕ − ϕc)β. (H and I) The structure factor S(q) and the orientation correlation function C(r) of active (magenta) and absorbing (cyan) states as marked by solid symbols in (E). (J) The measured critical exponent β from (G) as a function of R. For all the calculations, N = 10,000 and TR =0.

Absorbing-active transition

We first simulate systems with N = 10,000 and TR = 0. Under low packing fraction ϕ and small R condition, we find that the system falls into an absorbing or arrested state, in which each particle performs an independent circular motion without collisions and the mean squared displacement (MSD) 〈Δr2〉 of particles develops a plateau at long time (blue line in Fig. 1E). With increasing ϕ or R, the collisions between particles become more frequent, making the system unable to find a noninteracting state. Thus, the system remains at an active diffusive state with MSD 〈Δr2〉 ~ 4Dt (t → ∞) (red line in Fig. 1E). Here, D is the long-time diffusion constant. The phase behaviors of the system are summarized in the phase diagram Fig. 1B. In the following, we focus on the absorbing-active transition close to the boundary of two phases.

In Fig. 1F, we plot D as a function of ϕ − ϕc for different R. Here, we obtain ϕc by fitting with the critical power law D ~ (ϕ − ϕc)β, which determines the phase boundary in Fig. 1B. One can observe a sharp transition from the absorbing state [D(ϕ) = 0] to the active state [D(ϕ) > 0] when increasing ϕ at a small R = 1.75σ. The transition becomes smoother as R increases. In Fig. 1G, we show the log-log plot of D as a function of ϕ − ϕc. The obtained slope β is given in Fig. 1J. We find that the critical exponent β is about 0.15 for R = 1.75σ, which is substantially smaller than the values in classical absorbing transitions, i.e., β = 0.58 for directed percolation and β = 0.64 for conserved directed percolation (Manna type) (634). Such a small critical exponent is independent of system size (see fig. S1). With increasing R, we find that β increases to around 0.5 for R > 10σ. A similar increase of the critical exponent with increasing interaction range (in our case, R) has been reported (35). In our system, ϕc would decrease to zero with increasing R. Hence, β at large R cannot be directly obtained in our system due to the divergence of simulation time needed at the dilute limit.

To understand the physics behind the absorbing-active transition at small R, we choose a packing fraction ϕ = 0.20 close to the critical point (ϕc = 0.195) for system with R = 1.75σ. The MSD for the system started from random configuration is shown by the green line in Fig. 1E, in which one can see that the system ultimately falls into the absorbing state after staying at the active state for a long time. In Fig. 1 (C and D), we show typical snapshots for the active state and absorbing state before and after the absorbing transition, as indicated by the magenta and cyan symbols in Fig. 1E. Movies for these two states can be found in movies S1 and S2. One can notice a marked structural difference between these two states. In the absorbing state, particles with similar orientation form finite synchronized clusters, while the active state is more homogeneous without noticeable spatial heterogeneity. This structural difference is also reflected in the structure factor S(q) and orientation correlation function C(r) = 〈∑ij eiej δ(rijr)〉/ρN, as shown in Fig. 1 (H and I). Compared with the active state, S(q) for the absorbing state develops a pronounced peak at qσ ≃ 0.2, and the corresponding C(r) also shows a stronger orientation correlation. The synchronized clusters formation in our circle active particle system with isotropic circling-phase distribution shares a similar mechanism with the phase separation observed in an experimental bimodal phase-distributed system (see also fig. S4) (26). Both are a result of a crowding-induced attraction between particles with the same circling phase. We find that this distinct structural transformation during the absorbing transition is absent in the conventional absorbing transition (5, 6, 17, 34), suggesting that the small critical exponent measured in our system has a structural origin. With increasing R, the structural difference between two phases becomes weaker (see fig. S2), which occurs simultaneously with the increase of β. Further studies combining with finite-size analysis are necessary for determining whether the absorbing-active transition at small R is first order.

Hyperuniformity and large local density fluctuations

As shown in Fig. 1H, for the active state near the critical point in the system with R = 1.75σ, the structure factor of the system exhibits hyperuniform scaling S(q) ~ q(q → 0). In the random organization model aiming to mimic the colloidal suspension under periodic shearing, a similar hyperuniform scaling was observed near the critical point (5, 6). However, for relative large R, the critical q scaling shifts to a faster q2 scaling (see fig. S2). To explore this further, we simulate a system of N = 40,000 circle active particles. In Fig. 2A, we first plot the MSD for active state systems with different R at ϕ = 0.2. We find that the diffusivity in the system rises slightly with increasing R. By further checking the scalings of the density variance 〈δρ2〉 and S(q) for different R (see Fig. 2, B and C), we observe a strong hyperuniformity in the systems with R ≤ 25σ, as indicated by the asymptotic behaviors: 〈δρ2〉 ~ L−3 (L → ∞) and S(q) ~ q2(q → 0). From 〈δρ2〉, one can identify an R-dependent length scale LHU, above which the system becomes hyperuniform, while below which the system behaves like normal fluids, i.e., 〈δρ2〉 ~ L−2 (Fig. 2B). In S(q), a similar threshold qHU can be found, which features the end of the hyperuniform scaling. With increasing R, LHU increases and qHU decreases until R ≥ 50σ, where the hyperuniform scaling becomes less apparent due to the system finite-size effect as discussed below. This implies that R controls the length scales at which the system exhibits hyperuniformity.

Fig. 2 Dynamic hyperuniform state.

(A and D) MSD as functions of t for various R. (B and E) Density variances 〈δρ2〉 as functions of window size L for various R. The L−3 asymptotic line indicates the hyperuniform scaling, which is the same as in perfect crystals. The L−2 scaling is for normal fluids, while L0 is for clustering- or phase separation–induced large density fluctuations. (C and F) Structure factor S(q) for various R. The q2 asymptotic line indicates the hyperuniform scaling, while the q−2 line represents clustering- or phase separation–induced large density fluctuations. (G) Typical snapshots for systems at ϕ = 0.4 with various R. For all the calculations, N = 40,000 and TR.

Figure 2 (D to F) shows the result of an analogous investigation for higher-density systems with ϕ = 0.4. Compared with low-density systems, we observe pronounced enhancement of the diffusivity with increasing R (Fig. 2D). The high-density systems also show clear hyperuniform scaling at large length scales for R ≤ 50σ and the threshold LHU (or qHU) increases (or decreases) with increasing R (Fig. 2, E and F). However, at ϕ = 0.4, we find 〈δρ2〉 ~ L−λ with λ < 2 when LLHU and λ decreases with larger R. This implies that at length scales LLHU, the system exhibits large density fluctuations, whose strength and length scale are both controlled by R. This large fluctuation can also be identified by the scaling S(q) ~ q−2 for q > qHU as shown in Fig. 2F, which was reported as a signature of critical instability of active particle system undergoing MIPS (20). However, in most of our systems except R ≥ 100σ, S(q) ~ q−2 does not diverge at qHU ≃ 0 as in MIPS (20) but stops at a finite qHU. This suggests that the large local density fluctuations observed in our system are due to clustering or microphase formation. The crossover of two different scalings, i.e., large density fluctuations and hyperuniformity, creates a peak in S(q) at qHU, whose height increases with larger R. Actually, for R ≥ 100σ, we speculate that in much larger systems, one can still observe the peak at finite qHU for ϕ = 0.4 and the hyperuniform scaling for both ϕ = 0.2 and 0.4, as suggested by later theoretical analyses. In fig. S3, S(q) is shown for a larger system (N = 102,400) at ϕ = 0.4. The result agrees with our speculation. Typical snapshots of the system at ϕ = 0.4 with various R are shown in Fig. 2G (also in movies S3 to S7 and fig. S3), and one can see many finite-size clusters disappearing and reforming in the system. The average size of these dynamic clusters increases with increasing R, and at R = 1000σ, because of the finite-size effect, the clusters percolate the simulation box. These findings are intriguing, as hyperuniformity and large density fluctuations induced by dynamic cluster formation are two seemingly opposite phenomena, which coexist here in the same system at different length scales. In the following, we formulate dynamic mean-field theories to understand this new dynamic hyperuniform fluid with large local density fluctuations.

Dynamic mean-field theory

Starting from the N-body Smoluchowski equation for active Brownian particles (20, 21), one can prove (section S1) that in a homogeneous circle active particles system with vanishing orientation order parameter Q=eeT121, the time-dependent local density field ρ(r, t) and the local polarization field p(r, t) satisfytρ=[ve(ρ)pDeρ](3)tp=12[ve(ρ)ρ]+De2p+Ωr×p(4)where ve(ρ) = v0 + ζρ is a density-dependent effective velocity of particles with a negative ζ reflecting the motility-induced “self-trapping” effect. De is the effective diffusion constant originated from the “evasive” motion of particles due to the collisions with neighboring particles (21). Ωr=γr1Ω is the reduced torque. The isotropic homogeneous state, i.e., [ρ(r,t)=ρ¯,p(r,t)=0], is the solution to Eqs. 3 and 4. By making a weak perturbation around this state, i.e., [ρ(r,t)=ρ¯+δρ(r,t),p(r,t)=δp(r,t)], we obtain two linearized equations in the Fourier space with the first-order approximation(iω+Deq2)δρq,ω=iveqpq,ω(5)(iω+Deq2)pq,ω=Ωr×pq,ωiwqδρq,ω(6)where [δρq, pq,ω] = ∫dreiqrdt eiωt[δρ, p], and w=v0/2+ζρ¯ is the parameter indicating the strength of self-trapping effect. By solving Eqs. 5 and 6, one obtains the dispersion relationship, which includes a diffusive mode ω0 = iDeq2 and two nondiffusive modesω1,2=iDeq2±vewq2+Ωr2(7)

The growth rate of the mode is κ = Re(iω), whose sign determines whether the perturbation δρ ~ eiωt+iqr grows or decays. One can prove that the mode 1 always decays, while the mode 2 may grow for w < 0 withκ2={Deq2q<v0vewR1Deq2+vewq2Ωr2q>v0vewR1(8)

In Fig. 3A, we plot three typical situations of κ2 as functions of q at high density (w < 0) by varying R (or Ωr). One can see that for finite R, κ2 decreases from zero following the typical diffusive mode −Deq2 from q = 0 to v0vewR1, above which κ2 starts to increase and has the chance to be positive at finite q > 0. This implies that the homogeneous system becomes unstable as a result of the growing fluctuation of finite wavelength. In the mean-field picture, this typically indicates that the system undergoes a microphase separation. This scenario is different from the complete phase separation in which instability starts from the infinite wavelength at q = 0 (20, 21).

Fig. 3 Dynamic microphase separation.

(A) Growing rate κ2 as functions of q for systems at ϕ = 0.35 with various R obtained from the dynamic mean-field theory. (B) Measured heights of the first peak in S(q) for systems with different combinations of R and ϕ in computer simulations indicated by the color of the symbols. The dotted line is the fitting using Eq. 9 for the phase boundary. Inset shows the measured position of the first peak in S(q) in computer simulations (symbols) and the theoretical prediction (dotted line) based on the fitting in (B).

The instability point of the system is defined as κ2max(q*)=0 for q* > 0, from which we can obtain the relationship between the critical packing fraction ϕ* and q* (see section S2)(ϕ*ϕc1)(2ϕ*ϕc)=8Dev0R(9)q*=v0RDe(10)where ϕc ≃ 0.32 is the critical packing fraction for R = ∞ at which the simulation shows complete MIPS (q* = 0). This suggests that MIPS can be seen as a limit of the microphase separation in our system when R → ∞ at the mean-field level. For this special case, when ϕ approaches ϕc, S(0) is supposed to jump from a finite value to ∞, which marks the “spinodal” of MIPS. For systems with finite R, the first peak of S(q), which is located around q*, is also expected to jump up to a higher value when the system crosses the “microphase separation point.” In Fig. 3B, we plot the measured heights of the first peak in S(q) for various combinations of R and ϕ in computer simulations. We find that when R > 50σ, the height of the first peak in S(q) jumps from below 5 to above 10 with increasing ϕ over ϕ*. The transition becomes smoother for systems of smaller R. Therefore, we choose S(q) = 7 as the threshold to fit the phase boundary using Eq. 9 (the dotted line in Fig. 3B), which leads to De ≃ 0.62v0σ. In the inset of Fig. 3B, we plot the location of q* measured in simulations and the prediction of Eq. 10 with the same De. The simulation results agree almost quantitatively with the theoretical prediction. For small R cases (R < 50σ), there is a noticeable inconsistency, suggesting the possible breakdown of the mean-field theory. In these cases, the size of clusters in the system is comparable with the particle size. Therefore, microphase separation is no longer a proper concept to describe what happens in the system. Even for large R cases, which seem consistent with the mean-field microphase separation scenario, it still remains open whether there exist true sharp microphase separation points for finite R or if the sharpness of the transition diverges at R → ∞.

To summarize, this theoretical analysis rationalizes the instability in the circle active particle system as a result of motility-induced self-trapping effect (negative w in Eqs. 7 and 8). The existence of nonzero torque Ωr restricts the growth of density fluctuations within a finite wavelength proportional to the circle radius R, which is the underlying reason for the dynamic clusters formation and the resulting large local density fluctuations observed in our simulations.

Mechanism of global hyperuniformity

To understand the hyperuniformity at large length scales, we focus on the spatial distribution of the instantaneous circling centers of active particles, which we treat as effective particles with radius R. For particle i at position ri, the circling center is at rio=riσ2(Fip×Ω)/|Ω|2, which is a fixed position for an isolated circle active particle. The motion of these effective particles is due to collisions with other particles. This is similar to what happens in the random organization model (5, 6, 10), where only overlapped particles experience random kicks. As indicated in the previous theoretical analysis, the self-trapping and the growing density fluctuations are confined within the length scale of R. At larger length scales, fluctuations decay as the diffusive mode (Fig. 3A). These imply that the dynamic equation for the density field of these effective particles ρo(r, t) at large length scales can be described by the Fick’s law of diffusion ∂tρ = D2ρ in the q spacetρqo=Deoq2ρqo+ξq(t) for (q2πR)(11)with [ρqo,ξq]=dreiqr[ρo,ξ]. Here, Deo is the diffusion coefficient of effective particles and ξq(t) is an additional noise term. Since we neglect the thermal noise, ξq(t) comes from the chaotic multiparticle interaction (36), which obeys the center of mass conservation (CMC). As proven in (10), the noise with additional CMC appears as a double-space derivative in the diffusion equation, i.e., ξ(t)=ρ¯2η(t), with 〈η(r, t)η(r′, t′)〉 = A2δ(rr′)δ(tt′) and A as the strength of the noise. This is different from the conventional single-space derivative noise that arises solely from the particle number conservation (18, 37). Following (10), by adding such a double-space derivative noise term into Eq. 11 as the perturbation, one obtains the dynamic equation for the density fluctuation δρo in the Fourier spaceiωδρq,ωo=Deoq2δρq,ωoq2ρ¯ηq,ω for (q2πR)(12)where ηq is the Fourier transform of the noise η(r, t). From Eq. 12, we obtain the structure factorSo(q)=12πτmaxNδρq,ωoδρq,ωo*dω=A22Deq2 for (q2πR)(13)where τmax is the maximum observation time (see section S3). Equation 13 shows the same hyperuniform exponent in S(q) as observed in the simulations and also shows the strongest hyperuniform scaling for density variance 〈δρo2〉 ~ L−3 at the length scales larger than R (2). In Fig. 4A, we plot the scaled density variance 〈δρo2R2 versus L/R for systems with various R at ϕ = 0.2, and one can see all points collapse into a single curve. The curve consists of two distinct scalings with a strongly hyperuniform scaling L−3 at L > LHU ≃ 2R and a normal fluid-like scaling L−2 at L < LHU. Similar collapse is shown in Fig. 4B for So(q) versus qR, where the hyperuniform scaling So(q) ~ q2 stops at qHU ≃ 2π/LHU. For high-density systems, i.e., ϕ = 0.4, we obtain the same crossover length scale LHU ≃ 2R as shown in Fig. 4 (C and D). In this case, the density variance and So(q) for different R do not collapse into a single curve because of cluster formation–induced large density fluctuations at LLHU. After obtaining the hyperuniform scaling for these effective particles, one can prove the existence of the same hyperuniform scaling at similar length scales for the circle active particles (38).

Fig. 4 Global hyperuniformity.

(A and C) Scaled density variance 〈δρo2R2 as functions of L/R and (B and D) So(q) as functions of qR, for various R at ϕ = 0.2 (A and B) and 0.4 (C and D), respectively.

From these analyses, one can see that, in systems of circle active particles, there are two ingredients for the global hyperuniformity: (i) The circular trajectories of active particles localize their active motions, creating the Fickian diffusion condition for the active particles at large length scales; and (ii) the local density fluctuations originate from the chaotic multiparticle interactions, which obey the CMC. These arguments are based on (10), which produces the same hyperuniform scaling in S(q) using the random-organization model with CMC. Although the two systems share the similar mechanism of hyperuniformity at large length scales, there are also marked differences between them. First, the driving force in our system is persistent, and the dynamics of the active particles is overdamped and deterministic. Second, our system has a characteristic length R, which separates two completely different fluctuations. These cannot be described by the model in (10).

We notice that the two ingredients for the dynamic hyperuniformity above seem to be satisfied in a recent experimental system, in which spherical Janus colloids with additional soft repulsion perform almost perfect circling motion with the same chirality but two different phases, driven by external fields (26). The thermal noise was believed to be negligible in this system, and the interactions between these nearly athermal active particles produce an effective temperature, which controls both kinetic and phase behavior of the system (26). To test whether the hyperuniformity can exist in this circle active particle system with the bimodal circling-phase distribution, we simulate a low density system (ϕ = 0.2) with two different circle radii R = 1σ and 3σ using the model (see Methods), which produced almost the same experimental result in (26). Our result is shown in fig. S4, from which we find the bimodal-distributed system falls into a phase-separated absorbing state for R = 1σ but stays in an active mixing (lane) state for R = 3σ, consistent with the previous finding (26). In the active mixing state, we observe the same hyperuniform scaling S(q → 0) ~ q2. This result demonstrates the robustness of the hyperuniformity mechanism unveiled by this general model and suggests a high possibility of realization in experiments. Moreover, it also suggests that other active systems satisfying the two ingredients, such as active spinner systems (29) or size/interaction oscillating particle systems (39), may be used to produce the same hyperuniformity as well.

Effect of thermal noise

According to the definition, a perfectly hyperuniform state requires S(0) to exactly equal zero. However, in real experimental systems, thermal noise is unavoidable. On the basis of the fluctuation-compressibility relationshipS(0)=κTρkBT(14)any thermal equilibrated systems with the positive isothermal compressibility κT at finite temperatures cannot be strictly hyperuniform due to thermal excitation (13). In crystals, thermal excitation appears as phonon modes, which cause background scattering or thermal diffuse scattering, whose effect can be measured by the Debye-Waller factor (13, 40). Nevertheless, the wide application of crystal materials indicates that thermalization only weakens but does not destroy most physical properties of the ground-state crystal. Similarly, near-hyperuniformity (1, 13, 15) in disordered structures may also be enough to achieve some desired functions, e.g., isotropic photonic/electronic bandgaps (14, 15). In Fig. 5, we analyze the influence of thermal excitation on the ground-state hyperuniform system with ϕ = 0.2 and R = 3σ. We find that with a gradual increase of the reduced noise strength TR from zero, S(q → 0) begins to saturate at some nonzero value, which increases along with TR (Fig. 5A). In section S4, we introduce the thermal noise f = [fx, fy] in Eq. 11 asξ(t)=ρ¯[η(t)+f(t)](15)where 〈fi(r, t)fj(r′, t′)〉 = 2Dthermδijδ(rr′)δ(tt′) with Dtherm=kBTγt1(1+R2/σ2) the self-diffusion constant of effective particles due to the thermal Brownian motion (25). Here, we assume that the thermal noise is a first-order weak perturbation on the chaotic noise η(t), which leads to the decoupling of these two noise sources: 〈fi(r, t)η(r′, t′)〉 = 0. Then, we can estimate So(q) as a function of the reduced noise strength TR for low-density systems asSo(q)=v0σDe(1+R2σ2)TR+A22Deq2 for (q<qHU)(16)

Fig. 5 Effect of thermal noises on hyperuniformity.

(A) Structure factor So(q) under different reduced noise strength TR for ϕ = 0.2 and R = 3σ. Open symbols show the simulation data, while the dashed lines are the theoretical predictions of Eq. 16. (B) Normalized So(0) as functions of TR from theoretical prediction (dashed line) and the fitting of simulation results (solid symbols) for systems with ϕ = 0.2. For all the calculations, n = 40,000.

In Fig. 5A, we plot the theoretical prediction of Eq. 16 as dashed lines for different TR values at R = 3σ and ϕ = 0.2, by assuming them all across the same So(qHU) point. In Fig. 5B, we compare the So(0) from the theoretical prediction with the So(0) obtained from the fitting of simulation results (see Methods) for systems with R = 3σ and 10σ at ϕ = 0.2. One can find quantitative agreements between simulation and theoretical predictions in Fig. 5, which suggest that the susceptibility of hyperuniformity to the noise in our nonequilibrium system is similar to that in thermalized crystals, i.e., Eq. 14. However, we also emphasize the difference: In our nonequilibrium system, the saturated value of So(0) is determined by the driving force as well. Therefore, in experiments, it is possible to observe a large range of hyperuniform scaling in S(q) or density variance as long as the driving force is much larger than the thermal noise, i.e., TR ≪ 1.


In conclusion, by combining computer simulations with theoretical analyses, we investigate the dynamic phase behaviors in 2D systems of circle active particles. In the zero-noise limit, we find that with increasing the density of system or the radius of circular motion R, the system undergoes a transition from an absorbing state to an active fluid state, which is accompanied by a structural transformation for small R. In the active fluid state, we find a characteristic length scale LHU. For LLHU, the system exhibits strong hyperuniformity with the density variance scaling the same as in perfect crystals, while for LLHU, we observe normal random fluctuations at low density and large density fluctuations (cluster formation) at high density. To understand the mechanism of the phase behaviors of the system, we develop a dynamic mean-field theory. Linear stability analysis suggests that at the mean-field level, the large local density fluctuations at relative large R are a result of motility-induced microphase separation, which is confined within the length scale of R. For the global hyperuniformity, we attribute it to the interplay between the Fickian diffusion of active particles at large length scales and local particle collisions that conserve the center of mass (10). Our work demonstrates that two extreme fluctuations, i.e., large density fluctuations and hyperuniformity, can coexist in the same dynamic system at different length scales. We emphasize that this stable hierarchical hyperuniform fluid is conceptually different from the disordered hyperuniform solid or critical hyperuniform state. From a practical point of view, our results suggest that even for exotic disordered hyperuniform structures, there is still plenty of room at the “bottom,” i.e., one may construct arbitrary local complex structures (ordered or disordered) with extra functionalities without harming the global hyperuniformity. This provides large freedom in designing hierarchical disordered hyperuniform materials with unconventional properties.


In our simulations, we used a square simulation box with periodic boundary conditions in all directions, starting from initial configurations with random particle positions and orientations to make sure the initial structure factor S(q) ~ 1. The time unit was chosen to be the time that a particle moves a distance of σ in the dilute limit, i.e., τ0 = σ/v0. To mimic the excluded volume interaction between colloidal particles i and j, we used the WCA potentialU(rij)={4ε[(σrij)12(σrij)6+14](rij<21/6σ)0(rij>21/6σ)(17)where rij is the center to center distance between particles i and j, with σ as the diameter of particles. We chose ε = Fpσ/24, which gives the typical contact distance rc = σ between particles based on force balance Fp=U(rij)rij|rij=rc. For systems at TR = 0, to exclude the noise generated by discrete dynamic integrations, we used perfect convex polygons to approximate the closed circle trajectories of active particles. This is realized by finely tuning the integration step, which is around 10−3τ0.

For systems with bimodal distributed circling-phase, following (26), we added an additional soft repulsion Us = A (r/σ)−4 with cutoff distance rsc = 5σ and A = 7.5ε to model the isotropic dipole interaction between active particles. The self-propulsion force for this system was reset to Fp = 27/6A/σ to make the dipolar interaction balance the driving force at rc = 21/6σ (26).

The density variance 〈δρ2〉 of the system was calculated using a spherical window whose diameter is smaller than the half of simulation box to avoid the finite-size effect. Under the periodic boundary condition, the structure factor S(q) was calculated for some discrete q vectors, i.e., [qx,qy]=2πL0[i,j](i,j=1,2,3), with L0 as the size of cubic simulation box. The simulation time for the equilibrating and sampling processes depends on the system size and density. For hyperuniform system, the criterion for the system to reach equilibrium is whether the hyperuniform scaling of S(q) has fully extended to the smallest qmin = 2π/L0 without further change. The fitting function used in Fig. 5B is Eq. 16, with an adjustable TR.

Correction (3 April 2020): In an earlier version, Fig. 1F contained an incorrectly plotted orange dashed line. The figure has been updated with the right set of data and correct plots.


Supplementary material for this article is available at

Section S1. Derivation of the dynamic mean-field theory for 2D system of circle active particles.

Section S2. Linear stability analysis.

Section S3. Calculation of So(q) for system with CMC.

Section S4. Effect of thermal noise on So(q).

Fig. S1. Finite-size effect on the long-time diffusion coefficient D as a function of packing fraction ϕ for systems with R = 1.75σ.

Fig. S2. Structural comparison of active and absorbing state at R = 10σ near the critical point ϕc = 0.0194.

Fig. S3. Structure factor S(q) for large systems with different R at φ = 0.4 and TR = 0.

Fig. S4. Hyperunifomity in an experimentally realizable system (26) with bimodal circling-phase distribution.

Movie S1. Active state in Fig. 1C.

Movie S2. Absorbing state in Fig. 1D.

Movie S3. Active state with R = 1000σ in Fig. 2G.

Movie S4. Active state with R = 100σ in Fig. 2G.

Movie S5. Active state with R = 50σ in Fig. 2G.

Movie S6. Active state with R = 25σ in Fig. 2G.

Movie S7. Active state with R = 23σ in Fig. 2G.

This is an open-access article distributed under the terms of the Creative Commons Attribution license, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


Acknowledgments: We thank H. Zhang in Shanghai Jiao Tong University for fruitful discussions. We are grateful to the National Supercomputing Centre (NSCC) of Singapore for supporting the numerical calculations. Funding: This work is supported by Nanyang Technological University Start-Up Grant (NTU-SUG; M4081781.120), the Academic Research Fund from Singapore Ministry of Education [M4011616.120, M4011873.120, and MOE2017-T2-1-066 (S)], and the Advanced Manufacturing and Engineering Young Individual Research Grant (A1784C0018) by the Science and Engineering Research Council of Agency for Science, Technology and Research Singapore. Author contributions: Q.-L.L. and R.N. conceived the research. Q.-L.L. performed the research. R.N. directed the research. Q.-L.L., M.P.C., and R.N. analyzed the data and wrote the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are presented in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article