Breakdown of Ergodicity and Self-Averaging in Polar Flocks with Quenched Disorder

We show that quenched disorder affects polar active matter in ways more complex and far-reaching than believed heretofore. Using simulations of the 2D Vicsek model subjected to random couplings or a disordered scattering field, we find in particular that ergodicity is lost in the ordered phase, the nature of which we show to depend qualitatively on the type of quenched disorder: for random couplings, it remains long-range order as in the pure case. For random scatterers, polar order varies with system size but we find strong non-self-averaging, with sample-to-sample fluctuations dominating asymptotically, which prevents us from elucidating the asymptotic status of order.

Quenched disorder is known to be able to affect qualitatively the asymptotic properties of various systems [1][2][3][4][5][6][7][8][9][10]. Its influence on active matter has recently attracted interest, and rightly so since in many of the corresponding real situations active particles have to avoid obstacles, or move on a rough substrate or in a disordered mesh [11].
While some interesting results were obtained for scalar active matter [12][13][14][15][16][17][18], many of these studies have dealt with the case of dry polar flocks, in continuity with the seminal role played by the Vicsek model and the Toner-Tu theory [19][20][21][22][23][24]. Most efforts were devoted to the fate of the two-dimensional (2D) ordered liquid phase in the presence of quenched disorder. It was found that an optimal amount of noise or disorder can maximize polar order [25][26][27][28]. Experiments studied how flocks of Quincke rollers found in [29] are altered and eventually destroyed by quenched disorder [30,31]. Recently, Toner and Tu [32,33] extended their theory of the homogeneous ordered phase to take quenched disorder into account, predicting in particular quasi-long-range order in 2D. Numerical work has produced partial results compatible with these predictions [25,26,[33][34][35].
The study of disordered systems has a long history outside active matter. Important concepts in this context are ergodicity and self-averaging, which can both be broken by disorder. Ergodicity is lost when multiple configurations coexist for a given sample (realization of disorder). Systems for which spatial and sample averages are not equivalent in the thermodynamic limit are nonself-averaging [36][37][38][39][40][41][42]. It is also known that the type of quenched disorder can make a difference [43]. Somewhat surprisingly, ergodicity, self-averaging, and the influence of the type of quenched disorder have all been largely ignored in the active matter studies published so far [44].
In this Letter, we revisit the problem and show that quenched disorder affects polar active matter in ways more complex and far-reaching than believed heretofore. Using simulations of the 2D Vicsek model, we find that quenched disorder breaks ergodicity and rotational invariance in the ordered phase: several dynamical attractors coexist for a given realization of disorder. In the disordered phase, ergodicity is recovered, but the short correlation length dynamics are organized around an underlying sample-dependent skeleton best revealed in timeaveraged fields. The type of disorder applied influences the above properties, but it also fundamentally change the structure of the phase diagram, self-averaging, and the nature of the ordered phase: A random coupling-(or noise-) strength landscape does not alter the phase diagram nor the self-averaging long-range ordered phase. Random scatterers, on the other hand, deeply modify the layout of the phase diagram, and the Toner-Tu liquid is replaced by a non-ergodic and non self-averaging phase, in which 3 types of fluctuations compete (dynamical/thermal, sample-to-sample, but also between attractors existing for a given sample), a numerically and theoretically challenging situation which prevents us from elucidating asymptotic nature of this phase.
We study extensions of the standard discrete-time Vicsek model with angular noise [19,45,46]. Particles i = 1, ..., N with position r i move at constant speed v 0 and align their velocities with that of current neighbors. We use square domains of linear size L with periodic boundary conditions, divided into unit boxes in which quenched disorder variables are defined. Governing equations read: where n is the index of the unit box containing r t i , U[u] = u/|u| returns unit vectors, R i,t η [u] rotates vector u by a random angle drawn for each particle i at each timestep t from a uniform distribution inside an arc of length 2πη centered on u, and . j∼i is the average over all particles j within unit distance of i (including arXiv:2010.02356v1 [cond-mat.stat-mech] 5 Oct 2020 i). In the RS case (random scatterers, Eq. (1b)), R n ε rotates vectors by some angle defined forever on each box n and drawn from a zero-mean uniform distribution of width 2πε. In the RC case (random coupling-or noisestrength, Eq. (1c)), R i,t η(n) is similar to R i,t η , but the noise amplitude η(n), different in each box n, is drawn once and for all from a uniform distribution over [η, η + ε] [47] Finite-size phase diagrams. The standard Vicsek model possesses two main parameters, the global number density of particles ρ 0 and the (annealed) noise strength η. For large ρ 0 or small η, a polar liquid with true longrange order is observed, whereas only a disordered gas (with short-range correlations) exists at low densities and strong noise. In the (ρ 0 , η) plane, these two phases are separated by a coexistence domain in which dense, ordered bands travel in a sparse gas [24,46,48].
In Eqs. (1), a third parameter is present, the strength of quenched disorder ε. A detailed study of 3-parameter phase diagrams for each type of disorder is a very demanding task. Our efforts have led to the 'stylized' finite-size phase diagrams presented in Fig. 1 (the protocole followed to define them is detailed in [49]). The RC phase diagram looks identical to that of the pure case in the (ρ 0 , η) plane indicating that ε and η may play similar roles in this case ( Fig. 1(a,b)). On the other hand, quenched disorder modifies substantially the layout in the RS case ( Fig. 1(c,d)): the bands region remains present but its extent is bounded away from both low and high ρ 0 or η values at any finite ε. The ordered region is also bounded similarly. These RS results extend and clarify the findings of [25,26].
We now turn to the characterization of the encountered phases, focussing on ergodicity, self-averaging, fluctuations, and memory. We use the modulus and direction of the instantaneous global polar order, m = | e t j j | and θ m = arg e t j j , as well as coarse-grained fields calculated on the unit boxes on which quenched disorder is defined, notably the momentum field m(r, t). Numerical details can be found in [49]. All results presented below were obtained for ρ 0 = 1, as in Fig. 1(b,d).
Most phases are actually qualitatively different from those of the disorderless case. This is even true for RC disorder, even though its phase diagram resembles that of the pure case. The only exception is the coexistence phase: with any type of disorder, we found that its characteristic traveling bands retain their main properties, and notably lead to global long-range polar order [50].
Ergodicity is broken by quenched disorder in the globally ordered regimes found at finite size: for a large fraction of realizations of disorder, and for systems not too small, different initial conditions typically lead to different polarly ordered steady states. However, each sample leads only to a rather small number of attractors, each attracting many different initial conditions (Fig. 2). Each attractor is best characterized by the long-timeaverage of momentum field, m ∞ (r) = lim T →∞ m T (r) with m T (r) = 1 T t=t0+T t=t0 m(r, t). However θ m remains quasi-constant in time in most cases, and different from attractor to attractor, so that following it is sufficient to distinguish them. Quenched disorder thus fixes global order at particular angles, in contrast with the pure case, for which θ m wanders slowly in a diffusive manner ( Fig. 2(a,b)).
To be true, global order continues to wander in small systems with weak quenched disorder. For a given sample, there exists, within the ordered phase, an Ldependent region bordering the pure case inside which ergodicity is not yet broken, located below the dashed lines in Fig. 1(b,d). Increasing L, this region shrinks: 'non-steady', i.e. ergodic, samples dominate at small size, but their fraction quickly decreases, while more and more attractors are found on average ( Fig. 2(g,h)).
While the above observations hold for both RC and RS disorder, there are important differences, notably in the spatial structure of attractors that is much more homogeneous in the RC case than in the RS case (compare Fig. 2(c,d) and (e,f)). Moreover, the global angle of attractors is almost always along the 'easy axes' of the L×L domain, and their number is 4 for large enough L in the RC case ( Fig. 2(a,g)). With RS disorder, attractors have more varied angles, and most often 2 are found in the accessible L range, although their mean number slowly increases with L ( Fig. 2(b,h)).
Memory. Quenched disorder induces permanent memory of the underlying frozen landscape in both the ordered and disordered phases. This is best seen from the existence of well-defined non-trivial longtime averaged fields such as m ∞ (r) and the fact that Edwards-Anderson-like order parameters such as Q EA = lim T →∞ Q(T ) with Q(T ) = m(r, t) · m(r, t + T ) r,t take finite values. In the disordered phase, ergodicity holds and all initial conditions eventually lead to the same longtime dynamics and momentum field m ∞ (r) (Fig. 3(a,b)). In the steady state, Q(T ) converges to some small but finite, sample-dependent Q EA value, in contrast with the pure case disordered phase for which Q(T ) fluctuates around zero (Fig. 3(c)). In the ergodicity-broken ordered phase, Q EA takes rather large values that are not only sample-dependent, but also attractor-dependent in the RS case (not shown).
At a qualitative level, the ergodicity and memory properties presented so far apply to all types of quenched disorder considered here. However, the RC case stands out of the RS case as indicated by the phase diagrams and the structure and statistics of attractors. We now describe in depth how different are their polarly ordered phases.
Fluctuations in the ordered phases. The breakdown of ergodicity in polarly ordered phases implies to consider 3 sources of fluctuations: dynamical, sample-to-sample, and attractor-to-attractor, as we now illustrate in Fig. 4. For a given attractor of a given sample, m fluctuates in time, yielding an asymmetric probability distribution function (PDF(m)). Attractors of a given sample give near identical PDF(m) in the RC case, but not in the RS case ( Fig. 4(a,b)). The sample-and attractor-averaged PDF( m t ) is a very narrow Gaussian in the RC case, but is wide and asymmetric in the RS case (green symbols in Fig. 4(a,b)). This means that for the RC case not only attractor-to-attractor, but also sample-to-sample fluctuations are negligible compared to dynamical ones. In the RS case, on the other hand, neither source of fluctuations can be neglected a priori. To gauge which fluctuations will dominate and the nature of orientational order in the L → ∞ limit, we now turn to finite-size effects. In the RC case, dynamical fluctuations dominate, but all 3 types of fluctuations are in competition in the RS case, which leads to define the following 'connected' (dynamical), 'disconnected' (sampleto-sample), and 'attractor' susceptibilities [51]: where angle and square brackets respectively stand for averages over time and samples, while the upper bar denotes average over attractors. Fig. 4(c) shows how these susceptibilities vary with system size for the ordered RS system considered in Fig. 4(b). While χ con decreases and seems to level off with increasing system size, both χ dis and χ att grow like L α with α ∼ 0.7. This divergence means that the system is strongly non-selfaveraging [52] and that sample-to-sample fluctuations will dominate asymptotically (since χ att χ dis , assuming this behavior holds for L → ∞). As a result, for the system presented in Fig. 4(c), the total susceptibility χ tot = χ con + χ dis + χ att χ con + χ dis first decreases with L but then increases at large sizes when χ dis dominates.
Strong non-self-averaging in the RS case implies that estimating numerically the scaling of the main global polar order parameter M = [ m ] is numerically challenging. Fig. 4(d) shows M (L) for the same parameters as in the rest of the figure, after averaging over typically 1000 samples and recording m for millions of timesteps after transients (see numerical details in [49]). Whereas M (L) decreases slower than a powerlaw in the RC case, indicating true long-range polar order, it decreases faster in the RS case. The local slope (exponent) σ(L) of these loglog plots (Fig. 4(e)) goes to zero as L −2/3 in the RC case (similarly to the pure case [24,53]). In the RS case, σ first decreases slightly, levels off, but then increases: a simple quasi-long-range order (algebraic decay of M , constant σ) is excluded.
Asymptotic nature of the quasi-ordered phase in the RS case. Scanning the whole phase at fixed ρ 0 and η varying ε clarifies the situation described above at a single ε value without bringing definitive answers. Fig. 5(a,b) show σ(L) and χ tot (L) at various ε values. At very small ε, σ and χ tot first decay with L like in the pure case, then depart from this trend at a crossover scale which decreases with increasing ε. Confirming our data in Fig. 4 obtained for a particular ε value, we find no evidence at any ε of 'simple' quasi-long-range order in the range of scales studied: Once σ has stopped decreasing, it does not really plateau and starts increasing slowly.
At the largest ε values considered, σ and χ tot increase with L, then σ levels off at +1, the value characteristic of the short-range order of the disordered phase (M ∼ 1/L), while χ tot decreases. The maximum of χ tot is a measure of the correlation length, which seems to diverge when ε decreases. Correspondingly, χ tot L 2 exhibits a maximum as a function of ε, whose height scales as L γ/ν with γ ν 1.9, while its location ε * (L) seems to converge to a finite asymptotic value ( Fig. 5(d,e)). All this points to a continuous phase transition separating the disordered phase from the quasi-ordered one. But sample-to-sample fluctuations diverge with L all over the quasi-ordered regimes with the same exponent as reported above, χ dis ∼ L 0.7 (Fig. 5(c)). These fluctuations having barely started to dominate χ tot at the largest accessible scales (Fig. 5(b)), it is premature to conclude about the asymptotic nature of the phase corresponding to the quasi-ordered states observed at finite size, and a fortiori about the transition.
To summarize, quenched disorder affects polar active matter in more ways than believed so far. In particular, it breaks ergodicity in the ordered regimes observed at finite size, but not in the disordered phase, which only shows infinite memory of the frozen landscape. These results should be observable experimentally, e.g. in the Quincke roller system of [29,54,55], which is believed to be a realization of (effectively) dry polar active matter. However, our findings seem to contradict the conclusions of [31]: there the breakdown of polar flocks was argued to lead to a "dynamical vortex glass" with many coexisting attractors, while we have shown that sufficiently long averages reveal an ergodic disordered phase with infinite memory (Fig. 3). We hope that further experiments will clarify this important point.
We also showed that the nature of the ordered phase depends qualitatively on the type of quenched disorder: for random couplings, it retains the true long-range order well known in the pure case. For random scatterers (and random fields [50]), the ordered phase is strongly nonself-averaging, with sample-to-sample fluctuations dominating asymptotically. Unfortunately, this asymptotic regime is largely inaccessible numerically, so that we are unable to conclude. Theoretical work is thus necessary, which, in our opinion, should take ergodicity and selfaveraging issues into account.