Emergence of cooperative dynamics in fully packed classical dimers

We study the behavior of classical dimer coverings of the square lattice - a paradigmatic model for systems subject to constraints - evolving under local stochastic dynamics, by means of Monte Carlo simulations and theoretical arguments. We observe clear signatures of correlated dynamics in both global and local observables and over a broad range of time scales, indicating a breakdown of the simple continuum description that approximates well the statics. We show that this collective dynamics can be understood in terms of one-dimensional"strings"of high mobility, which govern both local and long-wavelength dynamical properties. We introduce a coarse-grained description of the strings, based on the Edwards-Wilkinson model, which leads to exact results in the limit of low string density and provides a detailed qualitative understanding of the dynamics in all flux sectors. We discuss the implications of our results for the dynamics of constrained systems more generally.


I. INTRODUCTION
Dimer models are archetypal systems for the study of the effects of strong local constraints [1,2]. Despite their simplicity, classical dimer models on bipartite lattices exhibit a number of interesting phenomena, such as macroscopic groundstate degeneracy, topological order, and deconfinement of monomers [3]. Their static properties are well understood, and are captured by an effective coarse-grained theory, involving a height field [4] in two dimensions (2D) or an effective gauge field in higher dimensions [5]. In either case, the result is a critical equilibrium phase, with power-law correlations between local degrees of freedom. This class of systems provides the simplest examples of "exotic" thermodynamic behavior purely determined by entropy [3].
Even when the thermodynamic properties of a system have a simple effective description, its dynamics can be more intricate and interesting [6]. It is natural to ask whether this is the case for dynamical extensions of the classical dimer model, about which much less is known. In particular, Ref. [7] considered the simplest extension of the coarse-grained description to dynamics, predicting simple relaxational decay of correlations, while Ref. [8] considered the dimer model with nonlocal loop dynamics. These works should be contrasted with studies of defect-driven dynamics in dimer models [9] and of monopole dynamics in spin ice [10,11].
In this work, we consider local stochastic dynamics in the defect-free square-lattice dimer model, using both simulations and theoretical arguments. As far as we are aware, what we report here are the first systematic simulation results for the natural dynamics in these systems, i.e., one of locally flipping plaquettes. We show that the simple continuum picture can fail to describe the true physics even over long time scales and that the phenomenology is, in fact, far richer than the simplicity of the model would suggest.
The first main contribution of this paper is to demonstrate, using simulations, significant deviations from exponential relaxation in global and local observables over a broad range of time scales. We argue that a simple continuum description fails because the dynamics is facilitated by local objects-in this case, one-dimensional strings [12][13][14]-and hence highly heterogeneous. The understanding of the importance of these objects, which has broad implications for the study of cooper-ative dynamical phenomena, is our second main contribution.
Close-packed dimer models obey a topological constraint [3] that amounts to conservation of strings, or, equivalently, of the flux of an effective magnetic field. Any local rearrangement of dimers conserves flux, and we exploit this by considering dynamics within a fixed flux sector. At large flux, the system is spanned by a low density of strings, whose fluctuations govern the relaxation. We introduce a coarse-grained description of the strings, based on the Edwards-Wilkinson model of fluctuating interfaces [15], from which we derive exact expressions for dynamical observables, constituting the third main contribution of this work. We confirm these results using simulations at large flux, to which they can be compared with at most one adjustable parameter.
We find that the behavior is qualitatively similar, and consistent with the string picture, for smaller flux, and even in the isotropic limit of vanishing flux. In these latter cases, strings can still be defined and the string picture remains illuminating, even though their density is so high that they cannot be treated as independent. By application of dynamical scaling theory, we furthermore predict a crossover to the Coulomb-phase results of Ref. [7] at a time scale that diverges on approaching the critical point at saturation flux (zero string density).

