INTRODUCTION
============

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 ([@R1]). 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 ([@R2]). 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 ([@R2]). 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 ([@R3]), avian photoreceptor patterns ([@R4]), and some nonequilibrium systems ([@R5]--[@R10]).

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 ([@R11]) can be opened at low dielectric contrast ([@R12]). 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 ([@R1], [@R13]--[@R16]). 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 ([@R1]), making their experimental realization highly challenging. An alternative approach is to use driven systems, e.g., self-organized colloidal suspensions under periodic shearing ([@R17]), to form nonequilibrium dynamic hyperuniform states, which can effectively avoid the dynamic trapping, and, in principle, have a higher resistance to thermal perturbations ([@R7]). Now, experimentalists only succeeded in generating weakly hyperuniform structures with λ ≃ *d* + 0.5 ([@R8]) using this method. Nevertheless, this is certainly a direction that deserves further investigation ([@R10]), as some self-driven systems, or active matter systems, have produced a number of unexpected emergent phenomena never found in corresponding equilibrium systems ([@R18]--[@R22]). However, in conventional active matter systems, i.e., active nematic systems and active Brownian particles systems, giant number fluctuations characterized by λ \< *d* ([@R19]) and motility-induced phase separation (MIPS) with λ ≃ 0 ([@R20]--[@R23]) 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) ([@R24]--[@R28]) or active spinner fluids ([@R29]), has attracted considerable attention. Both experiments and simulations have shown many interesting collective phenomena in these systems ([@R25], [@R30]--[@R33]). 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.

RESULTS
=======

Model
-----

