Hydrodynamics of Random-Organizing Hyperuniform Fluids

TL;DR

This paper develops a hydrodynamic theory for hyperuniform fluids in non-equilibrium systems, revealing how damping and activity lead to suppressed density fluctuations and hyperuniformity, with potential experimental realizations.

Contribution

It introduces a hydrodynamic framework for understanding fluidic hyperuniformity in driven-dissipative systems, linking damping and activity to hyperuniform density fluctuations.

Findings

Fluidic hyperuniformity arises from damping of stochastic oscillations in density fluctuations.

Transition from absorbing to hyperuniform to equilibrium fluid controlled by dissipation strength.

Proposes experimental realization in active spinner systems.

Abstract

Disordered hyperuniform structures are locally random while uniform like crystals at large length scales. Recently, an exotic hyperuniform fluid state was found in several non-equilibrium systems, while the underlying physics remains unknown. In this work, we propose a non-equilibrium (driven-dissipative) hard-sphere model and formulate a hydrodynamic theory based on Navier-Stokes equations to uncover the general mechanism of the fluidic hyperuniformity (HU). At a fixed density, this model system undergoes a smooth transition from an absorbing state to an active hyperuniform fluid, then to the equilibrium fluid by changing the dissipation strength. We study the criticality of the absorbing phase transition. We find that the origin of fluidic HU can be understood as the damping of a stochastic harmonic oscillator in $q$ space, which indicates that the suppressed long-wavelength density fluctuation in the hyperuniform fluid can exhibit as either acoustic (resonance) mode or diffusive (overdamped) mode. Importantly, our theory reveals that the damping dissipation and active reciprocal interaction (driving) are two ingredients for fluidic HU. Based on this principle, we further demonstrate how to realize the fluidic HU in an experimentally accessible active spinner system and discuss the possible realization in other systems.

Tables

Table 1. Table 1: Comparison of critical exponents measured in 2D systems with those from universality classes of DP, CDP ( 33 , 36 ) and their long-range mean-field values.

	$\beta$	$\alpha$	${\nu }_{⟂}^{\ast }$	$z^{\ast }$
DP	0.583(3)	0.450(1)	0.733(8)	1.766(2)
CDP	0.64(1)	0.52(1)	0.80(2)	1.53(2)
Mean-field (long-range)	1	1	1	1
${\tilde{\rho }}_c=0.12192\left(2\right)$	0.65(2)	0.54(2)	0.84(3)	1.49(3)
${\tilde{\rho }}_c=3.2350\left(5\right)\times {10}^{-5}$	1.02(3)	1.03(3)	0.98(3)	1.00(2)

Equations

$$m\frac{d{\mathbf{v}}_i(t)}{dt}=-\gamma {\mathbf{v}}_i(t),$$

$$l_d\equiv \sqrt{m\Delta E}/\gamma ,$$

$$W_{driv}\simeq f_a\Delta E,\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}{W}_{disp}\simeq {\overline{v}}^2\gamma ,$$

$$k_B{T}_k\equiv \frac{m\overline{v^2}}{d}\simeq \frac{m{\overline{v}}^2}{d}\simeq \frac{\Delta E}{4d}{\left(\frac{l_d}{l_m}\right)}^2.$$

$$\frac{\partial \rho }{\partial t}$$

$$\frac{\partial (\rho \mathbf{u})}{\partial t}+\nabla \cdot (\rho \mathbf{uu})$$

$$q^{-2}\frac{{\partial }^2\delta {\rho }_{\mathbf{q}}}{\partial t^2}=-(\tilde{\gamma }{q}^{-2}+{\nu }^{\parallel })\frac{\partial \delta {\rho }_{\mathbf{q}}}{\partial t}-c_s^2\delta {\rho }_{\mathbf{q}}+{\sigma }_{\parallel ,\mathbf{q}}^r,$$

$$⟨\frac{1}{2}{c}_s^2|\delta {\rho }_{\mathbf{q}}{|}^2⟩=\frac{k_B{T}_{bath}}{2}.$$

$$S(\mathbf{q})=\frac{1}{N{m}^2}⟨|\delta {\rho }_{\mathbf{q}}{|}^2⟩=\frac{q^2}{B{l}_d^{-2}+q^2},$$

$$\frac{c_s-\sqrt{c_s^2-\tilde{\gamma }{\nu }^{\parallel }}}{{\nu }^{\parallel }}<q<\frac{c_s+\sqrt{c_s^2-\tilde{\gamma }{\nu }^{\parallel }}}{{\nu }^{\parallel }}.$$

$$\frac{\partial \delta {\rho }_{\mathbf{q}}}{\partial t}=-D{q}^2\delta {\rho }_{\mathbf{q}}+q^2{\sigma }_{\parallel ,\mathbf{q}}^r.$$