Outline
In Section II, we introduce the dimer model and the local dynamics that we study in the remainder of the paper. We also briefly review, in Section II D, the coarse-grained theory introduced by Henley [7] to describe dynamics of the height field, which predicts exponential decay of correlations. The majority of our original results are presented in Section III, where we use an effective theory of string dynamics based on the Edwards-Wilkinson equation to derive results for correlations and for the persistence, a local probe of dynamics. We conclude in Section IV with a brief discussion of the broader significance of our results for dynamics in strongly constrained systems. Some technical details and additional simulation results are presented in appendices. We study a dimer model on an L × L square lattice with periodic boundaries. The occupation variable d µ (r) gives the number, 0 or 1, of dimers on the link joining sites r and r + δ µ , where µ ∈ {x, y} and δ µ is a lattice vector. A configuration is allowed only if every site is occupied by a single dimer. We refer to a plaquette as flippable when it contains two parallel dimers; the flippability f of a configuration is defined as the proportion of plaquettes that are flippable, We define the effective "magnetic field" , where ε r = ±1 on the two sublattices [3]. The constraint on dimer configurations then becomes Gauss' law, div r B = 0, where is the lattice divergence. The flux Φ corresponding to B can be defined by Φ µ = r B µ (r) = r ε r d µ (r). Because ε r+δ ν = −ε r , a pair of neighboring parallel dimers, of either orientation, gives zero net contribution to Φ, and so plaquette-flip dynamics conserves the flux.

A. Staggered configurations and strings
The flux is maximized by a staggered dimer configuration. For example, if d y (r) = 0 for all r, d x (r) = 1 for ε r = +1, and d x (r) = 0 otherwise, then Φ = +Φ maxx , where Φ max = 1 2 L 2 . The other three staggered configurations, related by symmetry, have flux of the same magnitude, |Φ| = Φ max , along other lattice directions. We define the flux relative to its maximum by φ = Φ/Φ max and, for flux along the x direction, the deviation from maximum θ = 1 − φ x .
To reduce the flux from maximum, one can shift a row of dimers spanning the system, which changes Φ x by −L. We refer to such a set of shifted dimers as a string [13,16]. After the shift, plaquettes along its length become flippable; flipping these deforms the string but conserves the flux. A possible path for a single string and the resulting configuration are shown in Fig. 1. With N s strings introduced into a staggered configuration, θ = 2N s /L; the linear density is therefore 1 2 θ.

B. Height mapping
The constraint div r B = 0 can be resolved by defining the height z on each plaquette [7], in terms of which B µ (r) = − 1 4 curl (r,µ) z, where the curl is the difference between the plaquettes on each side of a link. Global shifts of z do not affect B, corresponding to the gauge redundancy in 3D [3]. With an appropriate gauge choice, flipping a plaquette modifies z only on that plaquette [7]. If B has periodic boundary conditions, z(r + Lδ µ ) = z(r) + 4L −1 ν µν Φ ν , where is the Levi-Civita tensor. The spatial average of the derivative of the height (the tilt) is therefore intensive and proportional to φ. We define ζ(r) = z(r)−2 µν µν r µ φ ν , with periodic boundary conditions.

C. Dynamics
The most natural dynamics for the dimer system is one where individual plaquettes flip randomly. This dynamics is efficiently implemented numerically via continuous-time Monte Carlo (MC) [17], in which, when flippable, plaquettes flip according to a Poisson process with rate constant γ. Dynamics at equilibrium within a sector of fixed flux φ can be studied by starting from a fixed configuration and equilibrating using plaquette-flip dynamics. We will denote by · · · an average both over the equilibrium ensemble (where all allowed states with flux φ have equal weight) and, where applicable, over subsequent trajectories.
The dynamical correlation function of the height is defined by G ζ (q, t) = ζ (q, t)ζ(−q, 0) whereζ is the Fourier transform of ζ. We also consider the persistence p(t), the proportion of plaquettes that have not flipped at any point up to time t, which provides a local probe of the evolution.

