Hyperuniformity and phase separation in biased ensembles of trajectories for diffusive systems

Introduction – Non-equilibrium systems exhibit diverse collective behaviour and complex emergent phenomena, many of which have no counterparts at equilibrium. Even in simple interacting particle systems, one may encounter long-ranged correlations spohn83 , dissipative “avalanche” events with no typical size btw , and dynamical phase transitions bodineau2004 ; garrahan2007 . Theories that capture the universal aspects of these fluctuations are much sought-after, as a route to general descriptions of non-equilibrium phenomena. Here, we analyze non-equilibrium ensembles of trajectories bodineau2004 ; lecomte2005 ; garrahan2007 , defined through constraints on macroscopic observables such as the total current or activity within a given time period. Phase transitions within these ensembles occur when such a constraint leads to a qualitative change in macroscopic behaviour bodineau2004 ; garrahan2007 ; bodineau2008 ; hedges2009 . In diffusive systems bertini2001 ; bertini2005 ; tkl2007 ; hurtado2011-pnas ; bertini-revs , we demonstrate transitions into “hyperuniform” (HU) states torquato2003 , as well as transitions into the macroscopically inhomogeneous (“phase separated”) states that have previously been found bodineau2004 ; bodineau2008 . Hyperuniform states are characterised by anomalously small density fluctuations on large length scales torquato2003 ; florescu2009 ; zachary2011 ; berthier2011 ; man2013 ; chicken2014 ; levine-arxiv ; they have been identified in jammed particle packings berthier2011 ; zachary2011 and in biological systems chicken2014 . These systems are highly optimised in response to a global constraint (mechanical stability in jamming, optimal fitness in biology). The constrained dynamical ensembles that we consider in this study are also optimised: they are the maximally probable states consistent with the constraint. Our results (i) provide further evidence that hyperuniformity is generic, by demonstrating that it occurs in a new set of optimised non-equilibrium ensembles, and (ii) resolve the physical interpretation of some phase transitions that have been previously discovered in diffusive systems appert2008 ; bodineau2008 .

Models – We study biased ensembles of trajectories both computationally and analytically. For computational studies, we consider a one-dimensional model of $N$ diffusing hard particles in a periodic box of size $L$, with each particle having size $l_{0}=1$. This Brownian hard-particle model (BHPM) evolves by Langevin dynamics: the position $x_{i}$ of particle $i$ obeys $\partial_{t}x_{i}=-\beta\nabla_{i}U+\eta_{i}$ where the $\eta_{i}$ are independent white noises, $U$ is the potential energy, and $\beta$ the inverse temperature. The diffusion constant of an isolated free particle is $D_{0}$. We use a Monte Carlo (MC) dynamical scheme to simulate this system. Full system details are given in Appendix A.

We also consider lattice-based exclusion models where $N$ particles are distributed over $L$ lattice sites, again with periodic boundaries. At most one particle may occupy any lattice site. In the partially asymmetric simple exclusion process (PASEP), particles hop left with rate $\ell$ and right with rate $r$, provided their destination site is empty. The symmetric simple exclusion process (SSEP) is the case $\ell=r=1$. The steady states of the BHPM and the PASEP have no correlations between particles beyond hard-core exclusion. (Unlike the SSEP and BHPM, the PASEP does not obey detailed balance, but for periodic boundaries, it may still be shown that site occupancies are uncorrelated in the steady state.)

Biased ensembles of trajectories – Let $K=K[x(t)]$ be a measure of dynamical activity in a trajectory $x(t)$. For exclusion processes, $K$ is the total number of particle hops in a trajectory. For the BHPM, we follow hedges2009 : we choose a coarse-graining time $\tau_{0}$ and focus on trajectories of length $t_{\rm obs}=M\tau_{0}$, defining $K=\sum_{j=1}^{M}\sum_{i=1}^{N}|\hat{x}_{i}(t_{j})-\hat{x}_{i}(t_{j-1})|^{2}$ with $t_{j}=j\tau_{0}$. The position $\hat{x}$ is defined by subtracting the centre-of-mass motion (see Appendix A) which helps to minimize finite-size effects. We take $\tau_{0}=\ell_{0}^{2}/(2D_{0})$, in which time an isolated particle diffuses a distance comparable with its size. We fix the units of time by setting $\tau_{0}=1$.