$$S(\mathbf{q},\omega )$$

$$S_R(f)\equiv \frac{1}{m^2{\tau }_{max}}{\left|{\int }_0^{{\tau }_{max}}{\rho }_R(t)e^{2i\pi ft}dt\right|}^2,$$

$$S_R(f)\text{leavevmode}\sim \text{leavevmode}{f}^{-2}.$$

$$S_R(f)\sim \text{leavevmode}{f}^{-1/2}.$$

$$S_R(f)\sim f^0.$$

Full text

\dates

This manuscript was compiled on

\dropcap

Hyperuniform structures are characterized by vanishing long-wavelength density fluctuations with the structure factor $S(q\rightarrow 0)=0$ (1, 2), which include perfectly ordered structures, such as athermal crystals or quasi-crystals. In past decades, an increasing number of disordered structures were also found to exhibit hyperuniformity (HU), e.g., maximally random jammed states (3), animal photoreceptor cell patterns (4, 5), low-discrepancy sequences (6), eigenvalues of random matrices (7), early universe fluctuations (8) etc. (2). These disorder structures have shown even better properties than crystals usually exclusively have, like isotropic photonic bandgaps opened at low dielectric contrast (9, 10, 11), abnormal transparency (12, 13) and also has many applications in graphic anti-aliasing (14) and quasi-Monte Carlo sampling (15). Disordered HU implies a long-range direct correlation function in particle systems (1, 2). In equilibrium, this requires delicately designed long-range interactions (12, 16), which are challenging for experimental realization. Recently, HU has been found in some dynamic states of non-equilibrium systems, including the non-ergodic dynamic critical states (17, 18, 19, 20, 21, 22) and ergodic fluid states (23, 24). The latter shows the promise of realizing novel flowing functional materials, e.g., photonic fluids (25). Current studies of dynamic HU mainly focus on time-discrete random-organizing models, in which particles interacting with their neighbours perform activated random displacement (26, 17, 18, 19, 23, 27). However, the general mechanism of the non-equilibrium hyperuniform fluids and the fundamental difference from equilibrium fluids remains unknown, let alone a full hydrodynamic description of these exotic fluids. To address these questions, in this work, we propose a non-equilibrium hard-sphere fluid model system, which can smoothly transform from an equilibrium simple fluid to a non-equilibrium hyperuniform fluid. We construct a hydrodynamic theory for the hyperuniform fluids based on the generalized Navier-Stokes equations. The theory reveals the mechanism of fluidic HU intuitively through a damped stochastic harmonic oscillator in $q$ space. Based on the analysis, we obtain two important ingredients for the fluidic HU and further show the realization in an experimentally accessible 2D active spinner system. Our findings suggest possible ways to realize fluidic photonic materials that can heal from damage, and their optical properties can be externally controlled (28).

Results

Model

We first consider a minimal model consisting of $N$ hard spheres (See Fig. 1A). Each particle has a mass $m$, a diameter $\sigma$ and a random initial velocity. Particles undergo active collision with reciprocal center-to-center impulsions, in which an additional kinetic energy $\Delta E$ is injected through each collision between two particles. The kinetic energy is then dissipated through the (frictional) damping, and the equation motion of particle $i$ between two consecutive collisions is

[TABLE]

with $\gamma$ the damping coefficient. This model is a generalization of previous time-discrete random-organizing models (17, 18, 19, 23, 27) by incorporating the inertia of the particles and it has two characteristic lengths: i) the mean free path $l_{m}$, which at low density is proportional to ${\sigma}/{\tilde{\rho}}$, with $\tilde{\rho}=N\sigma^{d}/V$ the dimensionless particle number density, $V$ the volume of the system and $d$ the dimension of the system; ii) the dissipation length $l_{d}$, i.e., the typical distance that an isolated activated particle can travel,

[TABLE]

which represents the strength of inertia of the particles. Here we assume two particles collide with vanishing approaching velocities and $\sqrt{\Delta E/m}$ is the typical activated velocity to each particle. As our system only has two energy (time) -related quantities, i.e., $\Delta E$ and $\gamma$, the typical dissipation time associated with $\gamma$ can be written as $\tau_{d}=m/\gamma=(l_{d}/\sigma)\tau_{0}$ by assuming $\tau_{0}\equiv\sqrt{m\sigma^{2}/\Delta E}$ the unit of time. Therefore, $l_{m}$ and $l_{d}$ uniquely define the state of the system. We devise a new event-driven algorithm to simulate this underdamped system in both 2D and 3D with full periodic boundary conditions (see Method for the simulation details). As the conclusion we obtain is independent of dimensionality, in the following, we focus on the 2D systems and leave the results of 3D systems in SI Appendix.

