Emergent universal long-range structure in random-organizing systems

Domain: physics
Published: 2026-01-23
Editorial outcome: major_revision

Abstract
Self-organization through noisy interactions is ubiquitous across physics, mathematics, and machine learning, yet how long-range structure emerges from local noisy dynamics remains poorly understood. Here, we investigate three paradigmatic random-organizing particle systems drawn from distinct domains: models from soft matter physics (random organization, biased random organization) and machine learning (stochastic gradient descent), each characterized by distinct sources of noise. We discover universal long-range behavior across all systems, namely the suppression of long-range density fluctuations, governed solely by the noise correlation between particles. Furthermore, we establish a connection between the emergence of long-range structure and the tendency of stochastic gradient descent to favor flat regions of energy landscape—a phenomenon widely observed in machine learning. To rationalize these findings, we develop a fluctuating hydrodynamic theory that quantitatively captures all observations. Our study resolves long-standing questions about the microscopic origin of noise-induced hyperuniformity, uncovers striking parallels between stochastic gradient descent dynamics on particle system energy landscapes and neural network loss landscapes, and should have wide-ranging applications—from the self-assembly of hyperuniform materials to ecological population dynamics and the design of generalizable learning algorithms.

Introduction
While typically associated with disorder, noise can paradoxically drive the emergence of diverse forms of order, such as pattern formation1, self-organization2,3, suppression of chaos4, selection of ordered states5, and swarming6. Physical systems exhibit a broad spectrum of order: perfect crystals and ideal gases mark the extremes, while intermediate regimes display correlated disorder, such as hyperuniformity—where local disorder coexists with the anomalous suppression of long-range density fluctuations7. Hyperuniformity can emerge either at criticality8–10, or away from it11–14. Away from criticality, in equilibrium, hyperuniformity requires long-range interactions7, whereas out of equilibrium, it can emerge from long- or short-range, and even noisy interactions11–16. The process by which long-range spatial structure develops away from criticality—particularly in systems interacting solely via short-range, noisy dynamics—is a long-standing question that remains poorly understood.

Non-equilibrium particle systems with short-range noisy interactions—such as random organization (RO)2,9,17–20, biased random organization (BRO)10,12,14,19,21,22, and stochastic gradient descent (SGD)23—provide an ideal framework to investigate the noise-driven emergence of long-range spatial structure. These systems have been studied in a wide variety of contexts, such as sheared colloidal suspensions2,9,24,25, random close packing10,21, two-dimensional crystallization14, and self-supervised learning23,26. All systems undergo a phase transition as the particle volume fraction increases; from a low-density state where all motion ceases (absorbing state), to a high-density state where motion persists forever (active state)27–29. Irrespective of microscopic details, RO, BRO, and SGD belong to the same universality class, i.e., display the same critical behavior10,20,23. Away from criticality in the active phase, however, variations in microscopic interactions significantly influence the emergent long-range structure12,14,17,18. Despite extensive experimental, numerical, and theoretical research over two decades, a quantitative microscopic understanding of dynamics and structure far from criticality remains elusive2,9,10,12,14,17–19,21–24. Fundamental questions remain unanswered: How does macroscopic structure emerge from noisy interactions? Moreover, what universal principles govern the variability in emergent structures within and across different random-organizing systems? Finally, is the emergent long-range structure in SGD related to its ability to discover flat regions of energy landscape—a feature linked to robust generalization in machine learning30?

Here, combining particle and continuum simulations with hydrodynamic theory, we provide a quantitative, microscopic understanding of the active phases of RO, BRO, and SGD—random-organizing systems with distinct sources of noise. We discover universal long-range behavior across all three systems governed by a single parameter: the noise correlation coefficient between particles. All systems self-organize to suppress density fluctuations below a crossover length scale, which diverges as the noise becomes anti-correlated (reciprocal interactions), resulting in strong (class I7) hyperuniformity. Further, by directly coarse-graining the microscopic dynamics, we develop a fluctuating hydrodynamic theory for random-organizing systems that quantitatively predicts both the emergence of long-range structure and the crossover length scale across all systems. Finally, we demonstrate how noise correlation, batch size, and learning rate bias SGD towards flatter regions of energy landscape—a finding that aligns with empirical observations in machine learning30–33—and establish a connection between the emergence of long-range structure and the ability of SGD to generalize effectively. Our study underscore the critical role of noise correlations in facilitating long-range structure, and has wide-ranging applications—from providing a robust method for self-assembling hyperuniform structures, to offering insights into other systems having correlated noise, such as neural population activity in the brain34, ecological population dynamics35, and gene expression dynamics in cells36.

Results
Setup
Random-organizing systems—RO, BRO, and SGD—are discrete-time systems consisting of N interacting spherical particles of radius R in d-dimensional space. At any given time-step, the positions of isolated particles (blue in Fig. 1b–d) do not evolve, and those of overlapping (active) particles (red in Fig. 1b–d) evolve according to system specific rules designed to resolve particle overlaps. All systems undergo an absorbing-to-active phase transition as the particle volume fraction ϕ increases (Fig. 1a). Starting from a random initial configuration, for ϕ < ϕc, the system finds an absorbing configuration (fa = 0), whereas for ϕ > ϕc, the system never finds such a configuration and reaches a non-equilibrium steady state (fa > 0) (Fig. 1a). Here, fa is the fraction of active particles in the system. We report all our results in the active phase when the system has reached steady state (Methods).

Random organization
RO was originally introduced to model experiments on sheared colloidal suspensions at high Péclet number2,24. Consider a system undergoing periodic shear cycles. After every shear cycle, overlapping particles are given a “kick” in a random direction, and non-overlapping particles return to their original position2,9. RO was subsequently simplified to a model without external shearing, which retains all the essential properties of the original version17–19,29. Here, we study this simpler isotropic version of RO17–19,29.

In RO, the dynamics of the position of particle i at time-step m+1 (xm+1 i) is given by,
xm+1 i = xm i + ϵ ∑ j∈Γm i um ji ζm ji, (1)
where ϵ controls the magnitude of the pairwise kick given by particle j to i, um ji is a random number sampled from a standard uniform distribution (U[0,1]) at time-step m, ζm ji is a random unit vector sampled uniformly on the surface of a d-dimensional unit hypersphere at time-step m, and Γm i = {j | |xm j − xm i| < 2R, j ≠ i} is the set containing all particles that overlap with particle i at time-step m (Fig. 1b). We set ζm ji = −ζm ij so that particles effectively move away from (“repel”) each other. Finally, c ∈ [−1, 0] is the Pearson correlation coefficient between corresponding components of the complete noise vectors ωm ij,α = ϵum ji ζm ji,α and ωm ji,α, defined by ⟨ωm ij,α ωm ji,β⟩ = ⟨ωm ij,α⟩ ⟨ωm ji,β⟩ = cδ αβ. c = 0 denotes pairwise kicks which are uncorrelated, and c = −1 denotes pairwise kicks which are anti-correlated (equal in magnitude, um ij = um ji, and opposite in direction, i.e., conserve pairwise center of mass).

Biased random organization
BRO was originally introduced to study random close packing of spheres10,21. Similar to RO, overlapping particles are given a kick, however, the kick is “biased” along the direction of the line joining the centers of overlapping particles (Fig. 1c). An interesting feature of BRO is that for ϵ → 0, its critical point was proposed as an alternate definition of random close packing10,21, which, despite decades of research, still lacks a clear definition37–41.

In BRO, the dynamics of the position of particle i at time-step m+1 (xm+1 i) is given by,
xm+1 i = xm i + ϵ ∑ j∈Γm i um ji xbm ji, (2)
where ϵ and um ji are defined as in RO, and xbm ji = −(xm j − xm i)/|xm j − xm i| is the deterministic unit vector pointing from the center of particle j to i at time-step m (Fig. 1c). c ∈ [−1, 0] is the Pearson correlation coefficient between corresponding components of the complete noise vectors ωm ij,α = ϵum ji xbm ji,α and ωm ji,α. Notice that c = −1 corresponds to pairwise center of mass conservation.

Stochastic gradient descent
SGD is a widely used optimization algorithm, e.g., in artificial neural networks, to minimize a loss function composed of a sum of many terms42. Stochasticity in SGD comes from the random selection of a subset of terms in the sum at every step. In the context of interacting particles, E = ∑ i ∑ j≥i V(xi, xj), where V(xi, xj) is the pairwise potential. SGD corresponds to randomly selecting a subset of terms in E and updating the particle positions—in contrast to both simultaneously moving all active particles, which corresponds to (noiseless) gradient descent23. Notice that our approach of SGD applied on particle systems is unlike previous works relating neural networks to particle systems by treating parameters (weights) as interacting particles43.

In SGD, the dynamics of the position of particle i at time-step m+1 (xm+1 i) is given by,
xm+1 i = xm i − α ∑ j∈Γm i θm ji ∇ i Vm ji, (3)
where ∇ i Vm ji = ∇ xi V(xi, xj) is the pairwise interaction potential between particles i and j, and θm ji is a random number sampled from a Bernoulli distribution having parameter bf (batch fraction) at time-step m (Fig. 1d). bf represents the average fraction of active particle pairs (i, j) that move at any given time-step. c ∈ [−1, 0] is the Pearson correlation coefficient between corresponding components of the complete noise vectors ωm ij,α = −α θm ji ∂ i Vm ji,α. Notice that ∇ i Vm ji = −∇ j Vm ji is the pairwise interaction potential between particles i and j, and θm ij = θm ji. While V ji can be any short- or long-range potential in SGD, here, we consider a class of short-range, repulsive potentials given by
V ij(r) = E(1 − r ij/2R)p, if 0 < r ij < 2R, (4)
where r ij = |xj − xi|, E is the characteristic energy scale, and p controls the stiffness of the potential.

Universal active phase behavior
There are three distinct sources of noise in random-organizing systems: magnitude of kicks, direction of kicks, and selection of particles. Notice that the origin of stochasticity in RO, BRO, and SGD is different; (i) for RO, the magnitude and direction of kicks are both noisy, while the selection of particles is deterministic, (ii) for BRO, the magnitude of kicks is noisy, while both the direction of kicks and selection of particles is deterministic, (iii) for SGD, the selection of particles is noisy, while the magnitude and direction of kicks are both deterministic (Fig. 1b–d, Eqns. (1), (2), and (3)). We perform particle simulations for RO, BRO, and SGD in the active phase (ϕ > ϕc) and measure the long-range structure, quantified by the radially averaged structure factor S(k) and variance in number density δρ2(l) (Methods). l is the diameter of the hypershpere used for measuring density fluctuations and k = 2π/l is the wavenumber. We study all properties above a threshold length scale l0 = 2π/k0, below which we find short-range behavior (Fig. 2b, c). Here, we report our results for the long-range behavior (k < k0).

Despite having microscopically different dynamics, all random-organizing systems display a universal long-range structure, controlled solely by the pairwise noise correlation c, and qualitatively independent of all other parameters—be it the kick magnitude ϵ, volume fraction ϕ, or learning rate α. (Fig. 2a, b, c, Supplementary Information (SI) Figs. S2, S3, S4). Remarkably, the active phase behavior across all systems is independent of d, in contrast to their critical behavior, which is heavily dependent on d9,10,23 (Figs. S2c, S3c, S4c). All systems self-organize to suppress density fluctuations below a normalized crossover length scale lc; specifically, for length scales 1 < l < lc, the structure factor follows a power law (S(k) ∝ k2) and δρ2(l) ∝ l−(d+1), whereas for l > lc, the structure factor is constant, (S(k) ∝ const.) and δρ2(l) ∝ l−d (Fig. 2b, c). Further, the crossover length scale lc increases monotonically as c decreases: as the pairwise noise becomes more negatively correlated, density fluctuations are suppressed up to larger length scales (Fig. 2b inset). Consequently, the infinite wavelength density fluctuations, δρ2(l → ∞) ∝ S(k → 0), decrease monotonically as c decreases (Fig. 2c inset). Finally, when the noise is anti-correlated (c = −1), the crossover length scale lc → ∞ and the system becomes strongly hyperuniform, S(k → 0) ∝ k2, and δρ2(l → ∞) ∝ l−(d+1) (Fig. 2b, c).

