We present simulations of an equilibrium statistical-mechanics model that uniformly samples the space of quiescent states of a periodically sheared suspension. In our simulations, we compute the structural properties of this model as a function of density. We compare the results of our simulations with the structural data obtained in the corresponding non-equilibrium model of Corté et al. [Nat. Phys. 4, 420 (2008)]. We find that the structural properties of the non-equilibrium model are very different from those of the equilibrium model, even though the two models have exactly the same set of accessible states. This observation shows that the dynamical protocol does not sample all quiescent states with equal probability. In particular, we find that, whilst quiescent states prepared in a non-equilibrium protocol can be hyperuniform [see D. Hexner and D. Levine, Phys. Rev. Lett. 114, 110602 (2015); E. Tjhung and L. Berthier, Phys. Rev. Lett. 114, 148301 (2015); and J. H. Weijs et al., Phys. Rev. Lett. 115, 108301 (2015)], ergodic sampling never leads to hyperuniformity. In addition, we observe ordering phase transitions and a percolation transition in the equilibrium model that do not show up in the non-equilibrium model. Conversely, the quiescent-to-diffusive transition in the dynamical model does not correspond to a phase transition, nor a percolation transition, in the equilibrium model.
I. INTRODUCTION
Gibbsian statistical mechanics is based on the assumption that the time-averaged properties of a system are equal to the ensemble averages. Although this “ergodic” hypothesis is not valid for all classical many-body systems, it is generally believed that it is correct for thermal equilibrium systems. However, for non-equilibrium systems, there is no reason to expect ergodicity to hold. It is, however, not straightforward to test the breakdown of ergodicity. To carry out such a test, we need to know exactly what parts of the phase space are accessible to the non-equilibrium system, and then test whether those regions are sampled with the same frequency both in and out of equilibrium.
Here we describe simulations of a simple model system where it is possible to test the ergodicity of a non-equilibrium system. To be more precise, we consider both a simple non-equilibrium system and an equilibrium system that has exactly the same accessible states as the non-equilibrium system. The non-equilibrium model that we consider was first introduced by Corté et al.1 to model the experiments by Pine et al.3 on periodically sheared suspensions of large (non-Brownian) particles. These experiments explored the time evolution of a system of spherical plastic particles of diameter ca 0.23 mm, suspended in a viscous, density-matched fluid. The fluid was confined between coaxial cylinders that were sheared periodically. Pine et al. found that below a critical maximum shear amplitude, after a transient period, the particle trajectories are reversible, i.e., the particles return to their original positions at the end of each shear cycle.
The resulting “quiescent” arrangements of particles lack obvious geometrical order, yet they are clearly organised in such a way that particles avoid collision during a shear cycle. Corté et al.1 introduced a simple model of the shear experiment, which starts by placing N disks of diameter σ uniformly at random in a rectangular box of area A with periodic boundary conditions, giving a volume fraction of ϕ = Nπσ2/(4A). Each cycle consists in applying a shear transformation with maximum amplitude γ0 in, say, horizontal (x-) direction and recording for each particle the number of other particles it encounters during the cycle (see Fig. 1). Then the box is restored to its initial shape and each particle is displaced randomly once for each encounter recorded during the cycle. This model reproduces the threshold behaviour observed in experiments.1,3 Below the threshold shear amplitude, the model system reaches a quiescent state. Corté et al. coined the phrase “random organisation” to describe such states.1
Horizontal shear transformation of a box. The maximum shear amplitude γ0 is given by the ratio of horizontal displacement h over box height v. The model of Ref. 1 is shown on the left, the right shows the symmetric version considered here. By shearing in one cycle both, to the left and the right, by half the amplitude one obtains the same behaviour.
Horizontal shear transformation of a box. The maximum shear amplitude γ0 is given by the ratio of horizontal displacement h over box height v. The model of Ref. 1 is shown on the left, the right shows the symmetric version considered here. By shearing in one cycle both, to the left and the right, by half the amplitude one obtains the same behaviour.
For shear amplitudes above a certain volume-fraction-dependent threshold value, particles keep colliding and random organisation into quiescent states does not occur. As an illustration, we show a typical example of the simulation results obtained with Corté’s model: Figure 2 shows the fraction of particles that are displaced during each cycle as function of the number of cycles. At a volume fraction ϕ = 0.2, quiescent states are observed for γ0 ≲ 2.66.1 The existence of such a threshold separating reversible from non-reversible behaviour has also been reproduced in more sophisticated Stokesian dynamics simulations.3,4 However, the key non-equilibrium dynamics is already contained in the model of Ref. 1. Corté’s model has been argued to be in the universality class of conserved directed percolation.5
Dynamical model: fraction of active (colliding) particles in the shear cycle as function of the number of shear cycles, for different maximum shear amplitudes γ0, at volume fraction ϕ = 0.2, for N = 1000 particles. As reported in Ref. 1, the critical shear amplitude is γ0,c(0.2) ≈ 2.66. For γ0 < γ0,c, the activity dies out and particles arrange in quiescent states. Above the threshold, γ0 > γ0,c, the activity appears to persist indefinitely. At the critical point, the fraction of active particles decays as a power law with the number of cycles.1,2
Dynamical model: fraction of active (colliding) particles in the shear cycle as function of the number of shear cycles, for different maximum shear amplitudes γ0, at volume fraction ϕ = 0.2, for N = 1000 particles. As reported in Ref. 1, the critical shear amplitude is γ0,c(0.2) ≈ 2.66. For γ0 < γ0,c, the activity dies out and particles arrange in quiescent states. Above the threshold, γ0 > γ0,c, the activity appears to persist indefinitely. At the critical point, the fraction of active particles decays as a power law with the number of cycles.1,2
The properties of the model of Corté et al. have been studied extensively. In particular, recent studies have shown that the quiescent states that emerge at sub-critical shear amplitudes are hyperuniform, meaning that low-wavevector density fluctuations vanish in the limit that the wave-vector goes to zero.6–8 Whilst this observation is intriguing, it is not immediately obvious whether this behaviour is a property of the quiescent states, in the sense that it would show up if we were to sample uniformly over all quiescent states, or whether the non-equilibrium dynamics does not sample the quiescent states uniformly. In order to elucidate this question, we construct an equilibrium equivalent of the non-equilibrium model of Corté et al. As we argue below, our simulations show that different quiescent states are not sampled with equal probability by the dynamical protocol and that this “non-ergodicity” is at the origin of the observed hyperuniformity.
The remainder of this paper is organised as follows. Section II introduces the equilibrium model that we study. In this model, all allowed states are quiescent and all quiescent states are allowed. In Sec. III we report equilibrium Monte Carlo studies that sample quiescent states of the equilibrium model. Section IV considers the nature of density fluctuations in these states.
II. EQUILIBRIUM “BUTTERFLY” MODEL
The dynamical shear model of Ref. 1 requires the information whether or not two given disks A and B, at positions (xA, yA)T and (xB, yB)T, will overlap during the next shear cycle. This can be determined by a straightforward geometrical argument, since the relative shear displacement of disks A and B only depends on their vertical distance, Δy ≡ yA − yB. Parametrising the trajectories which the particles trace out in one shear cycle by t, such that t = 0 corresponds to no shear and |t| = 1 to the maximum shear amplitude, the vector pointing from the centre of B to the centre of A is
We note that Δx and Δy are the distances before applying shear and they are therefore independent of t. The particles, assumed to have the same diameter σ, collide during the shear cycle if |r(t)| ≤ σ for some −1 ≤ t∗ ≤ 1. Setting r2(t∗) = σ2 gives a quadratic equation for t∗ which has no solutions if |Δy| > σ. However if the vertical distance |Δy| ≤ σ, one finds
Given the requirement that |t∗| ≤ 1, one can check for collision between disks A and B as follows:
if Δx2 + Δy2 ≤ σ2, at least one collision occurs (even without application of shear);
if Δy2 > σ2, no collision occurs;
if Δy2 ≤ σ2, at least one collision occurs if .
Figure 3 visualises how this results in a butterfly-shaped domain around particle A such that A and B will collide if and only if the centre of B is in that domain.
Sketch of shear butterfly geometry for unit particle diameter, σ = 1. The central area (red, grey) indicates the particle, the surrounding ring (green, light grey) is excluded due to hard overlap at zero shear, and the surrounding butterfly domain (blue, dark grey) is excluded as a result of one full shear cycle of maximum shear amplitude .
Sketch of shear butterfly geometry for unit particle diameter, σ = 1. The central area (red, grey) indicates the particle, the surrounding ring (green, light grey) is excluded due to hard overlap at zero shear, and the surrounding butterfly domain (blue, dark grey) is excluded as a result of one full shear cycle of maximum shear amplitude .
One can then place a given number of “butterflies” in a box and randomly displace the overlapping ones to simulate the dynamical model of Corté et al.1 (see Sec. I). Note that the butterfly model is not a hard-core model in the conventional sense, as the butterfly around A excludes the centre of B and conversely. Hence, overlaps of butterfly “wings” are allowed, as long as the centre of a butterfly remains outside the body of the other butterfly.
III. EQUILIBRIUM SAMPLING OF QUIESCENT STATES
Configurations in which no butterfly domain contains the centre of mass of any other particle are by construction quiescent states. Therefore, for fixed N, ϕ, and γ0, we can sample with equal probability among all quiescent states using a straightforward Monte Carlo algorithm.10 The Hamiltonian in this Monte Carlo algorithm is such that quiescent states have zero potential energy, whereas configurations with at least one overlap, have infinite potential energy.
We measure the compressibility factor, Z = P/ρ (in this paper, kBT = 1), where P is the pressure and ρ = N/A the number density, as function of the volume fraction ϕ, for different values of the maximum shear amplitude γ0. The data are shown in Fig. 4 correspond to a system with γ0 = 2.66. The equation of state exhibits a flat region around ϕ ≈ 0.4, which seems correlated with the density-driven transition to a striped phase, see also Fig. 5. Of course, hard disks exhibit no such striped phase.9 The transition to the striped phase occurs at a substantially higher volume fraction than the value ϕc ≈ 0.2 above which quiescent states disappear in the dynamical model (see Fig. 2): in other words, the quiescent striped phase is never observed in the dynamical model. We find the same behaviour for a range of γ0 values (see summary in Fig. 7). Technical aspects of the pressure measurement are described in the supplementary material.9
Equation of state of the equilibrium butterfly fluid: The compressibility factor P/ρ is shown as function of the volume fraction ϕ at maximum shear amplitude γ0 = 2.66, for different number of butterflies N. Note that the equilibrium equation of state is smooth at ϕ = 0.2, the density where the transition between the quiescent and diffusive states takes place in the dynamical model.
Equation of state of the equilibrium butterfly fluid: The compressibility factor P/ρ is shown as function of the volume fraction ϕ at maximum shear amplitude γ0 = 2.66, for different number of butterflies N. Note that the equilibrium equation of state is smooth at ϕ = 0.2, the density where the transition between the quiescent and diffusive states takes place in the dynamical model.
Pair-distance distribution (γ0 = 2.66 and N = 64): (a) shows g(x, y) vs the horizontal, x, and vertical, y, centre of mass distances (see Fig. 3), at volume fraction ϕ = 0.4. (b) shows the vertically averaged data g(x, ϕ) vs x, for different ϕ. (c) shows the horizontal average g(y, ϕ) vs y, for the same ϕ as in (b).
Pair-distance distribution (γ0 = 2.66 and N = 64): (a) shows g(x, y) vs the horizontal, x, and vertical, y, centre of mass distances (see Fig. 3), at volume fraction ϕ = 0.4. (b) shows the vertically averaged data g(x, ϕ) vs x, for different ϕ. (c) shows the horizontal average g(y, ϕ) vs y, for the same ϕ as in (b).
Critical volume fractions ϕc of shear models as function of the maximum shear amplitude γ0. Squares indicate the continuum percolation threshold of interpenetrable butterflies.9 The location of the inflection point in the equation of state of the equilibrium butterfly fluid, as discussed in Sec. III (see Fig. 4), is indicated by circles. Triangles indicate the non-equilibrium critical point of the dynamical model of Ref. 1 (see, e.g., Fig. 2).
Critical volume fractions ϕc of shear models as function of the maximum shear amplitude γ0. Squares indicate the continuum percolation threshold of interpenetrable butterflies.9 The location of the inflection point in the equation of state of the equilibrium butterfly fluid, as discussed in Sec. III (see Fig. 4), is indicated by circles. Triangles indicate the non-equilibrium critical point of the dynamical model of Ref. 1 (see, e.g., Fig. 2).
Figure 5 shows the pair distribution function g(x, y) at γ0 = 2.66 for different volume fractions ϕ. We measure g(x, y) by counting the number of particle pair distances, binning them according to displacement in x- and y-directions and normalising by the number of pair distances expected if the particle distribution were homogeneous with density ρ. Horizontal stripe formation is observed at ϕ = 0.4 in Fig. 5(a). Vertical and horizontal averages of this data are shown in Figs. 5(b) and 5(c). As the figure shows, any horizontal ordering decays within a few particle diameters, while periodic ordering persists in the vertical direction.
IV. HYPERUNIFORMITY AND NON-ERGODICITY
Anomalous density fluctuations of sheared suspensions have been studied theoretically, for the model of Ref. 1 and related ones.6,7 Experimentally, the density fluctuations have been studied in microfluidics experiments by Weijs et al.,8 who also performed simulations of their experimental system. Over a narrow range of driving amplitudes, these studies all find evidence for hyperuniformity, i.e., the structure factor S(k) vanishes as |k| → 0,12,13 in some parameter regimes. Explicitly, the structure factor S(k), at wave vector k ≠ 0, for an ensemble of N points is defined by S(k) = 〈ρ(k)〉/N, where , with rj the point positions and the average is taken over different realisations of the pattern.14,15 Then hyperuniformity is equivalent to , see, e.g., Refs. 12 and 13.
For equilibrium systems, the long-wavelength behaviour of S(k) is related to the compressibility via , which, away from close packing, is finite. The results of Refs. 6–8 are therefore indicative of the non-equilibrium nature of quiescent, randomly organised states. However, what the existing tests cannot tell us is whether hyperuniformity is related to the average behaviour of quiescent states, or that it is due to the non-ergodic sampling of quiescent states. As we argue below, our simulations show that it is non-ergodicity that leads to hyperuniformity.
Figures 6(a) and 6(b) show S(k) measured for the equilibrium and dynamical models, for γ0 = 2.66, at the critical volume fraction of the dynamical model, ϕ = 0.2, as well as for smaller and larger ϕ. We also computed S(k = 0) from ∂ρ/∂P (see Fig. 4 9). Clearly, S(k = 0) is non-zero, and hence the equilibrium model is not hyperuniform. By contrast the critical absorbing states of the dynamical model seem to be hyperuniform, as reported in Ref. 6 (see also supplementary material9). However, in the high-k region S(k) is similar for both models. It seems that these short-wavelength similarities are also present in S(k), i.e., without radial averaging.9 We note that for the dynamical model at high volume fractions ϕ ≳ ϕc ≈ 0.2, a finite fraction of particles remains colliding, such that non-quiescent states are also found9 and enter the computation of S(k). Figure 6(c) compares the low-k region for both models with the results for hard disks at different ϕ and for the ideal gas. At low ϕ, the dynamical model behaves similarly to the ideal gas, increasing ϕ towards ϕc ≈ 0.2, S(k → 0) drops and it increases again for higher ϕ > ϕc.
(a) Structure factor S(k) as function of wave vector magnitude k for the dynamical and the equilibrium model with γ0 = 2.66 at ϕ = 0.2 (which is the critical volume fraction of the dynamical model1). The solid black line shows the result for S(k → 0), obtained by computing ∂ρ/∂P from our independent measurement of the equation of state of the equilibrium butterfly fluid. (b) compares S(k) for the dynamical model at different volume fractions, both, above and below ϕc. The solid arrows indicate the values of S(k → 0) in equilibrium, obtained from ∂ρ/∂P. (c) Structure factor low-k limit S(k → 0) as function of the volume fraction ϕ for four different models. For the dynamical and the equilibrium models, we plot S(kmin), where kmin is the smallest accessible wave-vector. The solid green line is computed from ∂ρ/∂P with the hard disk equation of state of Ref. 11. The dashed line indicates the ideal gas result S(k) = 1.
(a) Structure factor S(k) as function of wave vector magnitude k for the dynamical and the equilibrium model with γ0 = 2.66 at ϕ = 0.2 (which is the critical volume fraction of the dynamical model1). The solid black line shows the result for S(k → 0), obtained by computing ∂ρ/∂P from our independent measurement of the equation of state of the equilibrium butterfly fluid. (b) compares S(k) for the dynamical model at different volume fractions, both, above and below ϕc. The solid arrows indicate the values of S(k → 0) in equilibrium, obtained from ∂ρ/∂P. (c) Structure factor low-k limit S(k → 0) as function of the volume fraction ϕ for four different models. For the dynamical and the equilibrium models, we plot S(kmin), where kmin is the smallest accessible wave-vector. The solid green line is computed from ∂ρ/∂P with the hard disk equation of state of Ref. 11. The dashed line indicates the ideal gas result S(k) = 1.
V. CONCLUSIONS
In our study, we have compared an equilibrium model that uniformly samples all quiescent states of a model for a periodically sheared suspensions of non-diffusing disks with the corresponding dynamical model. We find that, in the dynamical model, quiescent states are sampled non-ergodically. The differences between the dynamical and the equilibrium model are so pronounced that neither the phase transition nor the continuum percolation threshold in the equilibrium model is related to the dynamical phase transition in the periodically sheared system.
The different critical volume fractions are summarised in Fig. 7.
We found that critical quiescent states are only hyperuniform if sampled by the non-equilibrium protocol of Ref. 1 and not if sampled ergodically. We note that, in principle, quiescent states can “nucleate” from diffusing states, even above the dynamical transition threshold. It would be interesting to study the pathway for such a nucleation process because our results indicate that the nucleation “barrier” is determined by non-ergodicity, rather than by the lack of availability of quiescent states.
Acknowledgments
This work was supported by ERC Advanced Grant No. 227758 (COLSTRUCTION), by EPSRC Programme Grant No. EP/I001352/1, and by the Swiss National Science Foundation (Grant Nos. P2EZP2-152188 and P300P2-161078). K.J.S. acknowledges useful discussions with Nuno Araújo, Tristan Cragnolini, Daphne Klotsa, Erik Luijten, Stefano Martiniani, Anđela Šarić, Iskra Staneva, Jacob Stevenson, and Peter Wirnsberger. Additional data related to this publication are available at the University of Cambridge data repository, https://www.repository.cam.ac.uk/handle/1810/252964.