Absorbing Phase Transition

In our model with random initial particle velocities, when $l_{d}\ll l_{m}$, a single active collision on average induces less than one further collision, which causes the exponential decay of activation events as a function of time. The system thus falls into an absorbing state with vanishing particle velocity, in some sense analogous to the liquid-solid transition (see SI Appendix Movie S1). On the contrary, when $l_{d}\gg l_{m}$, a single activation can trigger chain reactions, leaving the system in a dynamically active state (see SI Appendix Movie S2, S3). Similar absorbing phase transitions have been reported in many other systems (26, 29, 30, 31), which have deep connections with the celebrated self-organized criticality (32). In this work, we focus on the low density regime, to avoid the possible crystallization. In Fig. 1B, we plot the phase diagram of the 2D system in the representation of $\sigma/l_{d}$ and $\tilde{\rho}$. The solid line with open symbols shows the phase boundary of system obtained from simulations, and the dashed line shows $\sigma/l_{m}$ as a function of $\tilde{\rho}$. Two lines are rather close to each other especially at small $\tilde{\rho}$, indicating the transition occurs round $l_{d}\simeq l_{m}$. In SI Appendix, Fig. S1-2, by using the finite-size scaling analysis (33, 34), we investigate the criticality of this transition in 2D systems with weak inertia ($l_{d}=\sqrt{10}\sigma$) and strong inertia ($l_{d}=10^{4}\sigma$). We summarize the obtained critical exponents in Table 1 and compare them with those from the universality of directed percolation (DP) (33) and conserved directed percolation (CDP) (35). We find the transition belongs to CDP for systems with small inertia and relative high critical density $\tilde{\rho}_{c}=0.12192$. Increasing the inertia to $l_{d}=10^{4}\sigma$ makes the system approach the dilute limit ($\tilde{\rho}_{c}=3.2350\times 10^{-5}$) and induces a crossover from CDP to its long-range mean-field scenario due to the corresponding increase of $l_{m}$. Here, $l_{m}$ acts as the effective “interaction range" in the system, which diverges as the system approaches the dilute limit. Similar results for 3D systems can be found in SI Appendix, Fig. S3-4 and Table S1.

From Equilibrium Simple Fluids to Non-Equilibrium Hyperuniform Fluids

Our model is a driven-dissipative system. In the active state, the power of energy injection and dissipation per particle at the mean-field level can be written as

[TABLE]

respectively, where $f_{a}\equiv{\overline{v}}/2{l_{m}}$ is the average collision frequency per particle and $\overline{v}$ is the average particle speed. Since the dissipation is a quadratic function of $\overline{v}$, a driving-dissipation balanced state with $W_{disp}=W_{driv}$ can be reached, which leads to the kinetic temperature of the system:

[TABLE]

In Fig. 2A, we plot $k_{B}T_{k}/\Delta E$ measured in the simulation for 2D systems with different $l_{d}$ at $l_{m}=4.4\sigma$ ($\tilde{\rho}=0.1$) (solid symbol) and the systems with different $l_{m}$ at fixed $l_{d}=160\sigma$ (open symbols). We find the two sets of data collapsing onto a single curve when $l_{d}\gg l_{m}$, which can be written as $\frac{k_{B}T_{k}}{\Delta E}=\frac{A}{8}\left(\frac{l_{d}}{l_{m}}\right)^{2}$, with the pre-factor $A=1.56$. The distributions of particle velocity in $x,y$ directions, i.e., $P(v_{x})$, $P(v_{y})$, for systems at $l_{m}=4.4\sigma$ with different $k_{B}T_{k}/\Delta E$ are shown in Fig. 2B. We find that $P(v_{x})$ and $P(v_{y})$ are statistically indistinguishable from each other. With increasing $k_{B}T_{k}/\Delta E$, the velocity distribution gradually changes from non-Gaussian distribution (a signature of non-equilibrium fluids) to the Maxwell-Boltzmann (Gaussian) distribution: that is, $P(v_{i})=\left(\frac{m}{2\pi k_{B}T_{k}}\right)^{1/2}\exp\left[\frac{-mv_{i}^{2}}{2k_{B}T_{k}}\right],\leavevmode\nobreak\ (i=x,y)$, with which the system behaves essentially the same as an equilibrium fluid at temperature $T_{k}$. Such a behavior is independent of the density or $l_{m}$ (SI Appendix, Fig. S5). In Fig. 2C, we plot the structure factor $S(q)$ of systems at $l_{m}=4.4\sigma$ with different $l_{d}$. When $l_{d}=96\sigma$, we find $S(q)\sim 1$, close to that of low density equilibrium simple fluids (open symbol). With decreasing $l_{d}$, the $S(q)$ curve gradually bends down at small $q$ regime and develops a deep $q^{2}$ scaling, indicating a strong HU (23, 24). Therefore, at a fixed density, the strength of inertia represented by $l_{d}$ determines not only how far the system is out of equilibrium, but also the length scale, above which the hyperuniform scaling regime begins.