Our results are consistent with previous studies on RO, which focused on the specific case of uncorrelated noise (c = 0)17,18, and on BRO, which focused on the specific case of anti-correlated noise (c = −1)12,14. SI Fig. S4, however, shows that all three systems—RO, BRO, and SGD—exhibit the same long-range behavior for any c ∈ [−1, 0].

Generalized model
We now develop a continuous-time model of discrete-time random-organizing systems. Using the framework of stochastic modified equations (SME)23,44, we approximate the discrete-time dynamics by a differential equation (SDE) (SI Sec. I.A). In the resulting generalized model, the dynamics of the position of particle i (xi) is given by an overdamped Langevin equation,
dxi(t) = −γ−1 ∑N j=1 ∇ i V ji + ∑N j=1 √Λ ji ξ ji(t), (5)
where γ is the friction constant, V ji, Λ ji are short-range, pairwise interaction potential and diffusion matrix between particles j and i, respectively, Λ ji denotes a matrix square root, and Eq. (5) is interpreted in the Itô sense. ξ ji is a pairwise, Gaussian noise given by the particle j to particle i having mean ⟨ξ ji,α(t)⟩ = 0 and covariance matrix ⟨ξ ij,α(t) ξ kl,β(t')⟩ = δ(t − t')δ αβ (δ ik δ jl + cδ il δ jk), where c ∈ [−1, 0] is the Pearson correlation coefficient between ξ ij,α(t) and ξ ji,α(t). Eq. (5), supplemented with system-specific γ, V ij and Λ ij, is an SDE approximating RO, BRO, and SGD (SI Sec. I.A, Table S1). Notice that in the generalized model, the source(s) of noise for RO and BRO remain the same, while the selection noise in SGD becomes a noise on the magnitude of the kicks (Fig. 1e–g, SI Sec. I.A)23.

We perform particle simulations of the generalized model for RO, BRO, and SGD and find that the long-range structure is quantitatively the same as that for their discrete-time counterparts (Methods, Fig. 2b, c). Thus, the generalized model serves as an accurate continuous-time approximation for all random-organizing systems.

Fluctuating hydrodynamic theory
Equipped with the generalized model, we formulate a theory for the evolution of the density field ρ(x, t). Dean’s method, originally introduced for Brownian particles with additive noise45, and later extended to study systems with multiplicative noise46–48, is a well-known approach for directly coarse-graining microscopic dynamics. However, it has not yet been extended to systems where noise is both pairwise and correlated across components. Starting from Eq. (5), we extend Dean’s method45 and its subsequent generalizations48 to incorporate pairwise correlated noise, and derive the resulting fluctuating hydrodynamic equation for ρ(x, t) in arbitrary spatial dimension d to get (SI Sec. I.B.2),
∂ρ(x,t)/∂t = −∇ · {ρ(x) v(x)} + ∇ ∇ : {D(x) ρ(x)} − ∇ · jn(x), (6)
where ∇ = ∇ x, and : denotes a double dot product. The velocity v(x) in Eq. (6) is given by
v(x) = −γ−1 ⟨∇V(x,y)⟩ ρ(y), (7)
where ⟨a⟩ ρ(y) = ∫ a ρ(y) dy, and V(x,y) is the “continuous” version of V ji, given by replacing xi and xj by x and y in Eq. (4). v(x) originates from the deterministic (first) term in Eq. (5), and can be understood as the average force at point x due to the local interaction potential, −⟨∇V(x,y)⟩ ρ(y), divided by the friction coefficient γ. The diffusion tensor D(x) in Eq. (6) is given by
D(x) = 1/2 ⟨Λ(x,y)⟩ ρ(y), (8)
D(x) originates from the noise (second) term in Eq. (5), and can be understood as the averaged diffusion tensor over the local density. The stochastic flux jn(x) in Eq. (6) is given by
jn(x) = −√ρ(x) √ρ(y) √Λ(x,y) η(x,y) dy, (9)
where √Λ denotes a matrix square root, and Λ(x, y) is the “continuous” version of Λ ji. η(x,y,t) is a vectorial, two-point, Gaussian noise field having mean ⟨η α(x, y, t)⟩ = 0, and covariance matrix ⟨η α(x,y,t) η β(u,w,t')⟩ = δ αβ δ(t − t') [δ(x − u)δ(y − w) + cδ(x − w)δ(y − u)], where c ∈ [−1, 0] is the Pearson correlation coefficient between η α(x, y, t) and η α(y, x, t). jn(x) originates from the noise (second) term in Eq. (5). Notice that the noise at any spatial location x can be viewed as the sum of independent kicks from ny ∝ ρ(y) δV particles at y given to nx ∝ ρ(x) δV particles at x, where δV is an infinitesimal volume. Then, since the sum of n Gaussian noises yields a Gaussian noise with standard deviation ∝ √n/ρ, we have jn ∝ √ρ(x) √ρ(y). Further, similar to Eq. (5), Λ(x,y) acts as a projection tensor, making the noise anisotropic. Finally, integrating over y collects contributions from different locations to the total stochastic flux at x.

We perform finite-difference simulations of Eq. (6) for RO, BRO, and SGD, and find that the long-range structure is quantitatively the same as particle simulations (Methods, Fig. 2b, c). Thus, our coarse-grained theory quantitatively captures the long-range structure for all random-organizing systems.

To gain further analytical insights on Eq. (6), we linearize ρ(x,t) to first order and derive the analytical form of the static structure factor (SI Sec. I.B.2),
S(k) = (1+c) + [M(1+c) − c] k2, (10)
where M is a known system-dependent constant (Eqs. S64, S66, and S68, SI Sec. I.B.2). Eq. (10) directly relates c, the microscopic noise correlation coefficient between a pair of particles, to the long-range structure in the system, and quantitatively predicts S(k) for all random-organizing systems without free parameters (Fig. 2b). Further, since δρ2(l → ∞) ∝ S(k → 0)7, Eq. (10) also predicts the behavior of infinite wavelength density fluctuations for all random-organizing systems (Fig. 2c inset).

Equation (10) shows a competition between two terms. The first term (1+c) makes the long-range structure random, while the second term [M(1+c) − c] k2 makes the long-range structure strongly hyperuniform. The ratio of these two competing terms gives the normalized crossover length scale,
lc = √M(1+c) / −c, (11)
which quantitatively predicts the crossover length scale for all random-organizing systems (Fig. 2b inset).

Two remarks are in order. First, in a variety of other systems having local center of mass conserving dynamics12,49–52 (equivalent to c = −1 in random-organizing systems) and displaying hyperuniformity in the active phase, the coarse-grained noise term is of Laplacian form ∇2[ρ(x)ω(x)], where ω(x) is a spatially uncorrelated, one-point Gaussian noise field. In contrast, the noise term in our microscopically derived theory (Eqs. (6), (9)) is fundamentally different from “Laplacian noise” in two ways: (i) it enters Eq. (6) in the “divergence” form and, (ii) is spatially correlated. Second, the widely used Fokker-Planck coarse-graining approach gives a density evolution equation identical to Eq. (6), but without the stochastic flux term (jn(x) = 0), which, after density linearization, trivially predicts spatially uniform steady-state density (ρ(x) = constant and S(k) = 0) independent of c (SI Sec. I.B.1). This demonstrates the crucial role of noise and noise correlations even at the coarse-grained scale in determining the long-range structure, thereby underscoring the importance of our theoretical framework.

Flatness of energy landscape in stochastic gradient descent
SGD is widely used for training neural networks, not only for its computational efficiency but also for its remarkable ability to steer neural networks toward flat regions of their loss landscapes32,42,53. This attribute is crucial for “good” learning algorithms since flatter regions are strongly correlated with better generalization performance on unseen data30–33. The bias towards flat regions, due to the selection noise in SGD, highlights the vital role of noise in shaping learning dynamics in neural networks. This, then, raises two questions: (i) Does SGD-driven descent of energy landscapes in particle systems similarly bias the dynamics toward flat regions, akin to neural networks? If so, (ii) Can this bias be linked to the long-range structure observed in these particle systems?

We first examine how the flatness of the explored regions of the energy landscape varies with noise correlation c, batch fraction bf, and learning rate α in SGD (Eq. (3)). We add a small noise N to the steady-state configurations X obtained in particle simulations to get a perturbed configuration X + ΔX and measure the change in the total energy of the system, ΔE = ⟨E(X+ΔX) − E(X)⟩ N ∝ Tr(H), where E = ∑ i ∑ j>i V ij, and H is the Hessian matrix of the system (Eq. (4), Fig. 3a, Methods)23. ΔE is a measure of the flatness of the energy landscape: lower ΔE corresponds to flatter regions (Fig. 3a)23,30. We find that ΔE decreases as c decreases, meaning that increase in long-range structure leads to flatter regions (Figs. 3b, 2b). We now fix c = −1 (anti-correlated noise), and find that ΔE increases with bf23 and decreases with α, suggesting that lower batch fractions and higher learning rates lead to flatter regions of energy landscape—consistent with results on SGD dynamics in neural networks (Fig. 3c, d)30–32. Notice that c = −1 and bf = 1 correspond to (noiseless) gradient descent (Fig. 3c). All results are qualitatively independent of particle volume fraction ϕ > ϕc, the potential V ij, and the spatial dimension d, useful features since typical neural manifold dimensions are O(100)26 (SI Fig. S5).

Finally, we investigate the relationship between ΔE and the long-range structure formed by SGD. Combining the change in S(k) before and after adding a small perturbation to the system54 with our fluctuating hydrodynamic theory (Eq. (10)), we derive an expression for ΔE(c, bf, α) (SI Sec. I.C, Eqs. S79, S81, and S82)—in quantitative agreement with numerical simulations without free parameters (Fig. 3b–d). Thus, the emergent long-range structure in SGD provides a quantitative framework for understanding the flatness of energy landscape, a key feature linked to generalization. This unites two seemingly disparate domains—neural networks and interacting-particle systems—by revealing that the bias of SGD towards flat regions is a universal hallmark of high-dimensional loss (or energy) landscapes, regardless of the underlying system. Beyond its relevance to machine learning algorithms, our framework may guide the design of self-organizing materials with tunable energetic and structural properties.

To understand the characteristics of the flatter regions of energy landscape explored by SGD, we also measure the energy E of the steady-state configurations X in particle simulations. The dependence of E on key parameters—noise correlation c, batch fraction bf, and learning rate α—is opposite to that of the energy change for small perturbations ΔE (SI Figs. S5, and S6). This finding is in quantitative agreement with predictions from our fluctuating hydrodynamic theory (SI. Sec. I.C, Eqs. S72, S74, S75, and Fig. S6). Notably, the magnitude of these changes differs significantly: variation in system parameters leads to a substantial change in ΔE ~ 50%, but only a minor change in E ~ 8%. Therefore, flatter regions of the energy landscape explored by SGD are associated with nearly constant or slightly higher energy. A similar trend is also observed in SGD dynamics on neural networks, where flatter regions of loss landscape often correlate with constant or slightly higher loss30.

Methods
Particle simulations
Our system consisted of N hyperspherical particles of radius R in a d-dimensional hypercubic box of side length L with periodic boundary conditions. The unit of length, time, and energy were chosen as 2R, δt, and E, respectively. δt is the simulation time-step for continuous-time simulations and E = 1 is the characteristic energy scale of the potential given by Eq. (4). N = 318309 and R = 1 were kept fixed and the particle volume fraction ϕ = NVs/Vc was varied by changing L. Vs and Vc denote volumes of a d-dimensional hypersphere of radius R, and hypercube of side length L, respectively. Particles were randomly distributed in the simulation box at t = 0. All simulations were run until the system reached a steady-state and all measurements were performed at steady-state. All results were averaged over 100 steady-state configurations.