To investigate trajectories that are constrained to non-typical values of $K$, we define a biased ensemble of trajectories spohn1999 ; lecomte2005 ; garrahan2007 , via a formula for the average of an observable $O$:

Here $\langle\cdot\rangle_{0}$ represents an average in the (unbiased) steady state of the model, $\langle O\rangle_{s}$ is an average within the biased ensemble, and $\psi_{K}(s)=\log\langle\mathrm{e}^{-sK}\rangle_{0}/(L^{d}t_{\rm obs})$ is a ‘dynamical free energy’. For sufficiently large $t_{\rm obs}$, averages in the biased ensemble are equal to averages over trajectories in which the activity $K$ is constrained touchette2013 .

Numerical results for the BHPM – Fig. 1 has results for the BHPM, calculated using transition path sampling tps ; hedges2009 . Fig. 1(a) shows the mean activity $k(s)=(L^{d}t_{\rm obs})^{-1}\langle K\rangle_{s}$. For $s>0$ there is a first-order transition into a phase separated state bodineau2004 ; bodineau2008 ; lecomte2012 . For $s<0$, the activity appears to depend smoothly on $s$, but the system develops strong long-ranged correlations. These are measured by the structure factor $S(q)=\frac{1}{L^{d}}\langle\delta\rho_{q}(t)\delta\rho_{-q}(t)\rangle$ where $\delta\rho_{q}=\int\mathrm{d}r\delta\rho(r)\mathrm{e}^{-iq\cdot r}$ and $\delta\rho(r)=\rho(r)-{\overline{\rho}}$, with ${\overline{\rho}}$ the mean density. To see the relevant behavior most clearly, we transform co-ordinates so that the particles are treated as point-like (see Appendix A): defining $L_{0}=L-Nl_{0}$, the equilibrium ($s=0$) ensemble then has $S(q)=N/L_{0}$, independent of $q$. Fig. 1(b) shows that for $s<0$ and small $q$, the structure factor deviates strongly from this equilibrium value. The signature of a hyperuniform state is that $S(q)\sim q$ at small-$q$ torquato2003 : density fluctuations on large length scales are strongly suppressed. This means that particle positions necessarily have long-ranged correlations [otherwise, self-averaging of the density within large regions of the system implies $\lim_{q\to 0}S(q)>0$]. Analysis of the small-$q$ behaviour in numerical simulations is limited by the system size, but the results for $s<0$ are consistent with hyperuniformity.

Fluctuating hydrodynamics – The BHPM is representative of a general class of diffusive systems, which may be described by “fluctuating hydrodynamics” spohn83 ; eyink1990 ; bertini-revs . Within this theory, the time-evolution of the density $\rho(r,t)$ on large length and time scales can be approximated by a Langevin equation

where $\eta$ is a white noise, $D=D(\rho(r,t))$ and $\sigma=\sigma(\rho(r,t))$ are local measures of diffusivity and mobility, and $a$ is an asymmetric driving force. Details of the relationships between fluctuating hydrodynamics and the BHPM, SSEP and PASEP are given in Appendix B. Note that the fluctuating hydrodynamic theory is valid in all dimensions, not just $d=1$.

Hyperuniformity within fluctuating hydrodynamics – Consider a system described by (2) with $a=0$, and introduce a bias to larger-than-average activity, $s<0$. Averages within the biased ensemble are given by path-integral expressions: $\langle O\rangle_{s}=\mathrm{e}^{-\psi_{K}(s)L^{d}t_{\rm obs}}\int\mathcal{D}\rho\mathcal{D}\hat{\rho}\,O[\rho]\mathrm{e}^{-\int\mathrm{d}r\mathrm{d}t{\cal L}}$, where $\hat{\rho}$ is a (real-valued) response field, and