Hydrodynamic Theory for Hyperuniform Fluids

The random-organizing hard-sphere fluid has a well-defined kinetic temperature in the large $l_{d}/l_{m}$ limit. This motivates us to construct a fluctuating hydrodynamic theory for these hyperuniform fluids based on Navier-Stokes (NS) equations for isothermal fluids (37, 38). The modified NS equations for the local mass density $\rho$ and velocity $\mathbf{u}$ can be written as

[TABLE]

where $\tilde{\gamma}=\gamma/m$ is the reduced damping coefficient and $p$ is the local pressure. We assume a simple equation of state $p=c^{2}_{s}\rho$ for the system with $c_{s}$ the speed of sound. $\boldsymbol{\sigma}^{v}$ and $\boldsymbol{\sigma}^{r}$ are the classical momentum-conserved viscous stress tensor and random (noise) stress tensor for simple fluids, respectively (39, 40). Details about the hydrodynamic theory are given in SI Appendix. With the additional damping term $-\tilde{\gamma}\rho\mathbf{u}$, Eq. (5-6) describe a non-equilibrium fluid which violates the fluctuation-dissipation theorem (37). Using the standard hydrodynamic linearization procedure, we find the density fluctuation $\delta\rho_{\mathbf{q}}$ in $\mathbf{q}$ space satisfying

[TABLE]

with $\nu^{\parallel}$ the longitudinal kinematic viscosity. ${\sigma}^{r}_{\parallel,\mathbf{q}}$ is the longitudinal component of random noise in $\mathbf{q}$ space which satisfies $\langle{{\sigma}^{r}_{\parallel,\mathbf{q}}}(t){{\sigma}^{r}_{\parallel,\mathbf{q}}}(t^{\prime})\rangle=2\rho_{0}\nu^{\parallel}k_{B}T_{k}V\delta(t-t^{\prime})$ with $\rho_{0}$ the average mass density of the system. Eq. (7), in fact, describes a damped stochastic harmonic oscillator if we map $q^{-2}$ as the mass, $\delta\rho_{\mathbf{q}}$ as the displacement, $\tilde{\gamma}q^{-2}+\nu^{\parallel}$ as the effective damping coefficient, $c_{s}^{2}\delta\rho_{\mathbf{q}}$ as the restoring force and ${{\sigma}^{r}_{\parallel,\mathbf{q}}}$ as the driving noise. The average potential of the oscillator can be obtained based on the equipartition theorem (41) if we assume the damping from a viscous thermal bath at temperature $T_{bath}$, i.e.,

[TABLE]

Here, according to the Langevin equation, $T_{bath}=NmT_{k}\nu^{\parallel}/(\tilde{\gamma}q^{-2}+\nu^{\parallel})$. The magnitude of oscillation, or the static structure factor of the fluid, then can be written as

[TABLE]

with $B=\sqrt{{8}/{A}}$ for 2D systems. Here we simply assume $\nu^{\parallel}\simeq l_{m}\sqrt{k_{B}T_{k}/m}$ and $c_{s}\simeq\sqrt{k_{B}T_{k}/m}$ based on the kinetic theory of ideal gas, ignoring possible pre-factors in $\nu^{\parallel}$ and $c_{s}$. Eq. (9) can be also directly obtained using our hydrodynamic theory (see SI Appendix). From Eq. (8, 9), we find $S(q)=1$ independent of $q$ for low density equilibrium fluids ($\tilde{\gamma}=0$), consistent with the fact that the average potential of the oscillator at a fixed temperature is independent of its mass (Fig. 2E). For non-equilibrium fluids ($\tilde{\gamma}>0$), when $q\gg l_{d}^{-1}$, the $\tilde{\gamma}q^{-2}$ term in Eq. (7) can be neglected compared with $\nu^{\parallel}$, and the oscillator behaves the same as when $\tilde{\gamma}=0$. However, when $q\ll l_{d}^{-1}$, the $\tilde{\gamma}q^{-2}$ term dominates, leading to the decrease of oscillation magnitude according to $S(q)\simeq q^{2}l_{d}^{2}$ (hyperuniform scaling) independent of $l_{m}$. The oscillator thus has a vanishing oscillation magnitude at small $q$, i.e., $S(q\rightarrow 0)=0$, because of the infinitely large effective damping coefficient at this point (Fig. 2F). This mechanism of fluidic HU can be also understood as that the hyperuniform fluid has a $q$-dependent temperature which vanishes at $q\rightarrow 0$. Such a special temperature property may not be captured by the non-Gaussian distribution of particles’ velocity. In Fig. 2C, the results of Eq. (9) are shown as dashed lines, which well agree with the simulation results at small $q$. In SI Appendix, Fig. S6, we show the insensitivity of $S(q)$ with respect to $l_{m}$ for low density system far from critical point. Thus, only $l_{d}$ controls the length scale, above which the HU scaling occurs. In SI Appendix Fig. S7, we also show how $S(q)$ changes when the system crosses the phase boundary. We find that this $l_{d}$-dependency still holds when the 2D system approaches the critical point $l_{d}\simeq l_{m}$ exhibiting the critical hyperuniform scaling $S(q\rightarrow 0)\sim q^{0.45}$ (17, 18, 19).