Discrete-time simulations. For discrete-time simulations of RO, BRO, and SGD, the dynamics of particles were evolved according to Eqs. (1), (2), and (3), respectively. As evident from Eqs. (1), (2), and (3), positions of isolated particles (particles with no overlapping neighbors) do not evolve at any given time-step.

Dynamics for RO are controlled by four parameters, the kick magnitude ϵ, the particle volume fraction ϕ, the spatial dimension d, and the correlation of pairwise noise c. For the results reported in the main text, parameters were set as ϵ = 1, ϕ = 1.0, d = 2, and c ∈ [−1, 0]. The critical volume fraction ϕc ≈ 0.375 for this set of parameters.

Dynamics for BRO are controlled by four parameters, the kick magnitude ϵ, the particle volume fraction ϕ, the spatial dimension d, and the correlation of pairwise noise c. For the results reported in the main text, parameters were set as ϵ = 1, ϕ = 1.0, d = 2, and c ∈ [−1, 0]. The critical volume fraction ϕc ≈ 0.395 for this set of parameters.

Dynamics for SGD are controlled by six parameters, the learning rate α, the particle volume fraction ϕ, the batch fraction bf, the “stiffness” of the potential p, the spatial dimension d, and the correlation of pairwise noise c. For the results reported in the main text, parameters were set as α = 0.5, ϕ = 1.0, bf = 0.5, p = 1, d = 2, and c ∈ [−1, 0]. The critical volume fraction ϕc ≈ 0.615 for this set of parameters.

For the measurement of flatness of energy landscape in SGD, a configuration X ≡ {x1, x2, ..., xN} in the steady-state was perturbed by adding a Gaussian noise N to the positions of all particles, X + ΔX ≡ {x1 + N1, x2 + N2, ..., xN + NN}. All results were averaged over 5000 independent noise realizations. For the results reported in Fig. 3b, parameters were set as α = 0.5, bf = 0.5, ϕ = 1.0, p = 1, d = 2, σ = 0.01, and c ∈ [−1, 0]. For the results reported in Fig. 3c, parameters were set as α = 0.5, c = −1, ϕ = 1.0, p = 1, d = 2, σ = 0.01, and bf ∈ [0.1, 1.0]. For the results reported in Fig. 3d, parameters were set as bf = 0.5, c = −1, ϕ = 1.0, p = 1, d = 2, σ = 0.01, and α ∈ [0.5, 1.0].

Continuous-time simulations. For the continuous-time simulations of the generalized model of RO, BRO, and SGD, the dynamics of particles were evolved according to Eq. (5), supplemented with appropriate values of the friction constant γ, and the short-range interaction potential V ij and matrix Λ ij. The Euler-Maruyama method was used to solve Eq. (5), with a time-step δt = 1.0.

For the generalized model of RO, V ji = 0, and Λ ji,αβ = 1 (0,2R) (r ij)(ϵ2/dτ)δ αβ, where 1 (0, 2R) (r ij) is an indicator function such that 1 (0,2R) (r ij) = 1 ∀ r ij ∈ (0,2R) and 1 (0,2R) (r ij) = 0 otherwise. τ is the time scale quantifying the time elapsed in a discrete step. For the results reported in the main text, parameters were set as ϵ = 1, τ = 1.0, ϕ = 1.0, d = 2, and c ∈ [−1, 0].

For the generalized model of BRO, γ = Eτ = ϵR, Λ ji,αβ = (ϵ2R2/3τE2)∂ α V ji ∂ β V ji, and V ji is a short-range, pairwise, linear, repulsive potential (Eq. (4) with p = 1). For the results reported in the main text, parameters were set as ϵ = 1, τ = 1.0, ϕ = 1.0, d = 2, and c ∈ [−1, 0].

For the generalized model of SGD, γ = τ/αbf, Λ ji,αβ = (α2bf(1−bf)/τ)∂ α V ji ∂ β V ji, and V ji is given by Eq. (4). For the results reported in the main text, parameters were set as α = 1, bf = 0.5, τ = 1, p = 1, ϕ = 1.0, d = 2, and c ∈ [−1, 0].

Continuum simulations
Our system consisted of density ρ(x, t) evolving in a d-dimensional hypercubic box of side length L with periodic boundary conditions. The unit of length, time, and energy were chosen as R, δt, and E, respectively. δt is the simulation time-step, E = 1 is the characteristic energy scale of the continuous version of the potential given by Eq. (4), and 2R is the cutoff length of the potential. Space was discretized with a square grid of grid spacing δl = L/512. L = 200 and R = 1 were kept fixed and ∫ρ(x,t)dx = N was fixed at all times (density conservation). Density at all grid points was sampled from a Gaussian distribution with mean N/Vc, and standard deviation 0.01 N/Vc at t = 0. All other parameters for RO, BRO, and SGD were chosen to be the same as that of continuous-time particle simulations. All simulations were run until the system reached a steady-state and all measurements were performed at steady-state.

The finite-difference method combined with forward Euler time stepping was used to solve Eq. (6). Spatial derivatives for the stochastic flux and the diffusion term in Eq. (6) were approximated by the second-order central difference scheme. We split the velocity term (Eqs. (6), (7)) into two parts using chain rule,
γ−1 ∇ · {ρ(x) ∇⟨V(x,y)⟩ ρ(y)} = γ−1 ∇ρ(x) · ∇⟨V(x,y)⟩ ρ(y) + γ−1 ρ(x) ∇2⟨V(x,y)⟩ ρ(y) (12)
Spatial derivatives in the first and second term on the right hand side of Eq. (12) were approximated using the first-order forward difference, and the second-order central difference scheme, respectively. The finite-difference schemes were chosen to ensure strict density conservation and prevent checkerboard artifacts.

Structure factor and density fluctuations
The structure factor for a particle system is defined as S(k) = |ρb(k)|2/N, where ρb(k) = ∑N j=1 exp(−ik · xj), and ρ(x) = ∑N j=1 δ(x − xj). S(k) was calculated using the nonuniform fast Fourier transform68,69. The radial structure factor S(k) was then calculated by radially averaging S(k).