D. Continuous height-field theory
The static properties of the dimer model can be described by a continuum theory in terms of the coarse-grained height h(r) resulting from averaging ζ(r) over short length scales [3]. Apart from terms irrelevant at long distances, the effective dimensionless free energy is for φ alongx, implying correlations h (q)h(−q) = [ω(q)] −1 for the Fourier transformh, where ω(q) = K y q 2 x + K x q 2 y . The simplest extension of the continuum description to dynamical properties is the Langevin equation [ where the noise has correlations implying exponential decay at long time scales. We show below that this prediction can break down, even when the heightfield approach is accurate for the statics, due to cooperative effects that dominate the dynamics.

III. STRING DYNAMICS
The collective character of the dynamics can be uncovered by considering the behavior near maximal flux (i.e., FIG. 1. Left: A dimer configuration with a single string, relative to a fully staggered configuration with maximal flux Φ x = 1 2 L 2 along the horizontal direction. A single string spanning the system once reduces the flux Φ x by −L, the smallest possible amount, irrespective of its path. Flippable plaquettes, which appear when the dimers are shifted, are marked with stars ( ); the fully staggered configuration has none. Center and right: Configurations with φ x = 1 4 evolved for time t. Persistent plaquettes are white, while those that have flipped are blue; strings are yellow. The dynamics is spatially heterogeneous: Even at relatively long times, extended regions are unvisited by strings and hence persistent.
for small θ), where most of the system is unflippable. This regime can be understood in terms of a low density of wellseparated strings. We first consider the dynamics of a single string, using a continuum description based on the Edwards-Wilkinson equation, before turning to the consequences for the two classes of observables, correlation functions and the persistence.
For a single string traversing the system horizontally, let y(x, t) be the vertical position at horizontal position x and time t. By using a transfer matrix to enumerate all possible string configurations (see Ref. [14] and Appendix A), we find the following two exact results regarding the equilibrium distribution of a single string: (i) The mean number of flippable plaquettes is given by

A. Edwards-Wilkinson equation
At large length and time scales, we expect y(x, t) to obey the Edwards-Wilkinson equation [15,18], where , and Λ and D parameterize, respectively, the stiffness of the string and the strength of the noise. The average · · · 0 is taken over trajectories starting from a given initial configuration y(x, 0).
Accounting for the periodicity in the x direction (but not in y), the Green function for Eq. (7) is where kL/2π ∈ Z. To calculate the two-time correlation function in the equilibrium ensemble, we take both times to infinity with their difference finite, The typical width of a string in equilibrium can be characterized by the mean-square displacement between the points x and 0 at equal time, which is given, The short-time result, Λt L 2 , gives the dynamical scaling relation between the characteristic length in the y direction and time through the "growth exponent" β [18], At long times, Λt L 2 , the whole string can be treated as a random walker, with effective diffusion constant DΛ/L. These results, along with those from the string microscopics, fix the values of D and Λ. Comparison of V L (x, 0) with the exact result for the equal-time displacement variance gives D = 1/ √ 2. The mean rate of plaquette flips in equilibrium is γN f . Since each flip changes the mean vertical position by ±L −1 , the variance of the total shift is γtN f /L 2 in the longtime limit. Comparison of Eqs. (6) and (11)

B. Height correlation function
To calculate height correlations based on the coarse-grained string description, we write which treats the string as a step in the height plus a uniform gradient to preserve the boundary conditions. The height correlation function for q 0 can then be written, for a single string, as For small-enough density 1 2 θ, string contributions add inco- where the periodicity of V L (x, t) under x → x ± L has been used to shift the limits of integration. Asymptotic expressions for the correlations can be found in various limits. For the static correlations, G h (q, 0) ∝ θ/ω s (q), where ω s (q) = q 2 x + 1 8 q 4 y in the thermodynamic limit. For t L 2 , we find time dependence ; the proportionality constants are calculated exactly in Appendix C. A crossover from simple to stretched exponential therefore occurs at Λt ∼ q −1/β y , the time scale corresponding, according to Eq. (12), to the wavelength ∼ q −1 y . The stretching is the result of contributions from the continuum of modes of the string. These expressions, along with the full result found by numerical integration of Eq. (14), are compared with simulations in Fig. 2. We find close agreement, with no adjustable parameters, at large flux and small wavevector, and qualitative agreement at smaller φ x and larger q; see Appendix D.
The independent-string approximation should be valid for Λt θ −1/β , the time corresponding to a y displacement equal to the mean string separation. The stretched-exponential form therefore applies up to a time that diverges at θ = 0. For larger t, string interactions are important for the dynamics, and Eq. (15) is no longer valid. When Λt θ −1/β many strings contribute, their discreteness becomes unimportant, and so we expect a crossover to the Coulomb-phase behavior of Eq. (5). These crossovers can be understood via dynamical scaling theory, based on the critical point at θ = 0 [13,16], whose critical theory is that of hard-core bosons-or, equivalently, free fermions-in 1D. All critical exponents are therefore rational and follow from dimensional analysis.

C. Persistence
The persistence p(t) can similarly be understood in terms of the behavior of strings. At very short times, p(t) = e − f γt : plaquettes which are flippable at t = 0 will flip independently with rate γ. For times γt 1, each point x on a string performs a subdiffusive random walk, according to Eq. (11). As only plaquettes adjacent to strings are flippable, the mean persistence is equal to the probability that a plaquette is not reached by any string up to time t, and is therefore given by the survival probability for a stationary target in the presence of a density 1 2 θ of subdiffusive traps. Using the results of Ref. [19] relating the dynamic exponent to the persistence, we get This form ceases to apply for Λt ∼ θ −1/β , when string interactions become important, or for Λt ∼ L 2 , beyond which the long-time behavior in Eq. (11) applies and p(t) ∝ e −ct 1/2 . Our simulation results, Fig. 3, are in qualitative agreement with these arguments, showing an initial exponential followed by a stretched exponential, for all fluxes. The stretching exponent decreases continuously with φ x , approaching β = 1/4, in agreement with Eq. (16), as φ x → 1. At longer times, a faster decay, consistent with a single exponential, is observed.

IV. CONCLUSIONS
We have shown that the close-packed square-lattice dimer model, subject to local, plaquette-flip dynamics, displays emergent collective relaxation that is not anticipated by simple extensions of its static properties. Approximations to the dynamics based on free-energy gradients plus noise, such as Eq. (4), fail to capture the intrinsic heterogeneity: Due to the constrained nature of the system, motion is only allowed in the vicinity of strings, and relaxation is dominated by spatial fluctuations. In a sense, the noise that triggers rearrangements is not uniform in space and time; rather, its strength depends sensitively on the local configuration. Strings facilitate local rearrangements, dynamics is heterogeneous and collective (see Fig. 1), and relaxation functions are nonexponential.
This situation is reminiscent of glass-forming systems [20]: In a slowly relaxing material such as a glass former, "facilitation" indicates the fact that local relaxation can occur only near an already locally relaxing region [21]. Similarly, in the dimer model, plaquette moves are only possible in the vicinity of a string. This is the reason that the Langevin dynamics of Eq. (4) is not accurate for relaxation in regions where string density is low. The additive noise assumed in that approximation would allow rearrangements to occur anywhere in space. But if dynamics is facilitated, the noise that drives local rearrangements is not uniformly distributed, but rather concentrated near already mobile regions. A Langevin description, along the lines of Eq. (4), would therefore require a form of noise that is multiplicative and whose magnitude is strongly dependent on the local flippability.
Instead we have developed a string description of the dynamics, which directly incorporates the local nature of the relaxation. At low string density (high flux), we are able to make exact theoretical predictions for the correlations and persistence, which are confirmed by our simulations. We in fact find that much of the qualitative behavior is unchanged at smaller flux, where interactions between strings are certainly important. This indicates that the usefulness of the string picture, as well as the concept of facilitated and heterogeneous dynamics, extends well beyond the regime of high flux.
Besides their fundamental importance, our results are likely to be of relevance to spin ice, where closely analogous string excitations have been evidenced directly using neutron scattering [22], and where correlations with stretched-exponential decay have been noted [23]. Dynamical results for classical dimers are also relevant to the corresponding quantum dimer model at its Rokhsar-Kivelson point [24,25]. A string can be divided into four types of segment, illustrated in Fig. 4; the ensemble C s L of configurations for a single string is given by the set of ways in which these segments can be combined to produce a closed path of length L. One can write a transfer matrix such that T βα (k, µ) is nonzero only when a segment of type α can be followed by one of type β (labeled according to the order in Fig. 4). Each successive pair of segments is weighted by e −µ for every flippable plaquette it produces and by e −ik for every step in the positive y direction. We define the partition function Z o L as the weighted sum over all open paths of length L with any net vertical displacement y, where C o L denotes the ensemble of such paths and N f is the number of flippable plaquettes in the resulting dimer configuration. Summing over all sequences of path segments and (vertical) starting positions gives [? ] where the sum is over the set σ T(k,µ) of eigenvalues λ of T(k, µ). The allowed paths for a single string are those that return to their starting point after winding once around the system, and hence have net displacement y = 0. The partition function for such paths is For large L, the saddle-point approximation gives where is the largest eigenvalue (by magnitude) of T(k, µ) for any k.
The maxima are λ max (µ) = 1 + at the points k = 0 and ±π, and hence Setting µ = 0 gives the entropy of a single string, The mean number of flippable plaquettes in the presence of a single string is given by At flux φ = {1 − θ, 0}, the number of strings is N s = 1 2 Lθ. If the strings can be treated as approximately independent, as expected for sufficiently small θ, then the number of flippable plaquettes is simply N s N s f , and so the mean flippability in equilibrium is Numerical results, shown in Fig. 5, confirm Eq. (A16) in the limit of small θ and are in approximate agreement even for fairly large θ, suggesting that the independent-string picture is reasonable. (The logarithmic corrections modify the coefficient of θ, and are of the order of a few percent for L = 256.) The width distribution of a single string can be determined by a similar approach: The probability distribution for the net vertical displacement Y of a section of string with horizontal extent X is [? ] For large X, the ratio of sums can be found by expanding λ max (k, 0) in a Taylor expansion around its maxima, giving a pair of Gaussians of variance √ 2 X centred at k = 0 and ±π. Taking the Fourier transform, one finds that P s (X, Y) is given by a normal distribution of variance 1 √ 2 X when X and Y have the same parity, and vanishes otherwise (as required by the structure of a string).
If, for a single string traversing the system in the horizontal direction, we denote by y(x) the vertical position at horizontal position x, this result can be restated as follows: At length scales much larger than the lattice scale but smaller than the system size, 1 |x − x | L, the vertical displacement y(x) − y(x ) is normally distributed with zero mean and This result is confirmed by the small-displacement limit in Fig. 6.
Appendix B: Coarse-grained string picture The Edwards-Wilkinson equation, Eq. (7), has general solution (for t ≥ 0) (accounting for the periodicity in the x direction, but not in the y direction), where is the retarded propagator.
The two-time correlation function, within an ensemble of trajectories with fixed initial configuration, is therefore given by To calculate the equivalent correlation function in an equilibrium ensemble, · · · , we take both times to infinity while keeping their difference finite: In this limit, ∆(x, t 0 + t) = L −1 , and so the first term in Eq. (B3), which depends on the initial configuration, vanishes. The integrals in the second term can be performed to give where V L is given in Eq. (10).
For t = 0, the sum in Eq. (10) can be evaluated exactly, using for 0 ≤ θ ≤ 2π, to give giving  for Λt L 2 . As argued in the main text, Eq. (A16) can be used to fix Λ = 2( √ 2 − 1)γ. This result, including the value of Λ, is confirmed using MC results in Fig. 7.
At large t and x = 0, the first term in Eq. (10) dominates, and so for Λt L 2 .
Appendix C: Limiting forms of dynamical correlations In the thermodynamic limit, L → ∞, Eq. (10) can be replaced by where the Cauchy principal value is to be taken, and the correlation functions can be expressed in terms of For t = 0, one has V ∞ (x, 0) = D|x|, and so where ω s (q) = q 2 x + κ 4 y , with κ y = D 2 q y .
For small but nonzero t, consider the difference For Λt q −4 y , the integral over k can be replaced by 1 2 Λt × 2πδ(x), and so the result is Using Eq. (C3) gives For large time (but with Λt L 2 ), one can use the saddle-point approximation. The closest saddle points to the real line are at x = ±ix 0 , where and the resulting correlation function is C s (q, t) C s (q, 0) π 2 3/4 (Λt) 1/4 ω s (q) . (The precise condition for the validity of the saddle-point approximation also involves q x , but the quoted inequality is always sufficient, and also necessary except for very large q x /κ 2 y .)

Appendix D: Numerical results for dynamical correlations
Our analytical results for the correlations are based on a coarse-grained description of the strings, and so are expected to be quantitatively accurate only for small q. According to Eq. (C10), however, the stretched-exponential decay is visible only for Λt q −4 y , a time scale that grows rapidly as q y is decreased. The value q = { π 16 , π 4 } used in Fig. 2 of the main text is chosen to show the stretching most clearly on time scales accessible in the continuous-time MC simulations. (Because we have neglected the periodicity in the y direction, we also require Lq y 1.) Results for other values of q are shown in Fig. 8. As expected, the quantitative accuracy of the analytical results decreases as |q| is increased. Consistent with Eq. (C10), clear evidence of stretched-exponential decay is visible, as a decreased slope on a double-logarithmic scale, only for the larger values of q y . θ 0.2 (close, but not too close, to maximum flux). For larger θ, the density of strings is sufficiently high that their interactions become important, and the picture of independent strings breaks down. When θ L −1 , finite-size effects become important.