Dynamically, the oscillator represented by Eq. (7) is in the resonance state in regime:

[TABLE]

This means that the fluids can support acoustic (sound) modes (37). Due to the damping from $\gamma$, these sound modes are fast modes, which decay as $e^{-\kappa t}$ with damping rate $\kappa\simeq\frac{1}{2}(\tilde{\gamma}+\nu^{\parallel}q^{2})$. For low density systems with $l_{d}\gg l_{m}$, this acoustic regime can be approximated as $\frac{Bl_{m}}{2}l_{d}^{-2}<q<2l_{m}^{-1}$. With decreasing $q$ or increasing $\gamma$, the oscillator crosses from the resonance regime to the overdamped regime, in which the system only has the diffusive mode $\omega^{*}_{1}\simeq-iDq^{2}$ with the diffusion constant $D=c_{s}^{2}/\tilde{\gamma}$. In fact, in this overdamping limit, Eq. (7) can be written as a simple diffusion equation,

[TABLE]

This equation was employed to explain the HU in a time-discrete random-organizing model with the ‘center of mass conservation’, a discrete form of reciprocal interaction in overdamped system, which breaks down in continuous space-time when particles have finite inertia (23). In Fig. 2G, we plot the schematic diagram of density modes in the representation of $1/l_{d}$ and $q$. One can see that when $l_{d}<l_{m}$, the system stays in an absorbing state without any dynamic mode (blue region), while for $l_{d}>l_{m}$, the system has the diffusive mode (red region) and acoustic modes (green region) at relative small $q$, and non-acoustic ‘fast-decay mode’ at large $q$. The upper boundary of the acoustic regime is determined by Eq. (10), which becomes less accurate (dashed line) when the system approaches the critical point $l_{d}\simeq l_{m}$. The shadow area with bottom boundary $q\sim l_{d}^{-1}$ represents the $q$ space regime exhibiting the $q^{2}$ hyperuniform scaling. This area is mainly located at the diffusive regime, but also overlaps with the acoustic regime in ${l_{m}}l_{d}^{-2}<q\ll l_{d}^{-1}$, at which $S(q)\simeq q^{2}l_{d}^{2}\rightarrow 0$ for $l_{d}\rightarrow\infty$. Therefore, the suppressed long-wavelength density fluctuation in hyperuniform fluids can exhibit as either diffusive or acoustic modes, though at infinite wavelength ($q=0$), only diffusive mode exists.

Generally, different dynamic modes can be characterized by different shapes of dynamic structure $S(q,\omega)$ with $\omega$ the angular frequency. For acoustic modes, there would be two symmetric Brillouin peaks $\omega_{\pm}=\pm c_{s}q$ in $S(q,\omega)$, while for diffusive modes, only a single peak exists round $\omega=0$. Our hydrodynamic theory predicts

[TABLE]