The structure factor for a continuous density field ρ(x) is defined as S(k) = |δρb(k)|2/⟨ρ⟩, where δρ(x) = ρ(x) − ⟨ρ⟩, and ⟨ρ⟩ = ∫ ρ(x) dx / ∫ dx. Number density variance δρ2(l) in a hyperspherical window of diameter l was measured by exploiting the exact relationship between S(k) and δρ2(l)7,
δρ2(l) = ⟨ρ⟩ 2d Γ(1+d/2) / (π1/2 l)d ∫ 1 0 k S(k) J d/2 (kl/2)2 dk, (13)
where ⟨ρ⟩ = N/Vc, Γ is the Gamma function, and J is the Bessel function of the first kind. The integral in Eq. (13) was evaluated numerically using Simpson’s rule, based on the measured S(k).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The authors study the suppression of long-range density fluctuations (i.e., hyperuniformity) in stochastic many-body\ndynamics. Specifically, they compare three dynamics: random organization (RO), biased random organization (BRO), and\nstochastic gradient descent (SGD) [Fig. 2]. They argue that these dynamics can be mapped into a generalized model [Eq.\n(5)]. Using analytical coarse-graining, they predict the emergence of a lengthscale [Eq.(11)] controlling the emergence of\nhyperfuniformity [Fig. 2], and examine the flatness of energy minimum [Fig. 3].\nThe manuscript is an interesting addition to the existing literature on hyperuniformity. The main novelty stems from exploiting\nanalogies between distinct stochastic dynamics that describe either soft matter systems or algorithms in machine learning.\nThe analytical and numerical results are clearly non-trivial and carefully presented. My main concern regards the time-\ndiscretization convention of the multiplicative noise in the stochastic dynamics [Eq.(5)], see details below.\nOverall, I recommend publication of the manuscript in Nature Communications, provided that the comments below are\nproperly addressed.\n1. The authors exploit analogies between SGD in particle systems and neural networks. Here, they should acknowledge\nprevious works that have already examined such analogies, for instance [Comm. Pure Appl. Math 75, 1889 (2022)]. In the\nderivation of the fluctuating hydrodynamics, the authors should refer to other studies which have addressed the coarse-\ngraining of stochastic dynamics with pairwise multiplicative noise; for instance [J Chem Phys 140, 234115 (2014)].\n2. The authors should clarify the time-discretization convention of the stochastic dynamics with multiplicative noise [Eq.(5)] to\nlift any ambiguity regarding the associated Fokker-Planck equation. In fact, the numerical simulations of Eq.(5) should\nproperly take into account this convention (e.g., Ito or Stratonovitch). Moreover, it would be useful to determine whether the\nsteady-state distribution is actually given by the Boltzmann factor associated with the pairwise potential V_ij [Eq.(4)].\n3. The authors should acknowledge that the stochastic dynamics of the empirical density [Eq.(6)] actually contains the same\ninformation as the many-body dynamics [Eq.(5)]. In other words, there is not any assumption (in particular, the two-point\ndensity is not assumed to be given by the product of empirical densities) in deriving the fluctuating hydrodynamic theory. It\nwould be useful to clarify whether the functional steady-state distribution of the empirical density is given by the Boltzmann\ndistribution associated with the functional free-energy (as expected for equilibrium dynamics).\n4. The authors should clarify in which regime of parameters the linearization of density fluctuations that underlies the\nanalytical form of the structure factor [Eq.(10)], should be valid. In fact, one expects that this regime corresponds to high\ndensities and weak interactions; for instance, see previous works where a similar linearization has been used [New J. Phys.\n16, 053032 (2014); Phys. Rev. X 9, 041026 (2019); Phys. Rev. Lett. 133, 268002 (2024); review: arXiv:2411.13467].", "reviewer_2": "The main contribution is this work is the analysis of a range of microscopic models using simulations, continuous theory and\ncoarse-graining that display in some part of their phase diagram a non-equilibrium steady state characterized by the fact that\nthe explored configurations are hyperuniform. We find the paper very nice and timely, but would like to author to take our\nremarks into account before recommending publication.\n1. Some details of the model definitions are confusing. We understand the control of the degree of anticorrelation in the\ndisplacement components but the author do not explain how in practice they achieve this: do they pick one component, and\nthen draw the second one as a sum of two terms weighted by c? Please clarify.\n2. A number of contributions could be cited much better: the work of Harukuni Ikeda who studied a range of models with\nspecial noise correlations in several publication should be acknowledged better. Also Ran Ni and corworkers and also Jack\net al. (Phys. Rev. Lett. 114, 060601 2015) have developed coarse-grained analysis for hyperuniformity, please comment on\npossible link.\n3. We are not convinced by the whole discussion of \"flat minima\". First, that part feels totally disconnected from the main\npoint of the work, and the questions raised appear somewhat artificial. Second, we are unsure about some details of the\nanalysis as the manuscript mentions that the small perturbations added to the system are in \"steady state\" configurations,\nwhich are not minima at all - can the authors specify more precisely what they do, minima of what they consider and how? In\naddition, the derived connection between S(k) and Delta E relies on many approximations that are not clear to us: linearized\nDean equation and low-k behavior of S(k) whereas energy and interactions rather stem from large k (interparticle distance). The\nagreement displayed is suspiciously good. Anyhow, the causality between S(k->0) and \"flat minima\" is largely\noveremphasized, as best a correlation holds, but given the above caveat, we are unsure. In summary: We recommend to\nremove that part entirely.\n4. In the same vein, the dilute Dean equation is, strictly speaking, valid for dilute or weakly interacting systems. Please\ncomment on this, and why would this approximation be valid for the present simulated systems is unclear (and indeed\nunreasonable).", "reviewer_3": "In the manuscript by Anand, et al., the authors investigate the long-range correlated structures in random organising models\n(RO and BRO) with noise of different correlations, and they also tried to generalise it to the stochastic gradient descent\n(SGD), which is a commonly used method for machine learning. The results are interesting, but I have a few questions that\nneed to be addressed before I could make any recommendation.\n1. The central idea in this work is the influence of c, which is the correlation between the noise vectors acting on two active\nparticles i and j. They found that the system is hyperuniform (S(k) ~ k^2) until a length-scale, which increases with\ndecreasing c and diverges at c = -1. Actually, this is effect can be probably also seen as a combination of a noise conserving\nthe center of mass and another (white) noise. Essentially, the effect studied here is probably just the effect of thermal noise\nor temperature on non-equilibrium hyperuniform fluids as in Ref [44]. In case the (white) noise is very strong, c = 0, which\ncorresponds to the high temperature case. In [44], it was found that the effect of thermal noise (of strength T) can be\ndecoupled from the noise conserving the center of mass of the system, and S(k) ~ k^2 up to a length-scale, which diverges at\nthe zero thermal noise limit (T -> 0). When T << 1, S(0) ~ T, which recovers Eq. 10 in this work. Therefore, I would suggest\nthe authors to compare their theory and results with the effect of thermal noise as in Ref [44]. If they are essentially the same,\nthe findings in this work may not be universal, as it may depend on the exact form, of which the noise is injected into the\nsystem, e.g., the colour and/or function form of (thermal) noise.\n2. Following the question above, in [44], it was found that S(0) ~ T is only valid at the low temperature region (T << 1), which\nis similar to the effect of thermal noise on equilibrium crystals at low temperature, and I think Eq. 10 may also have the\nsimilar constraint, which should be valid only when c is near -1. Therefore, I would suggest the authors to clarify this.\n3. I don't understand the different between RO and BRO as shown in Fig. 2, as my understanding is that both models are the\nsame at the fixed c, in which the noise vectors have the same correlation c between i and j. I would suggest the authors to\nclarify the difference between the two models, if there are any."}, "author_response": "Response to the reviewers’ report\nWe report below a detailed response to the reviewers’ reports. The original reviewer comments are in italic, and our\nresponse appears in normal font. When addressing the reviewers’ comments, we also explicitly list the changes we\nmade in the revised manuscript and highlight those changes in red here and in the revised manuscript.\nResponse to Reviewer 1\nWe thank the reviewer for their positive comment and their recommendation for our work to be published. Below we\naddress the questions raised by the reviewer.\n1. The authorsexploitanalogies betweenSGDin particlesystemsandneuralnetworks. Here, theyshould acknowledge\nprevious works that have already examined such analogies, for instance [Comm. Pure Appl. Math 75, 1889 (2022)].\nIn the derivation of the fluctuating hydrodynamics, the authors should refer to other studies which have addressed\nthe coarse-graining of stochastic dynamics with pairwise multiplicative noise; for instance [J Chem Phys 140, 234115\n(2014)].\nWe thank the reviewer for pointing us to the relevant literature. Following the suggestions of the reviewer, we\nhave added both references in the main text. The modified “Stochastic Gradient Descent” section of the main text\nincludingreference[Comm. PureAppl. Math75,1889(2022)]nowreads-“NoticethatourapproachofSGDapplied\non particle systems is unlike previous works relating neural networks to particle systems by treating parameters\n(weights) as interacting particles [1].”\nFurther, the modified “Fluctuating hydrodynamic theory” section of the main text including reference [J Chem\nPhys 140, 234115 (2014)] now reads - “Dean’s method, originally introduced for Brownian particles with additive\nnoise [2], and later extended to study systems with multiplicative noise [3–5], is a well-known approach for directly\ncoarse-graining microscopic dynamics.”\n2. The authors should clarify the time-discretization convention of the stochastic dynamics with multiplicative noise\n[Eq.(5)] to lift any ambiguity regarding the associated Fokker-Planck equation. In fact, the numerical simulations of\nEq.(5) should properly take into account this convention (e.g., Ito or Stratonovitch). Moreover, it would be useful to\ndetermine whether the steady-state distribution is actually given by the Boltzmann factor associated with the pairwise\npotential V [Eq.(4)].\nij\nWethankthereviewerforpointingoutaboutthetime-discretizationconvention. Eq.(5)inthemaintextisinterpreted\nin the Itˆo sense. We have added the phrase: “... and Eq. 5 is interpreted in the Itˆo sense”, after Eq.5 in the main\ntext. For the numerical simulations of Eq.5, we used the Euler–Maruyama method, which we had mentioned in the\nMethods section “Continuous-time simulations” of the main text.\nFurther,intheSI,weshowthatsincerandom-organizingsystemsareout-of-equilibrium,thesteady-statedistribution\nfor the N-body Fokker-Planck equation associated with Eq.(5) is not given by the Boltzmann distribution. We have\nadded the following discussion in SI. Sec. I.B.1. -\nSteady-state distribution. Here we address the question of whether the steady-state solution for the evolution of\ntheN-particledensityρ isgivenbytheBoltzmanndistribution,asforinteractingBrownianparticlesatequilibrium\nN\n[6, 7] (Eq. S16). Plugging the steady-state distribution ρs in Eq. S16 and noting that ∂ρs /∂t=0, we get,\nN N\n   (cid:34)(cid:32) (cid:33) (cid:35)\n1 (cid:88) (cid:88) 1(cid:88)(cid:88) (cid:88)\nγ\n∂ xi,α∂\nxi,α\nV jiρs N+\n2\n∂\nxi,α\n∂\nxj,β\nδ\nij\nΛ\nki,αβ\n(t)+cΛ\nij,αβ\nρs\nN\n=0. (R1)\ni j i j k\nNoticethatEq.R1isdifficulttosolveforarbitraryΛ andnoisecorrelationc. However,substitutingΛ =aδ ,\nij,αβ ij,αβ αβ\nand c=0 in Eq. R1, where a is some constant, we get after some algebra,\n(cid:20) (cid:21)\n(cid:88) 1 Na\n∂ (∂ U )ρs + ∂ ρs =0, (R2)\nxi,α γ xi,α i N 2 xi,α N\ni\n2\n(cid:80)\nwhere U = V , and N is the number of particles. The solution of Eq. R2 is given by the Boltzmann distribution\ni j ij\nρB = (1/Z)e−(2/Naγ)E, where E = (cid:80) (cid:80) V , and Z is the normalization factor. This is not surprising since\nN i j≥i ij\nsubstituting Λ = aδ and c = 0 in Eq. S13 reduces it to a system of interacting Brownian particles at equi-\nij,αβ αβ\nlibrium, whose steady-state solution is known to be the Boltzmann distribution [6, 7]. Notice, however, that since\nΛ ̸=aδ and c̸=0, in general, for random-organizing systems (Table S1), the Boltzmann distribution is not the\nij,αβ αβ\nsteady-state solution for Eq. R1.\n3. The authors should acknowledge that the stochastic dynamics of the empirical density [Eq.(6)] actually contains the\nsame information as the many-body dynamics [Eq.(5)]. In other words, there is not any assumption (in particular,\nthe two-point density is not assumed to be given by the product of empirical densities) in deriving the fluctuating\nhydrodynamic theory. It would be useful to clarify whether the functional steady-state distribution of the empirical\ndensity is given by the Boltzmann distribution associated with the functional free-energy (as expected for equilibrium\ndynamics).\nWe agree with the reviewer that the stochastic dynamics derived using Dean’s approach formally contains the same\ninformationasthemany-bodydynamics. Toreflectthis,wehaveremovedtheassumptionaboutthetwo-pointdensity\nbeing given by the product of empirical densities from the main text. The modified “Fluctuating hydrodynamic\ntheory” section of the main text now reads - “Starting from Eq. 5, we extend Dean’s method [2] to account for\nmultiplicative, and pairwise correlated noise, and derive the time evolution of ρ(x,t) in arbitrary dimension d to get\n...”\nWehave alsoadded in the SISec. I.B.2 - “Note thattaking the ensemble average(over noiserealizations) of Eq. S42,\nand using the mean-field approximation ρ (x ,x )=ρ (x )ρ (x ), we recover the density evolution equation derived\n2 1 2 1 1 1 2\npreviously using the Fokker-Planck method (Eq. S19, Sec. I.B.I)”\nFurther,weshowintheSIthatthefunctionalFokker-Planckequationassociatedwiththestochasticdensityevolution\n(Eq.6)doesnothavetheBoltzmanndensityfunctionalasitssolution,sincerandom-organizingsystemsareatanon-\nequilibrium steady-state. We have added the following discussion in SI. Sec. I.B.2. -\nFunctional steady-state distribution. The functional Fokker-Planck equation for Eq. S42 is given by,\n∂P[ρ,t] (cid:90) δ (cid:40) (cid:34) 1 (cid:18) (cid:90) (cid:19) 1 (cid:18) (cid:90) (cid:19)(cid:35)\n= dx − ∂ ρ(x) ρ(y)∂ V(x,y)dy + ∂ ∂ ρ(x) ρ(y)Λ (x,y)dy P\n∂t δρ(x) γ x,α x,α 2 x,α x,β αβ\n1(cid:90) δ (cid:32)(cid:20) (cid:18) (cid:90) (cid:19)(cid:21) (cid:33)(cid:41)\n+ dy ∂ ∂ δ(x−y)ρ(x) ρ(u)Λ(x,u) du+cρ(x)ρ(y)Λ(x,y) P , (R3)\n2 δρ(y) x,α y,µ αµ αµ\nwhere P[ρ,t] is the probability distribution functional of the density ρ(x,t). Plugging the steady-state distribution\nfunctional Ps in Eq. R3 and noting that ∂Ps/∂t=0, we get,\n(cid:90) δ (cid:40) (cid:34) 1 (cid:18) (cid:90) (cid:19) 1 (cid:18) (cid:90) (cid:19)(cid:35)\ndx − ∂ ρ(x) ρ(y)∂ V(x,y)dy + ∂ ∂ ρ(x) ρ(y)Λ (x,y)dy Ps\nδρ(x) γ x,α x,α 2 x,α x,β αβ\n1(cid:90) δ (cid:32)(cid:20) (cid:18) (cid:90) (cid:19)(cid:21) (cid:33)(cid:41)\n+ dy ∂ ∂ δ(x−y)ρ(x) ρ(u)Λ(x,u) du+cρ(x)ρ(y)Λ(x,y) Ps =0. (R4)\n2 δρ(y) x,α y,µ αµ αµ\nEq.R4ischallengingtosolveforgeneralΛ(x,y) andc,andnoclosed-formsolutionisknown. However,substituting\nαβ\nΛ(x,y) =aδ , and c=0 in Eq. R4, where a is some constant, we get,\nαβ αβ\n(cid:90) δ (cid:40) (cid:34) 1 (cid:18) (cid:90) (cid:19) aN (cid:35)\ndx − ∂ ρ(x) ρ(y)∂ V(x,y)dy + ∂ ∂ ρ(x) Ps\nδρ(x) γ x,α x,α 2 x,α x,α\n(cid:32) (cid:33)(cid:41)\naN (cid:90) δ\n+ dy [∂ ∂ (ρ(x)δ(x−y))]Ps =0, (R5)\n2 δρ(y) x,β y,β\n(cid:82)\nwhere ρ(x)dx = N. Using integration by parts, and the fact that ∂ f(x,y) = −∂ f(x,y) for any function f\nx,α y,α\n3\nwhich depends only on |x−y|, we get,\n(cid:90) δ (cid:40)(cid:34) 1 (cid:18) (cid:90) (cid:19) aN (cid:35) aN (cid:20) (cid:18) δPs (cid:19)(cid:21)(cid:41)\ndx ∂ ρ(x) ρ(y)∂ V(x,y)dy + ∂ ∂ ρ(x) Ps+ ∂ ρ(x)∂ =0 (R6)\nδρ(x) γ x,α x,α 2 x,α x,α 2 x,α x,α δρ(x)\nThe solution of Eq. R6 is given by the Boltzmann distribution functional Ps =PB[ρ]≡(1/Z)e−(1/Naγ)F[ρ], where Z\nis a constant, and F[ρ] is the free energy functional given by,\n(cid:90) (cid:90)(cid:90)\n(cid:0) (cid:1)\nF[ρ]=Naγ ρ(x) lnρ(x)−1 dx+ ρ(x)ρ(y)V(x,y)dxdy. (R7)\nThis can also be seen by the fact that substituting Λ(x,y) = aδ , and c = 0 in Eq. S42 reduces it to a standard\nαβ αβ\nDean’s equation for a system of interacting Brownian particles at equilibrium, whose steady-state solution is known\nto be the Boltzmann distribution functional [2, 8]. However, since Λ(x,y) ̸= aδ , and c ̸= 0, in general, for\nαβ αβ\nrandom-organizing systems (Table S1), the Boltzmann distribution functional is not the steady-state solution for\nEq. S42.\nFIG. R1: Long-range structure in the absorbing phase (ϕ<ϕ ) of random-organizing systems, where ϕ denotes the critical\nc c\nparticle volume fraction. Radially averaged structure factor S(k) versus dimensionless radial wave number k×(2R/2π) in\ndiscrete-time particle simulations, where R is the particle radius, for (a) random-organization (RO); (b) biased\nrandom-organization (BRO); and (c) stochastic gradient descent (SGD). All parameters are the same as in the main text (see\nMethods), except ϕ, which is set to ϕ/ϕ =0.9 for all systems.\nc\n4. The authors should clarify in which regime of parameters the linearization of density fluctuations that underlies\nthe analytical form of the structure factor [Eq.(10)], should be valid. In fact, one expects that this regime corresponds\nto high densities and weak interactions; for instance, see previous works where a similar linearization has been used\n[New J. Phys. 16, 053032 (2014); Phys. Rev. X 9, 041026 (2019); Phys. Rev. Lett. 133, 268002 (2024); review:\narXiv:2411.13467].\nWe agree with the reviewer that our linearized fluctuating hydrodynamic theory is valid for high densities and\nweak interactions, as is well known for linearized Dean’s-type hydrodynamics theories. We have performed particle\nsimulations in the absorbing phase (ϕ < ϕ ), which show that our theory is not valid in that regime (Fig. R1). We\nc\nhave also cited previous works which show the long-range structure at the critical point (ϕ=ϕ ) is also different that\nc\nwhat is predicted by out theory. We have now added the following paragraph in the “Discussion” section of the main\ntext -\nUnder what conditions is the linearization approximation of our fluctuating hydrodynamic theory valid (Eqs. 6, and\n10)? It is well known that such linearization approximations are valid for high density and soft overlap potentials\n[8–11]. Indeed, our linearized theory does not explain the long-range structure either in the absorbing phase (ϕ<ϕ )\nc\nor at criticality (ϕ = ϕ ). In the absorbing phase, the long-range structure remains unchanged by the noise correla-\nc\ntion c (SI Fig. S8). At criticality, the long-range structure is still independent of c but becomes strongly dependent\non dimension d: for d < 4 the systems are hyperuniform (S(k → 0) ∼ kα), whereas for d ≥ 4 they are random\n(S(k →0)∼const.), with the exponent α varying with d [12–18].\nWe thank the reviewer for their positive comments. Below we address the questions of the reviewer point-by-point.\n1. Some details of the model definitions are confusing. We understand the control of the degree of anticorrelation\nin the displacement components but the author do not explain how in practice they achieve this: do they pick one\ncomponent, and then draw the second one as a sum of two terms weighted by c? Please clarify.\nWe apologize for the confusion. The degree of anti-correlation in RO, BRO, and SGD is controlled solely by scalar\nrandom variables.\nInRO,arandomvectorischosenfromthesurfaceofaunit-hyperspherewhichdenotesthekickdirectionfromparticle\ni to j, and the kick direction from particle j to i is always the opposite of that. The kick magnitudes for a pair of\nparticles are drawn from a uniform distribution such that they are correlated with a correlation coefficient −c. Then,\nsince the direction of kicks is always opposite, the correlation coefficient of the complete noise vector is c.\nFor BRO, the kick direction is deterministic and always opposite for a pair of particles. The kick magnitude is\ncorrelated in a similar way as RO giving rise to the complete noise vector being correlated with Pearson correlation\ncoefficient c.\nFor SGD, the kick direction and magnitude are both deterministic. The selection, determined by a Bernoulli random\nvariable, is correlated for a pair of particles, and correlated Bernoulli random variables are drawn for each pair.\nTo clarify this, we have added the following line in the section “Biased Random Organization” in the main text -\n“Notice that c = −1 corresponds to um = um, since xˆm = −xˆm.”. Further, we have also added the following line in\nij ji ji ij\nthesection“StochasticGradientDescent”inthemaintext-Noticethatsince∇ Vm =−∇ Vm, c=−1corresponds\ni ji j ji\nto θm =θm.\nij ji\n2. A number of contributions could be cited much better: the work of Harukuni Ikeda who studied a range of models\nwith special noise correlations in several publication should be acknowledged better. Also Ran Ni and corworkers and\nalso Jack et al. (Phys. Rev. Lett. 114, 060601 2015) have developed coarse-grained analysis for hyperuniformity,\nplease comment on possible link.\nWe thank the reviewer for pointing us to additional literature. Following the suggestions of the reviewer, we have\ncited two relevant works of Harukuni Ikeda in the main text section “Fluctuating hydrodynamic theory” - “First, in\na variety of other systems having local center of mass conserving dynamics [19–23] (equivalent to c=−1 in random-\norganizing systems) and displaying hyperuniformity in the dense phase, the coarse-grained noise term is of Laplacian\n(cid:112)\nform ∇2[ ρ(x)ω(x)], where ω(x) is a spatially uncorrelated, one-point Gaussian noise field.”. Further, we have also\ncitedJacketal. (Phys. Rev. Lett. 114, 0606012015)inthe“Introduction”ofthemaintext-“Awayfromcriticality,\nin equilibrium, hyperuniformity requires long-range interactions [18], whereas out of equilibrium, it can emerge from\nlong- or short-range, and even noisy interactions [19, 24–28]. ”. We had already cited three relevant works of Ran\nNi and co-workers in the main text, which are, Lei et al., PNAS 116, 22983 (2019), Lei et al., J. Chem. Phys. 159\n(2023), and Lei et al., Science advances 5, eaau7423 (2019).\nWhile these studies all perform coarse-grained analyses of hyperuniformity, they differ from our work in both the\nmodels they consider and in their scope. Each examines a particular non-equilibrium system, uncovers hyperuni-\nform scaling numerically, and then develops a corresponding continuum theory, phenomenological or microscopically\nderived. Yet none tackles random-organizing systems or offers a theory—microscopic or phenomenological—that\nexplains the universal onset of hyperuniformity in the active phase across models. Moreover, they restrict themselves\nto strictly hyperuniform regimes and provide no tunable microscopic parameter (for example, a noise-correlation\ncoefficient c) that smoothly bridges random and hyperuniform behavior within a single framework.\n3. We are not convinced by the whole discussion of ”flat minima”. First, that part feels totally disconnected from\nthe main point of the work, and the questions raised appear somewhat artificial. Second, we are unsure about some\ndetailsoftheanalysisasthemanuscriptmentionsthatthesmallperturbationsaddedtothesystemarein”steadystate”\nconfigurations, which are not minima at all - can the authors specify more precisely what they do, minima of what they\nconsider and how? In addition, the derived connection between S(k) and Delta E relies on many approximations that\nare not clear to us: linearized Dean equation and low-k behavior of S(k) whereas energy and interactions rather stem\nfrom large k (interparticle distance). The agreement displayed is suspiciously good. Anyhow, the causality between\nunsure. In summary: We recommend to remove that part entirely.\nWe appreciate the reviewer’s careful reading and have revised the “flat minima” section to address the reviewer’s\nconcerns, but we have chosen to retain the section because, as we argue below, we believe it is relevant.\nOurworkfocusesontheuniversallong-rangestructureexhibitedbyrandom-organizingsystemsintheiractivephase.\nSince SGD is a unique random-organizing system that has a well-defined energy (in contrast to RO and BRO which\nhavenoexplicit energy, seeEqs.1, 2, and3inthemaintext), exploringtheconnectionbetweentheenergy-landscape\nperspective—unique to SGD—and the measured structure factor—common across all random organizing systems—is\nbothnaturalandinsightful. Moreover, sinceSGDiswidelyusedinneuralnetworks, whicharealsohigh-dimensional,\nit is worth asking whether the fact that SGD leads to flat regions of the landscape in neural networks also holds true\nfor particle systems, which are completely distinct from neural networks. Finally, since the structure and energetic\npropertiesareintertwinedforparticlesystems,itisrelevanttodiscusstheminthesamestudy,ratherthanseparately.\nThis strong connection is evidenced by the following exact relationship between energy E and S(k) [29, 30],\nρN (cid:90) N (cid:90)\nE = V(r)dr+ Vˆ(k)(S(k)−1)dk. (R8)\n2 2(2π)d\nwhere V(r)=V(|r −r |) is the pairwise potential, ρ is the number density, and N is the total number of particles.\nj i\nWeagreewiththereviewerthattheuseofterm“flatminima”(borroweddirectlyfromthedeeplearningliterature[31–\n33]) is misleading since the configurations found by SGD when the system relaxes to steady-state are not necessarily\nminima of the potential energy surface. Hence, we have replaced the word “flat minima” by “flat regions of energy\nlandscape” throughout the main text to remove any ambiguity.\nThe relationship between ∆E and S(k) derived in the SI Sec. I.C does indeed rely on two main assumptions: (i)\nlinearized Dean’s equation used to derive the S(k); (ii) the low-k behavior of S(k) instead of the full S(k) (Eq. 10\nin the main text). There are two reasons for this. First, for short-range potentials, Vˆ(k) is typically dominated\nby the small-k region; hence we can approximate S(k) in Eq. R8 by the small-k (long-range) behavior. Second,\nwhile the reviewer correctly points out that energy should be dictated by large-k behavior, we emphasize that since\n(cid:82)\nS(k)dk =1, the low- and high-k behavior of S(k) are not independent of each other. In other words, the change in\nlow-k behavior is also necessarily reflected in a change in the average number and/or overlap of particles. To further\nmake our argument concrete, we also measure the energy E of the steady-state configurations in particle simulations,\nandfindthatE canalsobequantitativelyrelatedtoS(k)throughourtheory(seeSI.Sec. I.Cforderivation,Eqs.S68,\nS71, S70, and Fig. R2). This reinforces the fact that the prediction of ∆E from theory is not a coincidence. In fact,\nother properties like the total energy E, can also be reliably predicted from our theory. We have now added the\nfollowing paragraph in the main text in the section “Flatness of energy landscape in stochastic gradient descent”:\nTo understand the characteristics of the flatter regions of energy landscape explored by SGD, we also measure the\nenergyE ofthesteady-stateconfigurationsXinparticlesimulations. ThedependenceofE onkeyparameters—noise\ncorrelation c, batch fraction b , and learning rate α—is opposite to that of the energy change for small perturbations\nf\n∆E (SI Figs. S5, and S6). This finding is in quantitative agreement with predictions from our fluctuating hydro-\ndynamic theory (SI. Sec. I.C, Eqs. S68, S71, S70, and Fig. S6). Notably, the magnitude of these changes differs\nsignificantly: variation in system parameters leads to a substantial change in ∆E ∼50%, but only a minor change in\nE ∼ 8%. Therefore, flatter regions of the energy landscape explored by SGD are associated with nearly constant or\nslightly higher energy. A similar trend is also observed in SGD dynamics on neural networks, where flatter regions of\nloss landscape often correlate with constant or slightly higher loss [31].\n4. In the same vein, the dilute Dean equation is, strictly speaking, valid for dilute or weakly interacting systems.\nPlease comment on this, and why would this approximation be valid for the present simulated systems is unclear (and\nindeed unreasonable).\nOur linearized fluctuating hydrodynamic theory is valid for high densities and weak interactions, as is well known\nfor linearized Dean’s-type hydrodynamics theories [8–11]. We have performed particle simulations in the absorbing\nphase (ϕ<ϕ ), which show that our theory is not valid in that regime (Fig. R1). We have also cited previous works\nc\nwhich show that the long-range structure at the critical point (ϕ=ϕ ) is also different than what is predicted by our\nc\ntheory. We have now added the following paragraph in the “Discussion” section of the main text:\nUnder what conditions is the linearization approximation of our fluctuating hydrodynamic theory valid (Eqs. 6, and\n(SGD). Normalized energy E˜(c) versus noise correlation c for different particle volume fraction ϕ (a), spatial dimension d (b),\nand stiffness of the potential p (Eq. S11) (c). E(c) is normalized as: E˜(c)=[E(c)−E(c=0)]/[E(c=−1)−E(c=0)]. Black\nlines show prediction of Eq. S68 (SI Sec. I.C). Normalized energy E˜(b ) versus batch fraction b for different particle volume\nf f\nfraction ϕ (d), spatial dimension d (e), and stiffness of the potential p (Eq. S11) (f). E(b ) is normalized as:\nf\nE˜(b )=[E(b )−E(b =1.0)]/[E(b =0.1)−E(b =1.0)]. Black lines show prediction of Eq. S71 (SI Sec. I.C). Normalized\nf f f f f\nenergy change E˜(α) versus learning rate α for different particle volume fraction ϕ (a), spatial dimension d (b), and stiffness of\nthe potential p (Eq. S11) (c). E(α) is normalized as: E˜(α)=[E(α)−E(α=0.5)]/[E(α=1.0)−E(α=0.5)]. Black lines\nshow prediction of Eq. S70 (SI Sec. I.C). All parameters are the same as Fig. S5 (see caption). All symbols denote\ndiscrete-time particle simulations.\n10)? It is well known that such linearization approximations are valid for high density and soft overlap potentials\n[8–11]. Indeed, our linearized theory does not explain the long-range structure either in the absorbing phase (ϕ<ϕ )\nc\nor at criticality (ϕ = ϕ ). In the absorbing phase, the long-range structure remains unchanged by the noise correla-\nc\ntion c (SI Fig. S8). At criticality, the long-range structure is still independent of c but becomes strongly dependent\non dimension d: for d < 4 the systems are hyperuniform (S(k → 0) ∼ kα), whereas for d ≥ 4 they are random\n(S(k →0)∼const.), with the exponent α varying with d [12–18].\nResponse to Reviewer 3\nWe thank the reviewer for their comments. Below we address the questions of the reviewer point-by-point.\n1. The central idea in this work is the influence of c, which is the correlation between the noise vectors acting on two\nactive particles i and j. They found that the system is hyperuniform (S(k)∼k2) until a length-scale, which increases\nwith decreasing c and diverges at c=−1. Actually, this is effect can be probably also seen as a combination of a noise\nstructure factor S˜(k˜) versus normalized radial wave number k˜ for random-organization (RO) (a), biased random-organization\n(BRO) (b), and stochastic gradient descent (SGD) (c) in discrete-time particle simulations. All parameters are the same as\nFig.2inmaintext(seeMethods). S˜=S(k)/S (2π/L)whereS (2π/L)isthestructurefactorforc=0atk=2π/L,andLis\n0 0\nthe side length of the simulation box. k˜=k/k where k is the value at which S˜(k )=1 for anti-correlated noise (c=−1) of\n0 0 0\nthe same system. The solid colored lines show a combined best fit of Eq. S81 to various noise correlations c for each system,\nusing Dth as the single fitting parameter. Insets show Sˆ(k˜→0) versus thermal diffusion coefficient Dth.\nSˆ(k˜→0,c,Dth)=S˜(k˜→0,c,Dth)/S˜(k˜→0,c=0,Dth =0). Solid black lines in the insets denote 1+c+2Dth/ρ¯A\n1\n(prediction of Eq. S81 in the limit k˜→0). Gray shaded regions denote short-range behavior (k˜>1).\nconserving the center of mass and another (white) noise. Essentially, the effect studied here is probably just the effect\nof thermal noise or temperature on non-equilibrium hyperuniform fluids as in Ref [44]. In case the (white) noise is\nvery strong, c = 0, which corresponds to the high temperature case. In [44], it was found that the effect of thermal\nnoise (of strength T) can be decoupled from the noise conserving the center of mass of the system, and S(k)∼k2 up\nto a length-scale, which diverges at the zero thermal noise limit (T ∼0). When T <<1, S(0)∼T, which recovers Eq.\n10 in this work. Therefore, I would suggest the authors to compare their theory and results with the effect of thermal\nnoise as in Ref [44]. If they are essentially the same, the findings in this work may not be universal, as it may depend\non the exact form, of which the noise is injected into the system, e.g., the colour and/or function form of (thermal)\nnoise.\nFollowing the suggestions of the reviewer, we investigated the effect of thermal noise on random-organizing systems,\nwhich are otherwise inherently athermal. We find, both analytically and in particle simulations, that, similar to ref.\n[20], theadditionofthermalnoisechangestheS(k →0)behaviorinallsystems, withaconstantfactorthatincreases\nwith temperature T, in particular, S(k → 0) ∝ 1+c+aT, where a is a system-dependent constant (see SI Sec.I.D.\nfor derivation, Eq. S81). This is similar to how hyperunformity is destroyed in ref. [20], and also in equilibrium\ncrystals. In our work, we propose a more general and fundamental mechanism. We do not have two competing\nnoise sources. Instead, we have a single driving noise source whose intrinsic nature is controlled by the correlation\nparameter c, which directly tunes the degree of correlation in the pairwise interactions, from anti-correlated (c=−1)\ntocompletelyuncorrelated(c=0). Thisisincontrasttosystemshavingacenter-of-massconservingnoiseinaddition\nto a thermal bath, where two independent noise sources act on the system [20]. Adding a thermal noise to random-\norganizingsystems,however,increasestheinfinite-wavelengthdensityfluctuationsinanidenticalmannerforallnoise\n8\ncorrelations c. We have summarized this in the “Discussion” section of the main text -\nRandom-organizingsystemsareinherentlyathermal(Eqs.1,2,and3). Weinvestigatehowtheirlong-rangestructure\n√\nchangesatafinitetemperaturebyincorporatingthermalnoise, 2Dthξth(t),intothegeneralizedmodel(Eq.5),where\ni\nξth(t) is a standard Gaussian white noise, Dth = k T/γ, k is the Boltzmann constant, and T is the temperature.\ni B B\nWe assume that the thermal noise and the intrinsic athermal noise, ξ (t), are uncorrelated (Eqs. 5, and SI S79).\nji\nWe find, both analytically and in particle simulations, that the infinite wavelength density fluctuations, quantified\nby S(k → 0), increase with temperature, across all noise correlations c and systems (SI Sec. I.D, Eq. S81, Fig. S7).\nNotably, systems with anti-correlated noise (c = −1), which are hyperuniform when athermal, show S(k → 0) ∝ T,\nillustrating how hyperuniformity is affected by thermal fluctuations (SI Fig. S7 inset). This behavior, also observed\nin other non-equilibrium hyperuniform systems [20], is similar to how thermal fluctuations weaken hyperuniformity\nin equilibrium crystals, as described by the fluctuation-compressibility relationship S(k = 0) = κρk T, where κ is\nB\nthe isothermal compressibility [18, 30].\nMoreover, even the case of uncorrelated noise (c=0) in random-organizing systems cannot be mapped to a thermal-\ntype additive noise. This is because the athermal noise in random-organizing systems is multiplicative, which cannot\nbe simply replaced by an additive thermal noise even for c=0 (see Eq. 5 in the main text). Indeed, when we choose\nc=0and Λ =aδ toremovethemultiplicativenatureofthenoise, Eq.5canbeeffectivelyreducedtoasystem\nij,αβ αβ\nof interacting Brownian particles in a thermal bath at equilibrium (see the discussion on “Steady-state distribution”\nin SI Sec.I. B.1). However, Λ(x,y) ̸=aδ and c̸=0, in general, for random-organizing systems (see Table S1).\nαβ αβ\nFinally, the universal active phase behavior in random-organizing systems stems from the fact that while the indi-\nvidual noise types and sources are different in all systems, they can nevertheless be approximated by a generalized\nmodel(seeEq.5ofthemaintext, andSI.Sec. 1.Afortheprocedureofapproximatingthediscrete-timedynamicsby\nacontinuous-timeSDE).Hence, ourresultsareagnostictothespecifictypeorformofnoiseintheoriginaldynamics,\nas long as it can be mapped to the generalized model (Eq. 5).\n2. Following the question above, in [44], it was found that S(0) ∼ T is only valid at the low temperature region\n(T <<1), which is similar to the effect of thermal noise on equilibrium crystals at low temperature, and I think Eq.\n10 may also have the similar constraint, which should be valid only when c is near -1. Therefore, I would suggest\nthe authors to clarify this.\nWe stress that Eq. 10 is valid for all noise correlations c ∈ [−1,0] when the system is athermal, as confirmed by\nthe excellent agreement between our theory and particle and continuum simulations (Fig. 2). Upon the addition of\nthermal noise, we have shown that S(k →0)∝1+c+aT, hence, S(k →0)∝T is valid only for c=−1 (Fig. R3).\n3. I don’t understand the different between RO and BRO as shown in Fig. 2, as my understanding is that both models\nare the same at the fixed c, in which the noise vectors have the same correlation c between i and j. I would suggest\nthe authors to clarify the difference between the two models, if there are any.\nRO has two sources of noise - (i) direction of kicks, and (ii) magnitude of kicks. BRO has a single source of noise - (i)\nmagnitude of kicks. At a fixed c, the direction of the noise vector will be opposite in both systems. But the direction\nofthenoisevectorwillalwaysbealongthelinejoiningthecenterofparticlesinBRO,whereasinRO,itcanbeinany\nrandomly chosen direction. We had already pointed this out in the main text in the section “Universal active phase\nbehavior” - “There are three distinct sources of noise in random-organizing systems: magnitude of kicks, direction of\nkicks, and selection of particles. Notice that the origin of stochasticity in RO, BRO, and SGD is different; (i) for RO,\nthe magnitude and direction of kicks are both noisy, while the selection of particles is deterministic, (ii) for BRO,\nthe magnitude of kicks is noisy, while both the direction of kicks and selection of particles is deterministic, (iii) for\nSGD, the selection of particles is noisy, while the magnitude and direction of kicks are both deterministic (Figs. 1b-d,\nEqns. 2, 3, and 4).”\n[1] G. Rotskoff and E. Vanden-Eijnden, Communications on Pure and Applied Mathematics 75, 1889 (2022).\n[2] D. S. Dean, Journal of Physics A: Mathematical and General 29, L613 (1996).\n[3] E. Bertin, H. Chat´e, F. Ginelli, S. Mishra, A. Peshkov, and S. Ramaswamy, New journal of physics 15, 085032 (2013).\n9\n[4] A. P. Solon, M. E. Cates, and J. Tailleur, The European Physical Journal Special Topics 224, 1231 (2015).\n[5] A. Donev and E. Vanden-Eijnden, The Journal of chemical physics 140 (2014).\n[6] P.-H. Chavanis, Physica A: Statistical Mechanics and its Applications 390, 1546 (2011).\n[7] H. Risken, The fokker-planck equation (1996).\n[8] P. Illien, arXiv preprint arXiv:2411.13467 (2024).\n[9] P. Illien and A. Carof, arXiv preprint arXiv:2501.16206 (2025).\n[10] V. D´emery, O. B´enichou, and H. Jacquin, New Journal of Physics 16, 053032 (2014).\n[11] P.-H. Chavanis, Physica A: Statistical Mechanics and its Applications 387, 5716 (2008).\n[12] S. Wilken, A. Z. Guo, D. Levine, and P. M. Chaikin, Physical review letters 131, 238202 (2023).\n[13] G. Zhang and S. Martiniani, arXiv preprint arXiv:2411.11834 (2024).\n[14] E. Tjhung and L. Berthier, Physical review letters 114, 148301 (2015).\n[15] E. Tjhung and L. Berthier, Journal of Statistical Mechanics: Theory and Experiment 2016, 033501 (2016).\n[16] K. J. Wiese, Physical Review Letters 133, 067103 (2024).\n[17] X. Ma, J. Pausch, G. Pruessner, and M. E. Cates, arXiv preprint arXiv:2507.07793 (2025).\n[18] S. Torquato, Physics Reports 745, 1 (2018).\n[19] D. Hexner and D. Levine, Physical review letters 118, 020601 (2017).\n[20] Q.-L. Lei, M. P. Ciamarra, and R. Ni, Science advances 5, eaau7423 (2019).\n[21] F. De Luca, X. Ma, C. Nardini, and M. E. Cates, Journal of Physics: Condensed Matter 36, 405101 (2024).\n[22] H. Ikeda, Physical Review E 108, 064119 (2023).\n[23] H. Ikeda, SciPost Physics 17, 103 (2024).\n[24] Q.-L. Lei and R. Ni, Proceedings of the National Academy of Sciences 116, 22983 (2019).\n[25] M.Huang,W.Hu,S.Yang,Q.-X.Liu,andH.Zhang,ProceedingsoftheNationalAcademyofSciences118,e2100493118\n(2021).\n[26] L. Galliano, M. E. Cates, and L. Berthier, Physical Review Letters 131, 047101 (2023).\n[27] R. L. Jack, I. R. Thompson, and P. Sollich, Physical review letters 114, 060601 (2015).\n[28] Y. Kuroda and K. Miyazaki, Journal of Statistical Mechanics: Theory and Experiment 2023, 103203 (2023).\n[29] S. Torquato, G. Zhang, and F. H. Stillinger, Physical Review X 5, 021020 (2015).\n[30] J.-P. Hansen and I. R. McDonald, Theory of simple liquids: with applications to soft matter (Academic press, 2013).\n[31] S. Jastrzebski, Z. Kenton, D. Arpit, N. Ballas, A. Fischer, Y. Bengio, and A. Storkey, arXiv preprint arXiv:1711.04623\n(2017).\n[32] N. S. Keskar, D. Mudigere, J. Nocedal, M. Smelyanskiy, and P. T. P. Tang, arXiv preprint arXiv:1609.04836 (2016).\n[33] Z. Xie, I. Sato, and M. Sugiyama, arXiv preprint arXiv:2002.03495 (2020)."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have significantly changed the main text and the supplementary information to address the issues raised by the\nreferees. I still have some comments, as detailed below. Should these comments be properly addressed by the authors,\nI\nwould recommend the manuscript for publication at Nature Communications.\n1. In section 'Fluctuating hydrodynamic theory' of the main text, the authors discuss some previous works on coarse-graining\nstochastic dynamics with multiplicative noise, and mention\n'However, it has not yet been extended to systems with pairwise correlated noise between components'\nAs far as I understand, Ref.[48] actually addresses the case of Brownian particles with a multiplicative pairwise noise. It\nwould be useful if the authors clarified to which extent their derivation actually differs from that of Ref.[48]\n2. In the supplementary information, the authors demonstrate that their dynamics follows a Boltzmann steady-state\nmeasure at equilibrium, as it should. Here, equilibrium corresponds to choosing a specific noise statistics that amounts to enforcing\nsome additive, additive, diagonal noise correlations. In other words, the three models considered by the authors are generally out-of-\nequilibrium. Although 'Random Organization' and 'Biased Random Organization' are clearly inspired by nonequilibrium\nsystems, it could be useful to comment further on why 'Stochastic Gradient Descent' should be regarded as a nonequilibrium\ndynamics.\nIn section 'Stochastic Gradient Descent' of the main text, the authors mention that SGD is meant to minimize the total energy.\nIt would be insightful to show that the steady-state distribution of the stochastic dynamics in Eq.(5) is indeed highest at the\nminimum energy of the system. In other words, the authors should demonstrate that such a property holds for the specific\nchoice of diffusion matrix in table S1.", "reviewer_2": "I have reviewed the response and revisions of the author to the comments made in the first round and I am satisfied by their\nwork. I can recommend publication at this point.", "reviewer_3": "I appreciate the efforts that the authors have spent in revising the manuscript and the explanation they have offered to the\nquestions in my previous report. I agree with the authors that the influence of c is not exactly the same as the effect of thermal\nnoise in random organising hyperuniform fluids in the literature, as they come with slightly different mathematical forms. But\nI think the essential physics in the effect of c is the same as adding an extra noise, and in this work it is a white noise but not\nexactly the same way as being added in the active hyperuniform fluids previously. I believe that the effect of c can be also\ndecoupled as an extra noise term in addition to the noise conserving the center of the mass of the system. Then the finding\nin this work is not \"universal\" as claimed by the authors, as the effect should depend on the exact form of the extra noise,\nsuch as the color of the noise, the mathematical form of correlation in the noise, etc. Therefore, I don't think this manuscript is\nsuitable for publications in high profile journals like Nature Communications, but rather suggest the authors to resubmit it\n(with revision) to more specialised journals."}, "author_response": "Response to the reviewers’ report\nWe report below a detailed response to the reviewers’ reports. The original reviewer comments are in italic, and our\nresponse appears in normal font. When addressing the reviewers’ comments, we also explicitly list the changes we\nmade in the revised manuscript and highlight those changes in red here and in the revised manuscript.\nResponse to Reviewer 1\nWe thank the reviewer for their positive comments and for recommending publication of the work. Below we address\nthe additional questions raised by the reviewer.\n1. In section “Fluctuating hydrodynamic theory” of the main text, the authors discuss some previous works on coarse-\ngraining stochastic dynamics with multiplicative noise, and mention ‘However, it has not yet been extended to systems\nwith pairwise correlated noise between components’. As far as I understand, Ref.[48] actually addresses the case of\nBrownian particles with a multiplicative pairwise noise. It would be useful if the authors clarified to which extent their\nderivation actually differs from that of Ref.[48].\nWe agree with the reviewer that the ref.[48] (J Chem Phys 140, 234115 (2014)) studies the case of Brownian particles\nwithamultiplicative,pairwisenoise. However,thenoise,B inref.[48]actsindependentlyoneachparticle(⟨B B ⟩∼\ni i j\nδ ), whereas in our work, the noise term ζ is explicitly pairwise and correlated between components (⟨ζ ζ ⟩ ∼\nij ij ij kl\nδ δ + cδ δ ). This is the main feature of our work, and as we show, the correlated noise between particles,\nik jl il jk\ncontrolled by the noise correlation coefficient c, is essential for the emergence of long-range structure.\nFollowing the suggestions of the reviewer, we have clarified the distinction between “pairwise multiplicative” and\n“pairwise multiplicative and correlated” noise, and the modified “Fluctuating hydrodynamic theory” section of the\nmain text now reads - “However, it has not yet been extended to systems where noise is both pairwise and corre-\nlated across components. Starting from Eq. 5, we extend Dean’s method [1] and its subsequent generalizations [48]\nto incorporate pairwise correlated noise, and derive the resulting fluctuating hydrodynamic equation for ρ(x,t) in\narbitrary spatial dimension d (see SI Sec. I.B.2).”\n2. In the supplementary information, the authors demonstrate that their dynamics follows a Boltzmann steady-state\nmeasure at equilibrium, as it should. Here, equilibrium corresponds to choosing a specific noise statistics that amounts\nto enforcing some additive, diagonal noise correlations. In other words, the three models considered by the authors are\ngenerally out-of-equilibrium. Although ‘Random Organization’ and ‘Biased Random Organization’ are clearly inspired\nby nonequilibrium systems, it could be useful to comment further on why ‘Stochastic Gradient Descent’ should be\nregarded as a nonequilibrium dynamics.\nNoise in Stochastic Gradient Descent (SGD) comes from the selection of particle pairs. Since at thermal equilibrium,\nthereisno“selection”ofparticles,theSGDnoiseisinherentlynon-equilibriuminorigin. Asimplewaytounderstand\nthis is to look at the generalized model of random-organizing systems (Eq. 5), which approximates the dynamics of\nSGD. The noise term in Eq. 5 is multiplicative and correlated between particles—both of which are not the features\nof a stochastic system at thermal equilibrium. From an energetic perspective, energy is continuously added to and\nremoved from the system as a consequence of batch selection in SGD.\nFollowing the suggestions of the reviewer, we have added this explanation in the SI, and the modified SI Sec.I.A now\nreads - “For an overdamped system of interacting particles at thermal equilibrium, the noise is additive, and uncor-\nrelated across particles. In contrast, the noise in Eq. S13 is multiplicative and pairwise correlated between particles.\nThis distinction provides a clear microscopic explanation for why RO, BRO, and SGD dynamics are intrinsically out\nof equilibrium.”\n3. In section ‘Stochastic Gradient Descent’ of the main text, the authors mention that SGD is meant to minimize\nthe total energy. It would be insightful to show that the steady-state distribution of the stochastic dynamics in Eq.(5)\nis indeed highest at the minimum energy of the system. In other words, the authors should demonstrate that such a\nproperty holds for the specific choice of diffusion matrix in table S1.\nis different from gradient flow because of (i) a finite learning rate, and (ii) batch selection [2]. Thus, SGD does not\nminimize the total energy, rather the partial energy, as mentioned in “Stochastic Gradient Descent” section of the\nmain text, “SGD then corresponds to randomly selecting a subset of terms in E and updating the corresponding\nparticle positions—either one or both at once—to minimize the partial energy.”\nWe now show that SGD does not strictly minimize the total energy. The total energy E is,\n1(cid:88)(cid:88)\nE = V(x ,x ). (R1)\n2 i j\ni j̸=i\nThe variation in E can then be written as (up to second order in ∆x ),\ni\n1(cid:88)(cid:88)\ndE = [V(x +∆x ,x +∆x )−V(x ,x )]\n2 i i j j i j\ni j̸=i\n(cid:20) (cid:21)\n1(cid:88)(cid:88) 1 1\n= ∆xT∇ V +∆xT∇ V + ∆xT∇ ∇ V ∆x + ∆xT∇ ∇ V ∆x +∆xT∇ ∇ V ∆x .\n2 i xi ij j xj ij 2 i xi xi ij i 2 j xj xj ij j i xi xj ij j\ni j̸=i\n(R2)\nWe know that (cid:10) ζ (t)ζT(t′) (cid:11) =δ(t−t′)[δ δ +cδ δ ]I, and ∆x =x (t+∆t)−x (t) evolves according to the SDE\nij kl ik jl il jk i i i\nd d x t i =− γ 1 (cid:80) j̸=i ∇ xi V(x i ,x j )+ (cid:80) j̸=i (cid:112) Λ ji ζ ji . WecanthenwritethefollowingequationsundertheItointerpretation,\n1 (cid:88)\n⟨∆x ⟩=− ∇ V ∆t\ni γ xi ij\nj̸=i\n(cid:10)\n∆x\n∆xT(cid:11)\n=δ\n(cid:88)\nΛ ∆t+cΛ ∆t, (R3)\nik jk ij\nk̸=i\nwhere terms of order higher than ∆t are neglected. Note that Λ = 0 (no self-interaction). Using Eqs. R2 and\nii,αβ\nR3, we now get,\n(cid:12) (cid:12)2   \n(cid:12) (cid:12)\nd⟨E⟩ 1 (cid:88)(cid:12)(cid:88) (cid:12) 1(cid:88)(cid:88) (cid:88) c (cid:88)(cid:88) (cid:2) (cid:3)\ndt =− γ (cid:12) (cid:12) ∇ xi V ij(cid:12) (cid:12) + 2 Tr∇ xi ∇ xi V ij Λ ki+ 2 Tr ∇ xi ∇ xj V ij (Λ ij ) . (R4)\ni (cid:12)j̸=i (cid:12) i j̸=i k̸=i i j̸=i\nThe first term in Eq. R4 always minimizes the total energy, whereas the second and third terms do not. Hence, the\ntotal energy is not strictly minimized. However, energy minima are a fixed point of the equation since the right hand\nside of the evolution equation vanishes, i.e., d d x t i =0, when ∇ xi (cid:80) j̸=i V(x i ,x j )=0 for each particle i. In the regime\nof weak, perturbative noise (e.g., batch fraction b ∼1), the diffusion matrix is negligible compared to the drift term.\nf\nAs a result, the system primarily follows the gradient-descent drift, and relaxes toward an energy minimum. We\ninclude this discussion of the energy evolution in SI Section 1A.\nResponse to Reviewer 3\nWe thank the reviewer for their comments. Below we address the questions of the reviewer point-by-point.\n1. I appreciate the efforts that the authors have spent in revising the manuscript and the explanation they have offered\nto the questions in my previous report. I agree with the authors that the influence of c is not exactly the same as the\neffect of thermal noise in random organising hyperuniform fluids in the literature, as they come with slightly different\nmathematical forms. But I think the essential physics in the effect of c is the same as adding an extra noise, and in\nthisworkitisawhitenoisebutnotexactlythesamewayasbeingaddedintheactivehyperuniformfluidspreviously. I\nbelieve that the effect of c can be also decoupled as an extra noise term in addition to the noise conserving the center of\nthe mass of the system. Then the finding in this work is not “universal” as claimed by the authors, as the effect should\ndepend on the exact form of the extra noise, such as the color of the noise, the mathematical form of correlation in\nthe noise, etc. Therefore, I don’t think this manuscript is suitable for publications in high profile journals like Nature\nCommunications, but rather suggest the authors to resubmit it (with revision) to more specialised journals.\nWethankthereviewerforacknowledgingthattheinfluenceofthenoisecorrelationcoefficientcinourworkisnotthe\nsame as the effect of thermal noise random organizing hyperuniform fluids in the literature. However, we respectfully\ndisagree with the other two points raised and trust that our explanation below will clarify the merits of our work.\nThe two criticisms that we address are that (i) the essential physics in the effect of c is the same as adding an extra\nnoise,whichinourworkisawhitenoise,and(ii)theeffectofccanbedecoupledasanextranoiseterminadditionto\nthe noise conserving the center of the mass of the system. Then the finding in our work is not “universal” and should\ndepend on the exact form of the extra noise, such as the color of the noise, the mathematical form of correlation in\nthe noise, etc. Below we answer both of these concerns raised by the reviewer.\n(i)Inourwork,wedonothavetwonoisesources. Instead,wehaveasingle drivingnoisesourcewhoseintrinsicnature\niscontrolledbythecorrelationcoefficientc,whichdirectlytunesthedegreeofcorrelationinthepairwiseinteractions,\nfromcompletelyuncorrelated(c=0)toanti-correlated(c=−1). Weshowthattheemergenceoflong-rangestructure\ndepends specifically on the multiplicative and pairwise-correlated nature of this microscopically derived noise term.\nThis is in contrast to systems having a center-of-mass conserving (“Laplacian”) noise in addition to a thermal bath,\nwhere two independent noise sources act on the system [3].\n(ii) Previous work on random organizing hyperuniform fluids has looked at the effect of adding white noise to the\ncenterofmassconserving“Laplacian”noise[3]. Notethatthiscapturestheeffectofaddinganexternalthermalnoise\ntothecenterofmassconservingnoise. Asmentionedabove,inoursystem,wedonot havetwodifferentnoisesources.\nWhile one can phenomenologically construct a continuum description by adding white noise to the center of mass\nconserving “Laplacian” noise, such an approach fails to faithfully capture our system in two key respects. First, such\na description is purely phenomenological, rather than derived by directly coarse-graining the microscopic equations\nof motion, as we do in our work. Second, even at a phenomenological level, this continuum model is fundamentally\ninconsistent, as it implies the presence of two independent noise sources. Indeed, following the phenomenological\nroute renders these alternative frameworks non-universal and makes the choice of the “additional” noise term (its\ncolor, correlation structure, etc.) crucial. Conversely, our system is driven by a single stochastic source that we do\nnot choose because we derive it. In other words, we do not “just add the noise”.\nFinally, the universal active phase behavior in random-organizing systems arises from the observation that, although\nthe microscopic noise sources and their specific forms vary across systems (RO, BRO, SGD), their dynamics can\neffectively be approximated by a generalized model (see Eq. 5 of the main text, and SI. Sec. 1.A for the procedure of\napproximatingthediscrete-timedynamicsbyacontinuous-timeSDE).Consequently,ourresultsareagnostictothespecifictypeorformofnoiseintheoriginaldynamics,\nas long as it can be mapped to the generalized model (Eq. 5). Once again, this is not the result of a choice but the\nresult of a coarse-graining procedure applied to systems with different noisy dynamics that leads to the same SDE\n(Eq. 5). In addition, our results do not depend on choice of short range potential, dimensionality, volume fraction (as\nlong as in the active phase), and kick magnitude. As such, our claim that the behavior is universal is well founded.\n[1] G. Rotskoff and E. Vanden-Eijnden, Communications on Pure and Applied Mathematics 75, 1889 (2022).\n[2] D. S. Dean, Journal of Physics A: Mathematical and General 29, L613 (1996).\n[3] Q.-L. Lei, M. P. Ciamarra, and R. Ni, Science advances 5, eaau7423 (2019)."}, {"round": 3, "reviewer_comments": {"reviewer_1": "The authors have provided satisfactory answers to all my comments. Therefore, I recommend publication of the manuscript\nin its present form.", "reviewer_3": "I have read the response from the authors, which I disagree with. They insisted that there is only one noise in the system\ncontrolling but c with c=0 being uncorrelated and c=-1 anti-correlated, but as I mentioned that it can always be\ndecomposed into two noises with one anti-correlated responsible for the hyperuniformity and the other one responsible\nfor deviation from hyperuniformity. This picture actually agrees with what they have seen in their simulation, which is\nhyperuniform to certain lengthscale and can be seen as a distorted hyperuniformity. Therefore, I think that the authors\noverrated the significance of this work, and I do not think that the phenomena they saw is universal, which should\ndepend on the exact form of the other noise changing the hyperuniformity of the system. So I cannot recommend for\npublication in Nature Communications."}, "author_response": "Response to the reviewers’ report\nWe report below a detailed response to the reviewers’ reports. The original reviewer comments are in italic, and our\nresponse appears in normal font. When addressing the reviewers’ comments, we also explicitly list the changes we\nmade in the revised manuscript and highlight those changes in red here and in the revised manuscript.\nResponse to Reviewer 3\nWe thank the reviewer for their comments. Below we address the questions of the reviewer point-by-point.\n1. I have read the response from the authors, which I disagree with. They insisted that there is only one noise in the\nsystem controlling but c with c=0 being uncorrelated and c=-1 anti-correlated, but as I mentioned that it can always be\ndecomposed into two noises with one anti-correlated responsible for the hyperuniformity and the other one responsible\nfor deviation from hyperuniformity. This picture actually agrees with what they have seen in their simulation, which is\nhyperuniform to certain lengthscale and can be seen as a distorted hyperuniformity. Therefore, I think that the authors\noverrated the significance of this work, and I do not think that the phenomena they saw is universal, which should\ndepend on the exact form of the other noise changing the hyperuniformity of the system. So I cannot recommend for\npublication in Nature Communications.\nWerespectfullydisagreewiththetwopointsraisedbythereviewer. Asdiscussedindetailinourpreviousresponse,we\nagain clarify the merits of our work. The reviewer’s concerns are, (i) “...the noise can always be decomposed into two\nnoises with one anti-correlated responsible for the hyperuniformity and the other one responsible for deviation from\nhyperuniformity. This picture actually agrees with what they have seen in their simulation, which is hyperuniform to\ncertain lengthscale and can be seen as a distorted hyperuniformity”, and (ii) “...I do not think that the phenomena\nthey saw is universal, which should depend on the exact form of the other noise changing the hyperuniformity of the\nsystem”. Below, we address both concerns.\n(i) In our work, we do not have two noise sources. We have a single driving noise controlled by the correlation\ncoefficientc,whichcontinuouslytunesthepairwisecorrelationsfromuncorrelated(c=0)toanti-correlated(c=−1).\nThe presence of only one noise source is evident at all length scales—both at the particle level (the second term in\nEqs. 1,2,3, and 5) and at the continuum level (the third term in Eq. 6). If there were indeed two noise sources,\nEqs. 1−6 would contain two distinct noise terms instead of one, which they clearly do not.\nThe reviewer’s position likely arises from the observationthat one could hypothetically construct a phenomenological\nmodel involving, for example, a Laplacian noise term and a standard divergence noise term, which would produce a\nsimilar form for the static structure factor S(k). However, this line of reasoning is incorrect for two reasons. First,\nphenomenological descriptions are useful when a microscopic understanding is lacking; in our case, we do have a\ncomplete and explicit microscopic formulation. Using a phenomenological model is therefore unjustified and, in this\ncase, wrong. Second, even if one accepts such a phenomenological representation, it would misleadingly suggest the\nexistence of two independent noise sources, whereas our system contains only one (see Eqs. 1−6). Finally, a variety\nof phenomenological constructions could in principle produce the same S(k), with no principled way to determine\nwhich, if any, corresponds to the actual microscopic dynamics. Our microscopic description circumvents all of these\nissues associated with phenomenological descriptions.\nTo make the above distinction concrete, consider the example of anti-correlated noise (c = −1) with an additional\nthermal noise term (Eq. S83). This system truly has two independent sources of noise, which appear separately at\nboth the microscopic level (two distinct noise terms in Eq. S83) and the continuum level (two distinct noise terms\nin Eq. S84). Although the S(k) for this system qualitatively resembles that of a random-organizing system with\narbitraryc̸=−1butwithout thermalnoise(compareFigs.S7and2b), itwouldbeincorrecttotreatthetwosystems\nasequivalentmerelybecausetheirstructurefactorslooksimilar. Onesystemhastwodistinctnoisesources; theother\nhas only one. Their similarity in S(k) does not imply physical equivalence, since by this logic, all systems exhibiting\nsimilar S(k) versus k scaling would have to be considered the same, which is obviously not true.\n(ii) The universal active-phase behavior in random-organizing systems arises from the observation that, although\nthe microscopic noise sources and their specific forms vary across systems (RO, BRO, SGD), their dynamics can\nbe effectively approximated by a generalized model (see Eq. 5 of the main text and SI Sec. 1.A for the procedure used to\napproximate the discrete-time dynamics by a continuous-time SDE). Consequently, our conclusions hold irrespective\nof the specific type or form of the underlying noise (whether random kicks, biased random kicks, or selection noise),\nprovided that the system can be mapped onto this generalized model (Eq. 5). This universality is not a matter of\nchoice, but a direct consequence of their equivalence upon coarse-graining, which leads to the same SDE (Eq. 5).\nIf the results truly depended on the “exact form of the other noise,” as the reviewer suggests, then the long-range\nstructure for RO, BRO, and SGD would differ, since the microscopic origins of noise in the three cases are distinct\n(random kicks in RO, biased random kicks in BRO, and selection noise in SGD; see Eqs. 1−3). The fact that all\nthreesystemsexhibitthesamelong-rangestructureis,infact,strongevidencethattheresultsareindependentofthe\nmicroscopic details of the noise, so long as their dynamics can be effectively approximated by the generalized model.\nThisisnotahand-wavyassertion; itispreciselywhatweobservebothanalyticallyandinsimulations. Moreover,our\nfindings do not depend on the choice of short-range potential, dimensionality, volume fraction (as long as the system\nremains in the active phase), or kick magnitude. Thus, our assertion that the behavior is universal is well supported.\nToconclude,wethankReviewer3fortakingthetimetoreadmultiplerevisionsofourmanuscriptandforarticulating\ntheirconcerns. Weregretthatourpreviousclarificationswerenotpersuasive,butwehopethatthedetaileddiscussion\nabove makes the merits of our theory and the basis for our universality claims fully explicit, and that the manuscript\ncan now be allowed to move forward."}]

Decision letter