As illustrated in [Fig. 1A](#F1){ref-type="fig"}, we consider a 2D suspension of *N* active colloidal particles with diameter σ. Each particle experiences an in-plane force *F*^p^ with random initial orientation and a constant torque Ω perpendicular to the plane, which drive the particles to perform circular motion with the same handiness ([@R30], [@R31]). The dynamics of particle *i* at finite temperature *T* is governed by the over-damped Langevin equations ([@R26], [@R30])$${\overset{.}{\mathbf{r}}}_{\mathit{i}}(\mathit{t}) = \gamma_{\mathit{t}}^{- 1}\lbrack - \nabla_{\mathit{i}}\mathit{U}(\mathit{t}) + \mathit{F}^{p}\mathit{e}_{\mathit{i}}(\mathit{t})\rbrack + \sqrt{2\mathit{k}_{B}\mathit{T}/\gamma_{t}}\xi_{\mathit{i}}^{\mathit{t}}(\mathit{t})$$$${\overset{.}{\mathbf{e}}}_{\mathit{i}}(\mathit{t}) = \left\lbrack \gamma_{r}^{- 1}\Omega + \sqrt{2\mathit{k}_{B}\mathit{T}/\gamma_{r}}\xi_{\mathit{i}}^{r}(\mathit{t}) \right\rbrack \times \mathbf{e}_{\mathit{i}}(\mathit{t})$$where **r**~*i*~ and **e**~*i*~ are the position of particle *i* and its self-propulsion orientation, respectively, and *k*~B~ is the Boltzmann constant. γ~t/r~ is the translational/rotational friction coefficient. For simplicity, we set γ~t~ = γ~r~/σ^2^. $\xi_{\mathit{i}}^{t}(\mathit{t})$ and $\xi_{\mathit{i}}^{r}(\mathit{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 $\mathit{v}_{0} = \gamma_{\mathit{t}}^{- 1}\mathit{F}^{p}$. The reduced noise strength in the system is defined as *T*~R~ = *k*~B~*T*/(*F*^p^σ) that measures the strength of thermal noise compared with the self-propulsion. In the zero noise limit, i.e., *T*~R~ = 0, isolated active particles perform circular motions with fixed radius *R* = *F*^p^σ^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.

![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 *T*~R~ =0.](aau7423-F1){#F1}

Absorbing-active transition
---------------------------

We first simulate systems with *N* = 10,000 and *T*~R~ = 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) 〈Δ*r*^2^〉 of particles develops a plateau at long time (blue line in [Fig. 1E](#F1){ref-type="fig"}). 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 〈Δ*r*^2^〉 \~ 4*Dt* (*t* → ∞) (red line in [Fig. 1E](#F1){ref-type="fig"}). Here, *D* is the long-time diffusion constant. The phase behaviors of the system are summarized in the phase diagram [Fig. 1B](#F1){ref-type="fig"}. In the following, we focus on the absorbing-active transition close to the boundary of two phases.

In [Fig. 1F](#F1){ref-type="fig"}, 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](#F1){ref-type="fig"}. 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](#F1){ref-type="fig"}, we show the log-log plot of *D* as a function of ϕ − ϕ~c~. The obtained slope β is given in [Fig. 1J](#F1){ref-type="fig"}. 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) ([@R6]--[@R34]). 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 ([@R35]). 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](#F1){ref-type="fig"}, 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](#F1){ref-type="fig"} (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](#F1){ref-type="fig"}. 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*) = 〈∑~*i*≠*j*~ **e**~*i*~ ⋅ **e**~*j*~ δ(*r*~*ij*~ − *r*)〉/ρ*N*, as shown in [Fig. 1](#F1){ref-type="fig"} (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) ([@R26]). 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 ([@R5], [@R6], [@R17], [@R34]), 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](#F1){ref-type="fig"}, 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 ([@R5], [@R6]). However, for relative large *R*, the critical *q* scaling shifts to a faster *q*^2^ scaling (see fig. S2). To explore this further, we simulate a system of *N* = 40,000 circle active particles. In [Fig. 2A](#F2){ref-type="fig"}, 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](#F2){ref-type="fig"}, 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*) \~ *q*^2^(*q* → 0). From 〈δρ^2^〉, one can identify an *R*-dependent length scale *L*~HU~, above which the system becomes hyperuniform, while below which the system behaves like normal fluids, i.e., 〈δρ^2^〉 \~ *L*^−2^ ([Fig. 2B](#F2){ref-type="fig"}). In *S*(*q*), a similar threshold *q*~HU~ can be found, which features the end of the hyperuniform scaling. With increasing *R*, *L*~HU~ increases and *q*~HU~ 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.

![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 *L*^0^ is for clustering- or phase separation--induced large density fluctuations. (**C** and **F**) Structure factor *S*(*q*) for various *R*. The *q*^2^ 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 *T*~R~.](aau7423-F2){#F2}

[Figure 2](#F2){ref-type="fig"} (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](#F2){ref-type="fig"}). The high-density systems also show clear hyperuniform scaling at large length scales for *R* ≤ 50σ and the threshold *L*~HU~ (or *q*~HU~) increases (or decreases) with increasing *R* ([Fig. 2](#F2){ref-type="fig"}, E and F). However, at ϕ = 0.4, we find 〈δρ^2^〉 \~ *L*^−λ^ with λ \< 2 when *L* ≪ *L*~HU~ and λ decreases with larger *R*. This implies that at length scales *L* ≪ *L*~HU~, 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* \> *q*~HU~ as shown in [Fig. 2F](#F2){ref-type="fig"}, which was reported as a signature of critical instability of active particle system undergoing MIPS ([@R20]). However, in most of our systems except *R* ≥ 100σ, *S*(*q*) \~ *q*^−2^ does not diverge at *q*~HU~ ≃ 0 as in MIPS ([@R20]) but stops at a finite *q*~HU~. 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 *q*~HU~, 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 *q*~HU~ 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](#F2){ref-type="fig"} (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 ([@R20], [@R21]), one can prove (section S1) that in a homogeneous circle active particles system with vanishing orientation order parameter $\mathbf{Q} = \langle\mathbf{e}\mathbf{e}^{\mathit{T}} - \frac{1}{2}1\rangle$, the time-dependent local density field ρ(**r**, *t*) and the local polarization field **p**(**r**, *t*) satisfy$$\partial_{t}\rho = - \nabla \cdot \lbrack\mathit{v}_{e}(\rho)\mathbf{p} - \mathit{D}_{e}\nabla\rho\rbrack$$$$\partial_{t}\mathbf{p} = - \frac{1}{2}\nabla\lbrack\mathit{v}_{e}(\rho)\rho\rbrack + \mathit{D}_{e}\nabla^{2}\mathbf{p} + \Omega_{r} \times \mathbf{p}$$where *v*~e~(ρ) = *v*~0~ + ζρ is a density-dependent effective velocity of particles with a negative ζ reflecting the motility-induced "self-trapping" effect. *D*~e~ is the effective diffusion constant originated from the "evasive" motion of particles due to the collisions with neighboring particles ([@R21]). $\Omega_{r} = \gamma_{r}^{- 1}\Omega$ is the reduced torque. The isotropic homogeneous state, i.e., $\lbrack\rho(\mathbf{r},\mathit{t}) = \overline{\rho},\mathbf{p}(\mathbf{r},\mathit{t}) = 0\rbrack$, is the solution to [Eqs. 3](#E3){ref-type="disp-formula"} and [4](#E4){ref-type="disp-formula"}. By making a weak perturbation around this state, i.e., $\lbrack\rho(\mathbf{r},\mathit{t}) = \overline{\rho} + \delta\rho(\mathbf{r},\mathit{t}),\mathbf{p}(\mathbf{r},\mathit{t}) = \delta\mathbf{p}(\mathbf{r},\mathit{t})\rbrack$, we obtain two linearized equations in the Fourier space with the first-order approximation$$(\mathit{i}\omega + \mathit{D}_{e}\mathit{q}^{2})\delta\rho_{\mathbf{q},\omega} = - \mathit{i}\mathit{v}_{e}\mathbf{q} \cdot \mathbf{p}_{\mathbf{q},\omega}$$$$(\mathit{i}\omega + \mathit{D}_{e}\mathit{q}^{2})\mathbf{p}_{\mathbf{q},\omega} = \Omega_{\mathit{r}} \times \mathbf{p}_{\mathbf{q},\omega} - \mathit{i}\mathit{w}\mathbf{q}\delta\rho_{\mathbf{q},\omega}$$where \[δρ~**q**,ω~, **p~q,ω~**\] = ∫*d***r** ⋅ *e*^−*i***q**⋅**r**^∫*dt e*^−*i*ω*t*^\[δρ, **p**\], and $\mathit{w} = \mathit{v}_{0}/2 + \zeta\overline{\rho}$ is the parameter indicating the strength of self-trapping effect. By solving [Eqs. 5](#E5){ref-type="disp-formula"} and [6](#E6){ref-type="disp-formula"}, one obtains the dispersion relationship, which includes a diffusive mode ω~0~ = *iD*~e~*q*^2^ and two nondiffusive modes$$\omega_{1,2} = \mathit{i}\mathit{D}_{e}\mathit{q}^{2} \pm \sqrt{\mathit{v}_{e}\mathit{w}\mathit{q}^{2} + \Omega_{r}^{2}}$$

The growth rate of the mode is *κ* = *Re*(*i*ω), whose sign determines whether the perturbation δρ \~ *e*^*i*ω*t*+*i***q**⋅**r**^ grows or decays. One can prove that the mode 1 always decays, while the mode 2 may grow for *w* \< 0 with$$\kappa_{2} = \begin{cases}
{- \mathit{D}_{e}\mathit{q}^{2}} & {\mathit{q} < \frac{\mathit{v}_{0}}{\sqrt{- \mathit{v}_{e}\mathit{w}}}\mathit{R}^{- 1}} \\
{- \mathit{D}_{e}\mathit{q}^{2} + \sqrt{- \mathit{v}_{e}\mathit{w}\mathit{q}^{2} - \Omega_{r}^{2}}} & {\mathit{q} > \frac{\mathit{v}_{0}}{\sqrt{- \mathit{v}_{e}\mathit{w}}}\mathit{R}^{- 1}} \\
\end{cases}$$

In [Fig. 3A](#F3){ref-type="fig"}, 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 −*D*~e~*q*^2^ from *q* = 0 to $\frac{\mathit{v}_{0}}{\sqrt{- \mathit{v}_{e}\mathit{w}}}\mathit{R}^{- 1}$, 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 ([@R20], [@R21]).

![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](#E9){ref-type="disp-formula"} 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).](aau7423-F3){#F3}

The instability point of the system is defined as $\kappa_{2}^{\text{max}}(\mathit{q}*) = 0$ for *q*\* \> 0, from which we can obtain the relationship between the critical packing fraction ϕ\* and *q*\* (see section S2)$$\left( \frac{\mathit{\phi}*}{\phi_{c}} - 1 \right)\left( 2 - \frac{\mathit{\phi}*}{\phi_{c}} \right) = \frac{{8\mathit{D}}_{e}}{\mathit{v}_{0}\mathit{R}}$$$$\mathit{q}* = \sqrt{\frac{\mathit{v}_{0}}{{\mathit{R}\mathit{D}}_{e}}}$$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](#F3){ref-type="fig"}, 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](#E9){ref-type="disp-formula"} (the dotted line in [Fig. 3B](#F3){ref-type="fig"}), which leads to *D*~*e*~ ≃ 0.62*v*~0~*σ*. In the inset of [Fig. 3B](#F3){ref-type="fig"}, we plot the location of *q*\* measured in simulations and the prediction of [Eq. 10](#E10){ref-type="disp-formula"} with the same *D*~e~. 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](#E7){ref-type="disp-formula"} and [8](#E8){ref-type="disp-formula"}). 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 **r**~*i*~, the circling center is at $\mathbf{r}_{\mathit{i}}^{o} = \mathbf{r}_{\mathit{i}} - \sigma^{2}(\mathbf{F}_{\mathit{i}}^{\mathit{p}} \times \Omega)/{|\Omega|}^{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 ([@R5], [@R6], [@R10]), 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](#F3){ref-type="fig"}). 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*~ρ = *D*∇^2^ρ in the *q* space$$\partial_{t}\rho_{\mathbf{q}}^{o} = - \mathit{D}_{e}^{o}\mathit{q}^{2}\rho_{\mathbf{q}}^{o} + \xi_{\mathbf{q}}(\mathit{t})~\text{for}~\left( \mathit{q} \ll \frac{2\pi}{\mathit{R}} \right)$$with $\lbrack\rho_{\mathbf{q}}^{o},\xi_{\mathbf{q}}\rbrack = \int\mathit{d}\mathbf{r} \cdot \mathit{e}^{- \mathit{i}\mathbf{q} \cdot \mathbf{r}}\lbrack\rho^{o},\xi\rbrack$. Here, $\mathit{D}_{\mathit{e}}^{o}$ 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 ([@R36]), which obeys the center of mass conservation (CMC). As proven in ([@R10]), the noise with additional CMC appears as a double-space derivative in the diffusion equation, i.e., $\xi(\mathit{t}) = \sqrt{\overline{\rho}}\nabla^{2}\eta(\mathit{t})$, with 〈η(**r**, *t*)η(**r**′, *t*′)〉 = *A*^2^δ(**r** − **r**′)δ(*t* − *t*′) 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 ([@R18], [@R37]). Following ([@R10]), by adding such a double-space derivative noise term into [Eq. 11](#E11){ref-type="disp-formula"} as the perturbation, one obtains the dynamic equation for the density fluctuation δρ^o^ in the Fourier space$$\mathit{i}\omega\delta\rho_{\mathbf{q},\omega}^{o} = - \mathit{D}_{e}^{o}\mathit{q}^{2}\delta\rho_{\mathbf{q},\omega}^{o} - \mathit{q}^{2}\sqrt{\overline{\rho}}\eta_{\mathbf{q},\omega}~\text{for}~\left( \mathit{q} \ll \frac{2\pi}{\mathit{R}} \right)$$where η~**q**,ω~ is the Fourier transform of the noise η(**r**, *t*). From [Eq. 12](#E12){ref-type="disp-formula"}, we obtain the structure factor$$\begin{matrix}
{\mathit{S}^{o}(\mathbf{q}) = \int_{- \infty}^{\infty}\frac{1}{2\pi\tau_{\text{max}}\mathit{N}}\langle\delta\rho_{\mathbf{q},\omega}^{o}\delta\rho_{\mathbf{q},\omega}^{o*}\rangle\mathit{d}\omega} \\
{= \frac{\mathit{A}^{2}}{{2\mathit{D}}_{e}}\mathit{q}^{2}~\text{for}~\left( \mathit{q} \ll \frac{2\pi}{\mathit{R}} \right)} \\
\end{matrix}$$where τ~max~ is the maximum observation time (see section S3). [Equation 13](#E13){ref-type="disp-formula"} 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* ([@R2]). In [Fig. 4A](#F4){ref-type="fig"}, we plot the scaled density variance 〈δρ^o2^〉*R*^2^ 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* \> *L*~HU~ ≃ 2*R* and a normal fluid-like scaling *L*^−2^ at *L* \< *L*~HU~. Similar collapse is shown in [Fig. 4B](#F4){ref-type="fig"} for *S*^o^(*q*) versus *qR*, where the hyperuniform scaling *S*^o^(*q*) \~ *q*^2^ stops at *q*~HU~ ≃ 2π/*L*~HU~. For high-density systems, i.e., ϕ = 0.4, we obtain the same crossover length scale *L*~HU~ ≃ 2*R* as shown in [Fig. 4](#F4){ref-type="fig"} (C and D). In this case, the density variance and *S*^o^(*q*) for different *R* do not collapse into a single curve because of cluster formation--induced large density fluctuations at *L* ≲ *L*~HU~. 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 ([@R38]).

![Global hyperuniformity.\
(**A** and **C**) Scaled density variance 〈δρ^o2^〉*R*^2^ as functions of *L*/*R* and (**B** and **D**) *S*^o^(*q*) as functions of *qR*, for various *R* at ϕ = 0.2 (A and B) and 0.4 (C and D), respectively.](aau7423-F4){#F4}

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 ([@R10]), 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 ([@R10]).

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 ([@R26]). 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 ([@R26]). 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 ([@R26]). 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 ([@R26]). In the active mixing state, we observe the same hyperuniform scaling *S*(*q* → 0) \~ *q*^2^. 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 ([@R29]) or size/interaction oscillating particle systems ([@R39]), 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 relationship$$\mathit{S}(0) = \kappa_{T}\rho\mathit{k}_{B}\mathit{T}$$any thermal equilibrated systems with the positive isothermal compressibility κ~T~ at finite temperatures cannot be strictly hyperuniform due to thermal excitation ([@R13]). 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 ([@R13], [@R40]). 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 ([@R1], [@R13], [@R15]) in disordered structures may also be enough to achieve some desired functions, e.g., isotropic photonic/electronic bandgaps ([@R14], [@R15]). In [Fig. 5](#F5){ref-type="fig"}, 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 *T*~R~ from zero, *S*(*q* → 0) begins to saturate at some nonzero value, which increases along with *T*~R~ ([Fig. 5A](#F5){ref-type="fig"}). In section S4, we introduce the thermal noise **f** = \[*f*~*x*~, *f*~*y*~\] in [Eq. 11](#E11){ref-type="disp-formula"} as$$\xi(\mathit{t}) = \sqrt{\overline{\rho}}\nabla \cdot \lbrack\nabla\eta(\mathit{t}) + \mathbf{f}(\mathit{t})\rbrack$$where 〈*f*~*i*~(**r**, *t*)*f*~*j*~(**r**′, *t*′)〉 = 2*D*~therm~δ~*ij*~δ(**r** − **r**′)δ(*t* − *t*′) with $\mathit{D}_{\text{therm}} = \mathit{k}_{B}\mathit{T}\gamma_{t}^{- 1}(1 + \mathit{R}^{2}/\sigma^{2})$ the self-diffusion constant of effective particles due to the thermal Brownian motion ([@R25]). 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: 〈*f*~*i*~(**r**, *t*)η(**r**′, *t*′)〉 = 0. Then, we can estimate *S*^o^(*q*) as a function of the reduced noise strength *T*~R~ for low-density systems as$$\mathit{S}^{o}(\mathit{q}) = \frac{\mathit{v}_{0}\sigma}{\mathit{D}_{e}}\left( 1 + \frac{\mathit{R}^{2}}{\sigma^{2}} \right)\mathit{T}_{R} + \frac{\mathit{A}^{2}}{{2\mathit{D}}_{e}}\mathit{q}^{2}~\text{for}~(\mathit{q} < \mathit{q}_{\text{HU}})$$

![Effect of thermal noises on hyperuniformity.\
(**A**) Structure factor *S*^o^(*q*) under different reduced noise strength *T*~R~ for ϕ = 0.2 and *R* = 3σ. Open symbols show the simulation data, while the dashed lines are the theoretical predictions of [Eq. 16](#E16){ref-type="disp-formula"}. (**B**) Normalized *S*^o^(0) as functions of *T*~R~ 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.](aau7423-F5){#F5}

In [Fig. 5A](#F5){ref-type="fig"}, we plot the theoretical prediction of [Eq. 16](#E16){ref-type="disp-formula"} as dashed lines for different *T*~R~ values at *R* = 3σ and ϕ = 0.2, by assuming them all across the same *S*^o^(*q*~HU~) point. In [Fig. 5B](#F5){ref-type="fig"}, we compare the *S*^o^(0) from the theoretical prediction with the *S*^o^(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](#F5){ref-type="fig"}, 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](#E14){ref-type="disp-formula"}. However, we also emphasize the difference: In our nonequilibrium system, the saturated value of *S*^o^(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., *T*~R~ ≪ 1.

DISCUSSION
==========

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 *L*~HU~. For *L* ≫ *L*~HU~, the system exhibits strong hyperuniformity with the density variance scaling the same as in perfect crystals, while for *L* ≲ *L*~HU~, 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 ([@R10]). 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.

METHODS
=======

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~ = σ/*v*~0~. To mimic the excluded volume interaction between colloidal particles *i* and *j*, we used the WCA potential$$\mathit{U}(\mathit{r}_{\mathit{i}\mathit{j}}) = \begin{cases}
{4\varepsilon\left\lbrack \left( \frac{\sigma}{\mathit{r}_{\mathit{i}\mathit{j}}} \right)^{12} - \left( \frac{\sigma}{\mathit{r}_{\mathit{i}\mathit{j}}} \right)^{6} + \frac{1}{4} \right\rbrack} & {(\mathit{r}_{\mathit{i}\mathit{j}} < 2^{1/6}\sigma)} \\
0 & {(\mathit{r}_{\mathit{i}\mathit{j}} > 2^{1/6}\sigma)} \\
\end{cases}$$where *r*~*ij*~ is the center to center distance between particles *i* and *j*, with σ as the diameter of particles. We chose ε = *F*^p^σ/24, which gives the typical contact distance *r*~c~ = σ between particles based on force balance $\mathit{F}^{p} = {\frac{\partial\mathit{U}(\mathit{r}_{\mathit{i}\mathit{j}})}{{\partial\mathit{r}}_{\mathit{i}\mathit{j}}}|}_{\mathit{r}_{\mathit{i}\mathit{j}} = \mathit{r}_{c}}$. For systems at *T*~R~ = 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 ([@R26]), we added an additional soft repulsion *U*~s~ = *A* (*r*/σ)^−4^ with cutoff distance *r*~sc~ = 5σ and *A* = 7.5ε to model the isotropic dipole interaction between active particles. The self-propulsion force for this system was reset to *F*^p^ = 2^7/6^*A*/σ to make the dipolar interaction balance the driving force at *r*~c~ = 2^1/6^σ ([@R26]).

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., $\lbrack\mathit{q}_{\mathit{x}},\mathit{q}_{\mathit{y}}\rbrack = \frac{2\pi}{\mathit{L}_{0}}\lbrack\mathit{i},\mathit{j}\rbrack(\mathit{i},\mathit{j} = 1,2,3\ldots)$, with *L*~0~ 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 *q*~min~ = 2π/*L*~0~ without further change. The fitting function used in [Fig. 5B](#F5){ref-type="fig"} is [Eq. 16](#E16){ref-type="disp-formula"}, with an adjustable *T*~R~.

Supplementary Material
======================

###### http://advances.sciencemag.org/cgi/content/full/5/1/eaau7423/DC1

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.

Supplementary material for this article is available at <http://advances.sciencemag.org/cgi/content/full/5/1/eaau7423/DC1>

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 *S*^o^(*q*) for system with CMC.

Section S4. Effect of thermal noise on *S*^o^(*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 *T*~R~ = 0.

Fig. S4. Hyperunifomity in an experimentally realizable system ([@R26]) 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.