In Fig. 3, we show the measured $S(q,\omega)$ in simulation as a function of $\omega$ (symbols) for systems with fixed $l_{d}$ but different $q$ (Fig. 3A), and systems with fixed $q$ but different $l_{d}$ (Fig. 3B) at the same density $\tilde{\rho}=0.1$. One can clearly distinguish the acoustic modes in large $l_{d}$ (or $q$) cases and diffusive modes in small $l_{d}$ (or $q$) cases, respectively, in Fig. 3A (or Fig. 3B). The theoretical predictions of Eq. (12) are plotted as dashed lines by using the accurate expression of sound speed $c_{s}^{a}$ for the equilibrium fluid, which takes into account the effects of the adiabatic compression and finite density of system, i.e., $c_{s}^{a}=\sqrt{\Gamma/\rho_{0}\chi_{T}}$, with $\Gamma=1+2/d$ the adiabatic index for monoatomic fluids and $\chi_{T}=S_{\rm eq}(0)m(\rho_{0}k_{B}T_{k})^{-1}$ the isothermal compressibility (39). Here, $S_{\rm eq}(0)=0.72$ is obtained from the equilibrium system in Fig. 2C. We find that the theoretical predictions of Brillouin peaks match well not only with data from the equilibrium fluid but also with non-equilibrium hyperuniform fluids. This indicates that the sound speed and the compressibility in hyperuniform fluids with kinetic temperature $T_{k}$ should be close to that of the the equilibrium system with the same $T_{k}$ despite their huge difference in $S(0)$. For equilibrium fluids (energy conserved) or near-equilibrium fluids, due to the heat diffusion, the measured $S(q,\omega)$ shows an additional Rayleigh peak at $\omega=0$ (42, 39). Although our isothermal hydrodynamic theory without energy conservation can not capture this feature (see SI Appendix for the discussion), it predicts the bimodal acoustic modes and diffusive modes faithfully for hyperuniform fluids with small $l_{d}/l_{m}$ ratio where the heat diffusion is negligible. Especially, in the overdamping limit, both the simulation and theoretical results in Fig. 3B indicate $S(q,\omega)\rightarrow const.\leavevmode\nobreak\ (\omega\rightarrow 0)$, i.e., a fixed point for hyperuniform fluids, while for non-hyperuniform equilibrium Brownian particle system, $S(q,\omega)\propto\frac{q^{2}}{\omega^{2}+D^{2}q^{4}}\rightarrow(Dq)^{-2}\leavevmode\nobreak\ (\omega\rightarrow 0)$, i.e., a standard Lorentz distribution (42). The above dynamic structure features of random-organizing hyperuniform fluids facilitate their experimental detection using classical light scattering methods (42, 39).

Power Spectrum of Local Density Fluctuations

The aforementioned dynamic regimes in hyperuniform fluids can also be distinguished by the noise color reflected by the power spectrum of local density fluctuations  (32, 43), i.e.,

[TABLE]

where $f=\omega/2\pi$ is the frequency and $\tau_{max}$ is the maximum observation time. $\rho_{R}(t)\equiv mN_{R}(t)/\pi R^{2}$ is the local mass density with $N_{R}(t)$ the number of particles in a circular sub-domain of radius $R$. The calculation details and the complete expression of $S_{R}(f)$ based on our hydrodynamic theory can be found in SI Appendix. Here, we only summarize the main results. We first introduce three typical frequencies $\tilde{\gamma}$, $f_{1}$ and $f_{max}$, where $\tilde{\gamma}^{-1}$ is the typical dissipation time, $f_{1}^{-1}=R/c_{s}$ is the time with which the acoustic wave propagates a distance $R$ and $f_{max}=c_{s}/l_{m}$ is the typical cutoff frequency for the acoustic wave. Here, $\tilde{\gamma}/f_{max}\propto({l_{m}}/{l_{d}})^{2}$. For equilibrium simple fluids ($\tilde{\gamma}/f_{max}\rightarrow 0$), or more generally, $\mathrm{max}\left[f_{1},(\tilde{\gamma}f_{max})^{1/2}\right]\ll f\ll f_{max}$, the power spectrum exhibits a ‘red’ noise

[TABLE]

For overdamped hyperuniform fluids ($\tilde{\gamma}/f_{max}\sim 1$), or more generally, $\mathrm{min}(f_{1},f_{1}^{2}/\tilde{\gamma})\ll f\ll\tilde{\gamma}$, we have

[TABLE]

This special noise is different from the $f^{-3/2}$ noise found in non-hyperuniform equilibrium Brownian particle systems (32, 43). Moreover, there is also a crossover regime, i.e., $\tilde{\gamma}\ll f\ll(\tilde{\gamma}f_{max})^{1/2}$, in which the power spectrum behaves as ‘white’ noise

[TABLE]

