Communication: Evidence for non-ergodicity in quiescent states of periodically sheared suspensions

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.¹ to model the experiments by Pine et al.³ 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.¹ 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πσ²/(4A). Each cycle consists in applying a shear transformation with maximum amplitude γ₀ 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.¹ 

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 γ₀ ≲ 2.66.¹ 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.⁵ 

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.

The dynamical shear model of Ref. 1 requires the information whether or not two given disks A and B, at positions (x_A, y_A)^T and (x_B, y_B)^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 ≡ y_A − y_B. 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 r²(t^∗) = σ² 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 Δx² + Δy² ≤ σ², at least one collision occurs (even without application of shear);

• if Δy² > σ², no collision occurs;

• if Δy² ≤ σ², at least one collision occurs if $Δ x ± σ 2 − Δ y 2 ≤ γ 0 Δy$⁠.

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.

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.¹ (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.

We note that mapping from the asymmetric shear considered in Ref. 1 onto the symmetric butterfly described here involves a transform of the maximum shear amplitude and of the displacement region (see Fig. 1 ⁹).

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 γ₀, we can sample with equal probability among all quiescent states using a straightforward Monte Carlo algorithm.¹⁰ 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, k_BT = 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 γ₀. The data are shown in Fig. 4 correspond to a system with γ₀ = 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.⁹ 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 γ₀ values (see summary in Fig. 7). Technical aspects of the pressure measurement are described in the supplementary material.⁹ 

Figure 5 shows the pair distribution function g(x, y) at γ₀ = 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.

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.,⁸ 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 $ρ ( k ) = ∑ j = 1 N exp ( i k ⋅ r j ) 2$⁠, with r_j the point positions and the average is taken over different realisations of the pattern.^14,15 Then hyperuniformity is equivalent to $lim | k | → 0 S ( k ) =0$⁠, see, e.g., Refs. 12 and 13.

For equilibrium systems, the long-wavelength behaviour of S(k) is related to the compressibility via $lim | k | → 0 S ( k ) =∂ρ/∂P$⁠, 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 γ₀ = 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 ⁹). 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 material⁹). 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.⁹ 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 found⁹ 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.

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.

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.