in which $\kappa=\kappa(\rho)$ is the (density-dependent) local activity of the system. We assume $\kappa^{\prime\prime}(\rho)\leq 0$, which certainly holds for exclusion processes and may be expected to hold for generic particle systems; analysing the case with $\kappa^{\prime\prime}(\rho)>0$ is also straightforward kmp ; lecomte2010 ; hurtado-rev . The behavior of $\kappa(\rho)$ for the BHPM is shown in Appendix A.

Analysis of hydrodynamic behaviour requires a suitable rescaling of space and time co-ordinates. To avoid cumbersome notation we defer this procedure to Appendix B and quote our results in terms of the bare (unrescaled) parameters. Note, however, that these results apply only in the hydrodynamic limit. For $s\leq 0$, the path integral is dominated by trajectories where $\rho(r,t)\approx{\overline{\rho}}$ and $\hat{\rho}\approx 0$ so we write $\rho(r,t)={\overline{\rho}}+\delta\rho(r,t)$ and expand to quadratic order in $\delta\rho$ and $\hat{\rho}$. The result is

where we write $\kappa(\rho)=\kappa_{0}+\kappa_{0}^{\prime}\delta\rho+\kappa_{0}^{\prime\prime}\delta\rho^{2}/2+\dots$, with $\kappa_{0}=\kappa({\overline{\rho}})$, $\kappa_{0}^{\prime}=(d/d\rho)\kappa({\overline{\rho}})$ etc, and similarly for $D(\rho)$ and $\sigma(\rho)$.