In Fig. 2D, we plot $S_{R}(f)$ measured in the simulation for the same systems as shown in Fig. 2C but distinguish them by using $\tilde{\gamma}/f_{max}$. The theoretical predictions from our hydrodynamic theory multiplied by an adjustable coefficient $0.62$ are shown as dashed lines. We find that when $\tilde{\gamma}/f_{max}$ is small, the measured $S_{R}(f)$ shows the ‘red’ noise scaling in a broad range of frequency accompanied with some acoustic peaks around $f_{1}$. With increasing $\tilde{\gamma}/f_{max}$, the acoustic peaks fade away and the special $f^{-1/2}$ noise scaling gradually dominates for $f\ll\tilde{\gamma}$. The white noise plateau can only be seen for the cases of $\tilde{\gamma}/f_{max}=1.5\sim 6.0\times 10^{-2}$ due to the finite size of our system. Comparing Fig 2C with Fig 2D, we find that the development of hyperuniform scaling $q^{2}$ in $S(q)$ with decreasing $l_{d}$ is associated with the suppression of temporal fluctuation at small frequency regime with increasing $\tilde{\gamma}/f_{max}$. Actually, the ‘red’ noise regime corresponds to the non-overlapped acoustic region in Fig. 2G, while the white noise regime with upper-bound frequency $(\tilde{\gamma}f_{max})^{1/2}\sim c_{s}/l_{d}$ and the $f^{-1/2}$ noise regime with the upper-bound frequency $\tilde{\gamma}\sim c_{s}l_{m}/l_{d}^{2}$ correspond to the acoustic and diffusive hyperuniform regions in Fig. 2G, respectively. These special noise features provide us with an additional tool to identify and study the hyperuniform fluids experimentally in the frequency domain. We emphasize that all these structural and dynamic properties of our hyperuniform fluids are independent of dimensionality and we show the simulation data for 3D system in SI Appendix, Fig. S8, S9 and Movie S4.

Hyperuniform Active Spinner Fluids

The hydrodynamic theory above suggests two important conditions for forming non-equilibrium hyperuniform fluids: i) the damping dissipation characterized by $\tilde{\gamma}$, which has a dominant effect at small $q$, and ii) the active reciprocal interaction (driving), which conserves the classical formula of viscosity and the random noise. Based on this, we further consider a realistic system of $N$ active spinners on a frictional substrate (Fig. 4A) as an example to realize the hyperuniform fluid. Active spinners systems exhibit many interesting phenomena and have been extensively studied both theoretically (44, 45, 46, 47) and experimentally (48, 49, 50, 51, 52, 53). We employ the active spinner model from Ref. (45) (also see SI Appendix for details), in which each spinner is a dimer consisting of two spherical hard monomers. Driven by a constant torque $\Omega$, spinners perform self-rotating motion with the same chirality. The dynamics of both translational and rotational freedom degrees of spinners are underdamped. Thus, when two fast-rotating spinners collide, some part of their rotational kinetic energy is transferred to their translational freedom degree, inducing the active collision similar to the previous hard sphere model. With the fixed substrate friction coefficient, we find $\Omega$ in the spinners system plays the similar role of $\Delta E$ in the previous model to control the phase behaviour of the system. A larger $\Omega$ results in a faster rotation of spinners, which leads to more energy transfer in active collisions. In Fig. 4B and C, we calculate $S(q)$ and $S_{R}(f)$ for the active spinners system at $\tilde{\rho}=0.15$ with different $\Omega$ (solid symbols) compared with the equilibrium dimer fluid (open symbols). When $\Omega$ is small, we find the same hyperuniform scaling $S(q\rightarrow 0)\sim q^{2}$ and $S_{R}(f\rightarrow 0)\sim f^{-1/2}$ in the spatial and temporal domains as previously found. The response of $S(q)$ and $S_{R}(f)$ to the increase of $\Omega$ is also similar to the response to $\Delta E$ in the previous model. Nevertheless, we emphasize that in the active spinner system the interacting force is not center to center and the energy $\Delta E$ is not a constant. These differences demonstrate the robustness of fluidic HU, but also request a modified hydrodynamic theory for active spinner fluids to describe the ‘odd viscosity’ and topological waves in the fluids (46, 47). Finally, we note that the suppressed density fluctuation was detected in a recent experimental air-driven spinner system (48) and we expect further experimental realization of hyperuniform spinner fluids using different techniques (49, 50, 51, 52, 53).

Discussion and Conclusion