The structure factor may then be evaluated (see (bodineau2008, , Equ. (58)) and also Appendix B, yielding

We again emphasise that this result is valid only for small $q\ll 1$, and that $s,\kappa^{\prime\prime}_{0}<0$, by assumption.

Equ. (5) demonstrates a singular response to the field $s$. For $s=0$ and $q\to 0$, the structure factor approaches a non-zero constant $\sigma_{0}/(2D_{0})$, as expected in an equilibrium state with a finite compressibility. However, for any $s<0$, the large scale behaviour changes qualitatively: $S(q)=(q/2)\sqrt{\sigma_{0}/(s\kappa^{\prime\prime}_{0})}+O(q^{2})$. The numerical results of Fig. 1(b) are consistent with this theoretical prediction. Note that hyperuniformity is a large length scale phenomenon: the non-trivial behaviour in $S(q)$ appears only for small $q\lesssim\sqrt{s\sigma_{0}\kappa_{0}^{\prime\prime}}/D_{0}$.

We also calculate the mean activity $k(s)=\langle\kappa\rangle_{s}\approx\kappa_{0}+(\kappa_{0}^{\prime\prime}/2)\langle\delta\rho(r,t)^{2}\rangle_{s}$. Writing $\langle\delta\rho(r,t)^{2}\rangle=\int\frac{\mathrm{d}q}{(2\pi)^{d}}S(q)$, we see that the suppression of $S(q)$ at small $q$ acts to increase $k(s)$ [recall $\kappa_{0}^{\prime\prime}<0$]. Taking $s<0$ and $d=1$, we obtain (see appert2008 and also Appendix B):

which is valid to leading order in $|s|$. Since $k(s)=-\psi^{\prime}_{K}(s)$ where $\psi_{K}$ is the dynamical free energy, we identify this non-analytic behaviour in $k(s)$ with a second order dynamical phase transition. This singular behavior has been noted before appert2008 , but its link with hyperuniformity has not. In $d>1$, the suppression of $S(q)$ at small wavevectors leads to a singular contribution $K(s)-K(s=0)\sim(-s)^{d/2}$, with logarithmic corrections if $d$ is even (see Appendix B). As illustrated in Fig. 1, biasing to lower-than-average activity by choosing $s>0$ instead leads to phase separation bodineau2004 ; bodineau2008 ; lecomte2012 ; hurtado-rev .

Heuristic explanations for HU states – The origin of hyperuniformity in biased diffusive systems is the diverging hydrodynamic time scale associated with large-scale density fluctuations. To see this, consider linear response to the field $s$. Within a biased ensemble of trajectories, the probability of finding the system in a configuration $\mathcal{C}$ is $p_{\mathcal{C}}(s)=p_{\mathcal{C}}(0)[1-2s\int\mathrm{d}r\mathrm{d}t\langle\delta\kappa(r,t)\rangle_{\mathcal{C}}+O(s^{2})]$ where $\langle\delta\kappa(r,t)\rangle_{\mathcal{C}}$ is a “propensity” propensity , which is obtained by averaging the activity over trajectories that start in $\mathcal{C}$ at $t=0$, and comparing with typical equilibrium trajectories garrahan2009 ; jack2014-east .

If $\mathcal{C}$ has an unusual density fluctuation at a small wavevector $q\approx 1/R$, expanding $\delta\kappa$ to quadratic order in $\delta\rho$ gives $\int\mathrm{d}r\delta\kappa(r,t)\approx\frac{\kappa_{0}^{\prime\prime}}{2L^{d}}[|\rho_{q}(t)|^{2}-\langle|\rho_{q}(t)|^{2}\rangle_{0}]$. Diffusive scaling therefore indicates that $\langle\delta\kappa(r,t)\rangle_{\mathcal{C}}\approx\frac{1}{2}\kappa_{0}^{\prime\prime}[|A_{\mathcal{C}}|^{2}-S_{0}(q)]\mathrm{e}^{-t/\tau_{R}}$, where $\tau_{R}=R^{2}/D_{0}$ is a relaxation time, $A_{\mathcal{C}}=\rho_{q}(0)/L^{d/2}$ is the amplitude of the density fluctuation and $S_{0}(q)$ the structure factor of the unbiased state. In hyperuniform states we expect $|A_{\mathcal{C}}|^{2}\ll S(q)$ yielding $(d/ds)\left.p_{\mathcal{C}}(s)\right|_{s=0}\propto\kappa_{0}^{\prime\prime}S_{0}(q)\tau_{R}$: the time scale $\tau_{R}\sim R^{2}$ diverges for large $R$ and $\kappa^{\prime\prime}<0$, so these $p_{\mathcal{C}}$ are strongly enhanced for $s<0$. On the other hand, phase-separated configurations have $R\approx L$ and $A\sim L^{d/2}$ so $(d/ds)\left.p_{\mathcal{C}}(s)\right|_{s=0}\propto-\kappa_{0}^{\prime\prime}L^{d}\tau_{L}$: these configurations have strong (divergent) enhancement for $s>0$ and as $L\to\infty$. Hence, this perturbative analysis reveals an instability of the small-$q$ modes to changes in $s$: we argue that this is the origin of the HU and phase-separated states when $s\neq 0$. The diverging diffusive time scale $R^{2}/D_{0}$ is central to this analysis, similar to other cases where diverging time scales lead to phase transitions in biased ensembles jack2010rom ; garrahan2009 .

Biased ensembles based on the total current – So far, we have considered ensembles of trajectories biased according to their activity $K$. In fact, HU states also appear in ensembles of trajectories where the total current is biased. We define the total current $J$ as the sum of all (directed) particle displacements in a trajectory. (For exclusion processes, this is the difference between the numbers of right- and left- hops.) For generality, we consider jointly-biased ensembles where the activity is biased by a field $s$ and the current $J$ is biased by a field $h$. The analogue of (1) is $\langle O\rangle_{s,h}=\mathrm{e}^{-\psi_{\rm KJ}L^{d}t_{\rm obs}}\langle O\mathrm{e}^{hJ-sK}\rangle_{0}$ (see Appendix B). Within fluctuating hydrodynamics and assuming a time-reversal symmetric unbiased state ($a=0$) one finds that the response to the bias depends only on the quantity $B=s\kappa_{0}^{\prime\prime}-\frac{1}{2}h^{2}\sigma^{\prime\prime}_{0}$ (see Appendix B). For $B=0$ then $S(q)$ has a finite (non-zero) limit as $q\to 0$; for $B>0$ the system is HU while for $B<0$ one has phase separation. The resulting dynamical phase diagram is shown in Fig. 2: the fluctuating hydrodynamic analysis holds only for $s,h\ll 1$, but in the absence of additional phase transitions one expects the same structure to hold throughout the $(s,h)$-plane. We discuss this conjecture below, using results from exclusion processes. We also note in passing that the centre of mass in the BHPM undergoes free diffusion so that the distribution of the current $J$ is trivial in that case: in the fluctuating hydrodynamic theory, this means that $\sigma\propto\rho$, so $\sigma^{\prime\prime}_{0}=0$.

The condition $B=0$ [solid line in Fig. 2(a)] recovers a homogeneous state with $S(q\to 0)>0$: we use the term “normal fluctuations” for this case, in contrast to hyperuniform or phase-separated states. In fact, biased ensembles with $B=0$ are identical to (unbiased) steady states of models in which time-reversal symmetry is broken ($a\neq 0$). This may be verified directly from the hydrodynamic Lagrangian (3) but a clearer interpretation of this result can be obtained by analyzing exclusion processes, as we now discuss.

Mappings between biased ensembles for exclusion processes – We analyse exclusion processes via operator representations of their master equations stinch-review . Starting from the master equation for the SSEP, we write a biased generator spohn1999 ; lecomte2005 that describes the jointly biased ensemble. This operator has a representation in terms of Pauli spin matrices tkl2008 (see Appendix C)

where $n_{i}=\sigma^{+}_{i}\sigma^{-}_{i}$. If we consider instead a PASEP biased by its current, the relevant operator is $W_{\rm A}(\ell,r,\tilde{h})=\sum_{i}r\mathrm{e}^{\tilde{h}}\sigma^{-}_{i}\sigma^{+}_{i+1}+\ell\mathrm{e}^{-\tilde{h}}\sigma^{-}_{i}\sigma^{+}_{i-1}-(r+\ell)n_{i}(1-n_{i+1})$, where $\ell,r$ are hopping rates and $\tilde{h}$ the biasing field. For appropriate parameters (including always the normalization $r+\ell=2$), we may have $W_{\rm S}=W_{\rm A}$, which means that the trajectories of the two biased ensembles are identical. Defining the hopping asymmetry $a_{0}=(r-\ell)/(r+\ell)$, equality between $W_{\rm S}$ and $W_{\rm A}$ requires $a_{0}=\tanh h$ and $\mathrm{e}^{s}=\cosh(h-\tilde{h})$. Any jointly-biased SSEP with $s>0$ leads to two solutions for $\tilde{h}$, which correspond to two possible current-biased PASEPs. These two possibilities are related by a Gallavotti-Cohen symmetry gallavotti-cohen ; spohn1999 . In the specific case $\tanh(\tilde{h}/2)=-a_{0}$ [so that $\mathrm{e}^{\tilde{h}}=\ell/r$], one has $W_{\rm A}(\ell,r,\tilde{h})=W_{\rm A}(r,\ell,0)$: the Gallavotti-Cohen symmetry relates biased and unbiased PASEP states, both of which have “normal” fluctuations.

The SSEP with $\cosh h=\mathrm{e}^{s}$ is another special case, because it leads to $\tilde{h}=0$: the steady state of a jointly-biased SSEP corresponds exactly to that of an unbiased PASEP. This is the microscopic interpretation of the condition $B=0$ in the fluctuating hydrodynamic analysis [for small $h,s$, we obtain $s=h^{2}/2$ which is consistent with $B=0$, because $\sigma(\rho)=\kappa(\rho)$ for the SSEP]. The unbiased steady state of the PASEP has normal fluctuations”, so we conclude that fluctuations are also normal along the line $\cosh h=\mathrm{e}^{s}$ in the dynamical phase diagram of the jointly-biased SSEP [solid line in Fig. 2(a)].

There are a family of mappings between biased SSEPs and PASEPs (see Appendix C): for example, an activity-biased PASEP may also be mapped to a jointly-biased SSEP. The resulting situation is shown in Fig. 2(b,c) where we show the PASEP dynamical phase diagrams that correspond to the (conjectured) phase diagram in Fig. 2(a). The hypothesis is that all points in Fig. 2(a) are either phase separated (PS) or hyperuniform (HU), except for the normal line $\cosh(h)=\mathrm{e}^{s}$. We provide several arguments in support of this picture. The fluctuating hydrodynamic analysis establishes these results in the small bias regime, $|h|,|s|\ll 1$, since the condition $\cosh h=\mathrm{e}^{s}$ then reduces to the case $B=0$ discussed above. The question is therefore whether some other phase transition might intervene and destroy the PS or HU correlations when $h,s=O(1)$.

We are not able to rule out this possibility but several exact results indicate strongly that there is no such phase transition. (i) For $s\to-\infty$, the density correlations of the PASEP are known popkov2011 : independently of the asymmetry $a_{0}$, there is a logarithmic effective interaction potential between particles which renders this state hyperuniform. This implies that the jointly-biased SSEP is HU as $s\to-\infty$ (for all $h$). (ii) For $h=0$, a variational argument garrahan2007 indicates that the SSEP phase separates for all $s>0$: the system can then access configurations where the total number of available hops remains finite as $L\to\infty$. (iii) Phase separation has been shown analytically for the totally asymmetric exclusion process ($a_{0}=1$): this transition corresponds to the appearance of “shocks” in response to a bias on the current bodineau2005 . For the SSEP, this establishes phase separation for all $B>0$ in the limit $h\to\infty$. Combining these results establishes that the proposed phase diagram of Fig. 2(a) is correct in all of the shaded regions: we cannot rule out other phase diagrams that are consistent with these constraints but this simple picture is the most likely scenario. If the proposed Fig. 2(a) is correct, the phase diagrams in Fig. 2(b,c) follow from the exact mappings between biased SSEP and PASEPs.

Conclusion – Fig. 2 indicates that exclusion processes respond very strongly to biases $h$ and $s$, which almost always lead to either phase-separated or hyperuniform states. The normal fluctuations that are familiar from equilibrium systems occur only under special high-symmetry conditions, such as $B=0$. These results provide another example zachary2011 ; chicken2014 ; levine-arxiv of hyperuniformity emerging in non-equilibrium states, and they show that the dynamical phase transition identified in appert2008 corresponds physically to the appearance of hyperuniformity. More generally, the theory of fluctuating hydrodynamics indicates that these dynamical phase transitions should be generic (“universal”) in systems with locally-conserved hydrodynamic variables such as energy or density. The interplay between these phase transitions and the “glass transitions” found previously in biased ensembles of trajectories garrahan2007 ; garrahan2009 ; hedges2009 merits further study – diffusive large-scale behaviour is not a necessary condition for those glass transitions garrahan2007 ; elmatad2010 , but the analysis presented here indicates that phase-separated states may compete with homogeneous glassy states in systems that are biased to low activity.

We thank Fred van Wijland, Vivien Lecomte, and Juan P. Garrahan for many useful discussions. RLJ and IRT were supported by the EPSRC through grant EP/I003797/1.