In conclusion, a general hydrodynamic mechanism of fluidic HU is discovered by investigating a random-organizing hard-sphere fluid model and formulating a fluctuating hydrodynamic theory based on the generalized Navier-Stokes equations. We find that the state of the system is determined by two characteristic lengths: the dissipation length $l_{d}$, reflecting the strength of the particle inertia, and the mean free path $l_{m}$ associated with the particle density. A smooth transition from an equilibrium hard sphere fluid to a non-equilibrium hyperuniform fluid occurs with decreasing $l_{d}$ from infinity to $l_{m}$, below which the system falls into an absorbing state. We determine this absorbing phase transition belonging to the universality of conserved directed percolation (CDP) for systems with low particle inertia ($l_{d}$). Increasing the inertia induces the crossover of the universality class from CDP to its long-range mean-field scenario. These hyperuniform fluids have a kinetic temperature controlled by ratio $l_{d}/l_{m}$, which reflects how far the system is out of equilibrium, while $l_{d}$ alone sets the length scale, above which the hyperuniform scaling occurs. Importantly, we demonstrate that the mechanism of fluidic HU can be understood as the damping of a $q$-space stochastic harmonic oscillator, and the suppressed long-wavelength density fluctuation in the fluids exhibits as either acoustic (resonance) or diffusive (overdamped) modes. These dynamic modes are identified in the dynamic structure factor and the power spectrum of local density fluctuations. Our theory suggests two important ingredients for the random-organizing fluidic HU, i.e., (under/over)damped dynamics (dissipation) and reciprocal active interaction (driving). Accordingly, we show the fluidic HU with controllable hyperuniform length scales in an experimentally realizable active spinner system. Generally, the damping can be induced by the substrate friction, viscous solvent/gas or ‘optical molasses’ (laser cooling) (54), while the reciprocal active interaction can be realized through energy transfer between different degrees of freedom (51, 55), active deformation (56, 57), localized motion (24, 58), time-oscillatory potentials (59) or dielectric heating (60). Thus, we anticipate that other driven-dissipative systems studied previously can also produce the fluidic HU, e.g., the vibrated granular monolayer which allows the kinetic energy transferring from vertical to horizontal freedom degree (55), actively deformable micro-robots (56)/cells (57) on a substrate or in solvent, or molecules gas in optical molasses (61) with vibrational/rotational excitation (60, 62). These dynamic hyperuniform fluids, if realized, would be capable of self-healing from damage and self-adapting to the change of external driving. Thus, compared with solid hyperuniform materials (3, 9, 10, 11), they behave essentially like a living functional material (28), which may have potential applications in optics, acoustics and micro-fluidic engineering.

\matmethods\showmatmethods

In our simulations, we use square (2D) or cubic (3D) simulation boxes to calculate the $S(q)$ and $S_{R}(f)$, while in calculate $S(q,\omega)$ elongated boxes with aspect ratio $1:8$ (2D) and $1:1:8$ (3D) are used. The mean free path $l_{m}(\tilde{\rho})$ in the hard sphere model is calculated using randomly-generated static configurations to ensure that $l_{m}$ is only a function of $\tilde{\rho}$ and independent of the dynamic state of the system. We have further checked that HU has no significant influence on $l_{m}(\tilde{\rho})$. In the calculation of $l_{m}$, particles’ velocities are all set to zero. We then give a randomly-selected particle $i$ a random direction and obtain the maximum distance the particle can travel in a straight line (without damping) before colliding with another particle $j$. Then, the particle $i$ is put back to its original position and particle $j$ is selected as the next moving particle. The structure factor is estimated by the scattering intensity $S({\mathbf{q}})=\frac{1}{N}\left<\left|\sum_{j=0}^{N}e^{i{\mathbf{q}}\cdot{\mathbf{r}}_{j}}\right|^{2}\right>$ with ${\mathbf{q}}=[q_{x},q_{y},q_{z}]=[i\frac{2\pi}{L},j\frac{2\pi}{L},k\frac{2\pi}{L}]$ ( $i,j,k=1,2,3,\cdots$) and $L$ the box size. The vector $\mathbf{q}$ is then projected to the scalar $q$. For 2D systems, $k=0$. The dynamic structure factor is estimated as $S({\mathbf{q}},\omega)=\frac{1}{N\tau_{max}}\left|\sum_{j=0}^{N}\int_{0}^{\tau_{max}}e^{i[{\mathbf{q}}\cdot{\mathbf{r}}_{j}(t)-\omega t]}dt\right|^{2}$ with angular frequency $\omega=2\pi j/\tau_{max}$ ($j=1,2,3,\cdots$). The direction of $\mathbf{q}$ is chosen to be along with the elongated box. In the calculation of $S_{R}(f)$, the frequencies are chosen similarly and the sub-domain radius is set to $R=L/4$. To increase the efficiency, we sample 100 random-located sub-domains at the same time.

Details about the event-driven simulations, finite-size scaling analysis of absorbing phase transitions, hydrodynamic theory, theoretical derivations of $S_{R}(f)$ in 2D and 3D, the molecular dynamic simulation of active spinners are included in the SI Appendix Supplementary text. A table list of symbols used in this work can be found in SI Appendix Table S2. The original data for Fig. 1-4 is provided in SI Appendix Database S1-S9.

\acknow

The authors are grateful to Profs. Dov Levine and Hao Hu for helpful discussions. This work is supported by Nanyang Technological University Start-Up Grant (NTU-SUG: M4081781.120), the Academic Research Fund from Singapore Ministry of Education (M4011873.120), 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.

\showacknow