Anomalous relaxation of density waves in a ring-exchange system

We present the analysis of the slowing down exhibited by stochastic dynamics of a ring-exchange model on a square lattice, by means of numerical simulations. We find the preservation of coarse-grained memory of initial state of density-wave types for unexpectedly long times. This behavior is inconsistent with the prediction from a low frequency continuum theory developed by assuming a mean-field solution. Through a detailed analysis of correlation functions of the dynamically active regions, we exhibit an unconventional transient long ranged structure formation in a direction which is featureless for the initial condition, and argue that its slow melting plays a crucial role in the slowing-down mechanism. We expect our results to be relevant also for the dynamics of quantum ring-exchange dynamics of hard-core bosons and more generally for dipole moment conserving models


I. INTRODUCTION
The field of dynamics in isolated quantum systems has recently received an increasing amount of attention thanks to the discoveries of a plethora of interesting nonequilibrium behaviors [1][2][3], and of versatile experimental platforms to realize the same [4,5]. These studies have been partially motivated by the desire to achieve the protection of quantum information from scrambling caused by Hamiltonian dynamics or environmental noise. This has lead to the rapid development of the field of many body localization [6][7][8][9], which relies on strong disorder to provide a safeguard against scrambling many body dynamics, and the nascent field of Hilbert space fragmentation [10][11][12][13][14][15], which results from the highly constrained configuration space of Hamiltonians with a large number of strong local constraints and/or highly frustrated interactions [16,17].
A milder version of the total arresting of dynamics generated by the phenomena mentioned above is realized by systems which approach equilibrium in a manner which is qualitatively different from standard diffusion. An example of this has recently been explored for the spin-1/2 Heisenberg chain as well as its classical version at intermediate energies where evidence of super-diffusion has recently been seen both theoretically [18,19] and experimentally [20], leading to connections to surface growth dynamics studied by Kardar-Parisi-Zhang [21,22]. An approach to equilibrium which is slower than that expected from diffusion has also been realized in systems which conserve higher moments (such as dipolar and octupolar) of the spin configuration [23][24][25][26][27]. For two-and higher-dimensional systems, this has also lead to the realization of exotic fractonic phases [28][29][30]. Remarkable advances have also been made on probing experimental realizations of Hilbert space fragmentation and/or higher-moment conservation and associated subdiffusion [31][32][33].
While progress has been made on the analytical description and experimental detection of the phenomena mentioned above, there has been a dire need for numerical simulations on microscopic models to lend support to many of the predictions. This is in general a difficult task as powerful methods to simulate large-scale quantum dynamics on the needed long-time scales are relatively few and capable of handling only specific regimes [34,35]. It was found insightful [23,36] to adopt methods from stochastic classical mechanics in the framework of cellular automata [23,37,38], which ignore part of the quantum phase fluctuations and have been able to successfully describe the long-term dynamical behavior of strongly interacting quantum systems. This intuition arises from the expectation that for generic systems with sufficiently large Hilbert spaces and for times long compared to the microscopic energy scales of the Hamiltonian, the dynamics does not show quantum coherence, thus reducing to a classical dynamics problem. As mentioned above, wellknown exceptions to this exist, but as simulation of exact quantum dynamics is out of reach using current methods, a study of the classical equivalent becomes of particular interest. This can also serve as a natural starting point to understand the complete quantum dynamics. The language of cellular automata also lends itself naturally to a hydrodynamic treatment, which identifies an equivalence between quantum many-body dynamics at late times and classical transport of globally conserved quantities [38][39][40]. Slow thermalization is often expected in integrable systems described by generalized hydrodynamics [41,42]. Quantum equivalents of cellular automata, which may be expected to capture the dynamics more accurately, have also been found to share important characteristics of integrable systems [43,44].
Following up on the studies of constrained systems, we consider in this work the case of a simple hard-core bosonic model living on a square lattice, undergoing ringexchange dynamics. This model has already been studied from the perspective of cuprates, as they serve as promising candidates to realize high-Tc superconductivity [45][46][47]. Although traditionally most studies have focused on the possible exotic ground state features of this arXiv:2211.16788v2 [quant-ph] 15 Apr 2023 model [48,49], some recent works [23,50] have considered the constrained dynamics generated by the ring-exchange mechanism, including starting from random configurations [23]. However, the relevance of fragmentation to generic low-momentum states which are only described by macroscopic patterns has not been investigated.
In this work we address this question using a classical approach based on stochastic dynamics. We find that structured initial configurations in the form of a boson density wave retain their coarse grained structure for a time which grows as a tunable power of the wavelength, with an eventual melting which is approximately described by a continuum model derived from a simple Taylor expansion. We study the detailed structure of the melting process via spatial correlation functions and find that the dynamics proceed through the development of strongly correlated large active regions which merge and destroy the initial modulated pattern.
The detailed plan of the paper is as follows. In Sec.II, we present the model, discuss the quantities conserved under the dynamics, and elucidate the general profile of stripe-like configurations which are perfectly frozen under the dynamics. The bulk of this section discusses the effects of small perturbations on these exactly frozen patterns, and the preservation of the memory of the initial state to infinite times as illustrated by simulation of exact quantum and stochastic classical dynamics on a small system size. The close agreement between exact quantum dynamics and its stochastic classical equivalent seen in this section motivates our approximation and we focus purely on the classical system for the following sections.
Sec.III recalls the expected continuum field theory based on Taylor expansions by first considering the simpler case of correlated random walkers, along with numerical checks of the equations developed. This is followed by a treatment of the hard-core model using a mean-field assumption.
For Sec.IV, we move to more general configurations which take the form of boson density waves, and show the persistence of the memory of the initial state for unexpectedly long periods of time. We also compare the prediction of the continuum field theory developed in Sec.III with the initial dynamics.
We follow up in Sec.V with a detailed analysis of the evolution of the dynamical active regions, and present a phenomenological picture of the mechanism leading up to the melting of the initial density wave configuration.
We summarize our results in Sec.VI and present possibilities for future follow ups via direct simulation of the quantum many body dynamics.
In a following appendix A, we briefly discuss the longtime momentum space profile of correlation functions, show that it is consistent with a mean-field treatment obtained in a previous work, and that it cannot be a basis of the anomalous scaling observed in the present work (as could have been deduced from a recent analysis [51] of a related model).

II. RING-EXCHANGE BOSON DYNAMICS, CONSERVED QUANTITIES AND FROZEN PATTERNS
We consider a system of hard-core bosons living on a square lattice, which evolve stochastically in time using only ring-exchange dynamics where bosons hop by pairs around a 4-sites plaquette of the lattice if and only if a single diagonal of the plaquette is occupied by two bosons ("flippable" plaquettes). This only allowed dynamics is shown in Fig. 1. This dynamic rule trivially conserves the total particle number as well as the number of particles in each individual column and in each individual row [23,[48][49][50]. For the rest of this work, we restrict ourselves to half-filling, i.e. L 2 /2 sites occupied by bosons, on a periodic L × L lattice, and to the sector where each column and each row has exactly L/2 occupied sites. One expects it to be the sector with the largest number of configurations, as it is maximally symmetric. The total number of configurations in this sector can be computed on large lattices using combinatorial techniques [52].
We find that the conserved quantities discussed above do not in themselves completely describe the dynamically connected sets of configurations. This can be seen by considering a "perfect" stripe configuration as shown in Fig. 1, where an alternating pattern of filled and empty sites does not leave a single flippable plaquette, making the configuration frozen under ring-exchange dynamics. By varying the widths and locations of the filled and empty stripes, we can create many similar frozen configurations. One may also naively expect that this restriction on the number of accessible configurations extends to the case where we do not have perfect stripes. To see this, we consider a configuration generated by interchanging the two neighboring diagonals on the edge between a filled and an empty region. This creates a diagonal made of flippable plaquettes, where the influence of this active region may be expected to only extend to a few lattice spacings around this diagonal. Note that similar arguments are expected to hold more general dynamics which conserve the dipole moment. To understand this intuitively, we can consider larger plaquette dynamics, for example 2×1 or 2×2 plaquettes; once again the perfect stripe configuration is frozen and spatial perturbations around it can be expected to at most soften the boundary dynamically, but still leave the stripes intact provided they are wide enough. A simple example of the effect of the 1×1 plaquette dynamics on a boundary between a Néel like (highly active) and a fully filled (inactive) region is shown in Fig. 2. Here one can see that although it is easy to transfer a hole from diagonal d 1 to d 2 , doing the same for d 3 would require holes at the circled locations, and thus it is not possible to propagate our excitation into the inactive region without sourcing another flippable plaquette from within the active region.
To test this intuition and illustrate this effect, we exactly enumerate all the configurations connected to an initial condition of the form described above (slight ran- dom perturbation to the perfect stripe configuration) for a 16 × 16 lattice (shown in Fig. 1). We find that the total number of configurations, which we call N c which belong to this "fragment" is 27, 990. To show that they retain some structure of the representative which we have chosen, we compute an overlap as defined below, where n c/seed x,y = 0, 1 is the number of bosons at site of coordinates (x, y) in the (c-th/initial) configuration. If we use all possible configurations (in which case N c would be the total Hilbert space size), we would expectŌ = 0 due to the symmetry under n → 1 − n. For the restricted Hilbert space belonging to the fragment being considered we find0 = 0.844003..., showing that for a dynamic simulation restricted to this sector, the initial and late time states would still retain significant over- Illustration of the melting of a Néel-filled boundary: Asterisks mark flippable plaquettes, and the transformation from left to right takes place by an exchange around the plaquette marked by the asterisk in bold type. di mark diagonals around the domain boundary.
lap. To confirm this, we also run a stochastic classical automaton (SCA) simulation, where at each time step, we propose L 2 random plaquette flips, and if the chosen plaquette is flippable, we flip it with probability 1/2. The resulting overlap with the initial state is shown in Fig. 1, and we see that it quickly approaches the value expected from exact enumeration, and retains it indefinitely. To study the accuracy with which the SCA reproduces the exact quantum dynamics, we also perform an exact diagonalization for this Hilbert space fragment, and compute ψ(t)|O|ψ(t) for |ψ(t) = e −iHt |ψ 0 using where the initial condition is the same occupation-basis state as we initialized the SCA with and O now denotes an operator which is diagonal in the occupation-basis and measures the overlap with the initial state. This operator can be generated directly from the expression above for the classical case by promoting n c x,y to the number operator, while retaining the integer status of n seed x,y . As shown in Fig. 1, we find that the overlap in the quantum dynamics closely traces the SCA during the initial decay away from an overlap of 1.
In the following sections, we quantify the relaxation to equilibrium using the Fourier transform n k (t) with frequency k along thex +ŷ direction, at time t. To this end, we measure the Fourier ratio R(t) = n k (t)/n k (0), and make a similar comparison as done for the spatial overlap above. This is presented in Fig. 1 as well, and once again we find a close agreement between the exact quantum dynamics and it's classical equivalent in the way they approach the equilibrium value of R(t) within the sector of interest. This result suggests that the quantum dynamics within sectors matches the classical automaton upon coarse graining past time scales of O(1), suggesting that quantum phases do not play a substantial role in the aspects of dynamics which we want to study, and that a classical automaton approach could be sufficient to study the effect of the kinetic constraints on large scale features.
Following the arguments above and noting that exact quantum dynamics for larger sizes is not possible with current computational capabilities, we directly study the long time behavior of our stochastic dynamics simulations to gather information about large system sizes in the following sections.

III. EXPECTED HYDRODYNAMIC DESCRIPTION
Here we discuss the coarse grained description in continuous space, of the microscopic dynamics which we have introduced in the previous section. We begin by relaxing the hard-core constraint and considering the limit of average number of particles per site being 1. This allows us to reduce the problem to non-interacting correlated random walkers, and leads to an analogue of the diffusion equation which encodes vanilla sub-diffusion. We show numerical evidence for the validity of the same. The expected continuum dynamical behavior of this type of ring-exchange in 2d was first presented in Ref. [23], but we recall it for completeness as well as to understand how it should be perturbed to take into account the hard-core nature. Next we move to an equivalent model of a realvalued field on the lattice which makes the connection to the hydrodynamic limit more transparent, and we again provide support with numerical simulations.
Lastly, we return to the hard-core model presented in the previous section, and show that the continuum theory describing the model must include non-linear terms in addition to the sub-diffusive term, and we show the possibility for a quantitative change in the behavior of the system due to the non-linearities.

A. Large particle number limit
We relax the hard-core constraint of the stripe configurations described in the previous section and assume the pattern to exist over a featureless background of an average density of n d 1 particles per site. The dynamics can now be understood in terms of correlated random walkers in the following way. First, we label each particle in the system as an independent walker. A move is defined now as first picking a walker (called walker a) at random, picking one of its four next nearest neighboring sites with probability 1/4, and moving one of the walkers on the chosen site (called walker b) in tandem with walker a in a ring exchange type manner. Due to the large number of walkers per site, we expect to always be able to find walker b.
To write down the number of particles n at a particular position (x, y) at time t + 1 as a function of the values at time t, we must consider all processes which can change n x,y . These processes are (1) choosing a walker at (x, y) and moving it away (∝ n x,y ); (2) choosing a walker at one of the four nearest neighbor sites and moving it or its partner to (x, y) (∝ 1/2(n x±1,y ) or 1/2(n x,y±1 )); (3) choosing one of the walkers at one of the next nearest neighbors and moving it in tandem with a walker at (x, y), thus reducing the number of walkers at (x, y) by one (∝ 1/4(n x±1,y±1 )). Putting these terms together, we can write the change in n x,y , which is an whole number, from time t to time t + 1 as (note that all terms in the expression below are at time t, and we do not explicitly mention it for ease of representation) This equation was obtained using similar arguments in Ref. [23]. If the system is initialized over a background of n d particles per site using the stripe configuration described in the previous section (shown in Fig. 1), n (x,y) can be seen as a step function switching periodically between n d and n d + 1. We would naively expect that the correlated dynamics discussed above would quickly eliminate the sharp boundaries of the n (x,y) texture and lead to a smooth function once we average over many realizations of the stochastic dynamics.
For a function which varies slowly as a function of (x, y) (note that this implies that the stripes in the initial condition should be wide compared to lattice spacing), we can perform a Taylor expansion of the expression above. We find that all terms to fourth order cancel, and the only term at fourth order yields where c = 1 if following the treatment above. For further convenient reporting of the wave-vector k in units of π, and taking into account in addition a factor of 4 coming from the acceptance probabilities in our numerical implementation, we consider instead a different normalization with c = π 4 /4. Rescaling the x-axis as in Fig. 3, this allows to recover a match to e −x for the fit in Fig. 3. For the rest of this manuscript, we maintain this convention for c.
As the stripe initial condition is a periodic square wave in (x + y), it is convenient to rewrite the above equation in the Fourier basis, and consider only the lowest harmonic (largest wavelength) of the square wave transform. This reduces the dynamical equation to ∂ t n kx,ky = −ck 2 x k 2 y n kx,ky . For diagonal stripes, k x = k y = k, and n k can be exactly reduced to n k (0) exp (−ck 4 t), where n k (0) is the value at t = 0. We can now numerically verify this behavior by calculating the Fourier ratio R k (t) = n k (t)/n k (0), and looking for a data collapse onto a single exponential for various values of k. We find that n d = 3 provides a sufficiently large background for a good data collapse, and show the same for a 128 × 128 lattice for various stripe widths (encoded in k) in Fig. 3 and averaged over 80 realizations of the stochastic dynamics for each k. For smaller values of n d , we find an increasing discrepancy between different values of k, with the divergence growing with decreasing n d .

B. Discretization of continuum equation
Here we study a discretization of Eq. (3) on the lattice by defining the real valued field n x,y , and attempt to recover the limit of the hard-core model. From the analysis of the correlated random walkers above, we expect that a ring exchange dynamic on a plaquette should lead to the ∂ 2 x ∂ 2 y form. To this end, we define the "activity" of a plaquette whose left bottom site is (x, y) to be a = [n x,y+1 − n x,y ] − [n x+1,y+1 + n x+1,y ]. This is evaluated at time t, and the fields living on the plaquette are updated at time t + 1 by adding a to sites (x, y) and (x + 1, y + 1) and subtracting a from the other two sites. The constant is included to ensure that the field remains ∈ [−1, 1], and plays only a quantitative role in the scaling study. To deduce the dynamic rule for the activity at a single lattice site, we must consider the four neighboring plaquettes around it which can affect the field on the chosen site via plaquette updates. By summing a for the four plaquettes with equal weight, and carrying out a careful grouping of terms, we see that the equation reduces exactly to the dynamical equation discussed in the previous subsection. Using a Taylor series expansion once again leads to a dynamical equation of the form given in Eq. (3). We attempt a data collapse for the decay of R(t) as defined in the previous subsection for stripe configurations using the dynamic rule mentioned above. As shown in Fig. 4, we once again find a satisfactory data collapse to a single exponential for a large range of k values.
Before we turn to the case of the hard-core model, we must note that the dynamics described in this subsection differ in one crucial way from the hard-core model. We can see this by considering all the plaquette configurations which yield a non-zero value of a (these are listed in Fig. 5), and observing that only the two completely "flippable" (or equivalently with the largest magnitude of a) contribute to dynamics in the hard-core case. This is not the case if one considers the dynamical rule used here and in the previous subsection (all plaquettes with all values of a are updated). We argue in the next subsection that this leads to a strong violation of Eq. (3).

C. Continuum theory under the hard-core constraint
Now we turn to the model described in the previous section, and restrict particle number to at most one per site, and dynamics only to proceed via exactly flippable plaquettes. Note that now we cannot define an updated state using just the definition of a as in the previous subsection. We require a function which evaluates to unity only for the two flippable plaquette configurations and zero for all others. For the purpose of this subsection, it is more convenient to set n (x,y) = 1 for an occupied site and −1 for an unoccupied site. It is not apparent if there is a unique function which achieves this, and one of the simplest functions which we were able to find on the plaquette to satisfy these constraints is where Σ p n = n x,y + n x,y+1 + n x+1,y+1 + n x+1,y is the sum of all n belonging to the plaquette. A careful consideration of the expression above reveals that it generates a value of ±1 for the flippable plaquettes shown in Fig. 5a, while returning zero for all other configurations (including those in Fig. 5b), thus satisfying our requirements for a hard-core ring exchange. The second term in the above expression is formed by noticing that the only configurations which violated the assignment of values we desired had a difference in the types of pair arrangements on opposite edges. Before performing approximations on this expression to derive a continuum theory, we find that it is convenient to expand the product of the last two bracketed terms discussed above as (2n x+1,y n x,y+1 − 2n x,y n x+1,y+1 ), where we have used n 2 x,y = 1 for all (x, y).
For the evolution ∆ t n x,y of density of a single site, we must once again consider the four plaquettes in which it participates. The resulting expression has a term identical to the dynamic equation in the previous subsections, but has an important addition in the form of the sum of 1 2 (Σ p n)(n x ,y n x,y − n x,y n x ,y ) over the four plaquettes.
To get the continuum limit as done previously, we assume initial conditions with several flippable plaquettes and a smooth density profile which varies slowly spatially. An average over stochastic dynamics and initial conditions thus allows us to replace the terms linear in n with the differential form −∂ 2 x ∂ 2 y g, where g is a real valued field living in continuous space-time aiming at replacing n (we explicitly distinguish it from n for this subsection due to the correlations possibly generated by the products of n).
To understand the behavior of the averaged value of terms such as n x+1,y n x,y+1 , an important assumption about the correlations has to be made. In the following of this subsection, we work in a mean-field picture ignoring correlation effects and assume that such terms can be rewritten as the product of g at the two points. In doing so, we will be able to derive the complete dynamics only in terms of the field g. This mean-field like assumption is true for the initial condition we work with, as it is drawn from an uncorrelated ensemble, but it is not a priori evident if the dynamics maintains this uncorrelated nature or rapidly builds up correlations. The non-linear term f is expressed as ( x,y n)(n x+1,y n x,y+1 − n x,y n x+1,y+1 ) where the subscript below the sum indicates the plaquette over which the sum is performed indexed by its left bottom site. Once again we expand the products and ensure that we replace all occurrences of n 2 x,y with unity for all (x, y). The equation above is thus reduced to a linear combination of single body and three body terms. This can be further reduced by performing the mean field decoupling n x,y n x ,y n x ,y = n x,y n x ,y n x ,y = g x,y g x ,y g x ,y , followed by a Taylor expansion around the relevant site in the derivatives of g. As we began our analysis by considering a sum over the four plaquettes surrounding the site (x, y), symmetry restrictions apply to the terms generated by the Taylor series, which require that only terms with non-zero coefficients have an even number of derivatives in x and y, and are symmetric with respect to the x → y transformation. Using this constraint and after some algebra, we find that the only surviving terms arise at fourth order in the derivatives and that the complete dynamical equation reduces to the following expression (ignoring a global factor of 1/4), given by The presence of non-linearities, derived even under a crude mean-field approximation, suggests that at leading order in the dynamics, the hard core constraint indeed plays an important role, and may invalidate the expectation that the coarse grained dynamics are equivalent to those of vanilla sub-diffusion.

IV. LONG TIME PERSISTENCE AND EVENTUAL MELT OF APPROXIMATE STRIPE CONFIGURATIONS
We have seen in the Sec. II that configurations which can be viewed as small perturbations around a perfect stripe configuration may maintain the memory of the initial state indefinitely. As these configurations are highly specific, it would seem unrealistic to choose one of these as the initial state for the dynamics of large systems. This motivated us to study "approximate" stripe patterns, which are chosen to be boson density waves with a wave-vector k = (k x , k y ).

A. Initial state preparation
To prepare such configurations, we first generate a target distribution using the function where A takes continuous values between 0 and 1. We cannot generate a configuration which has exactly this pattern as a boson configuration can take only ±1 (filled/empty) at each site. Thus we perform a Monte Carlo simulated annealing starting from a charge-density wave state with k t = k tê x + k tê y using an energy defined as E = l (D l − Lf k t ( r)) 2 , where the sum is over all diagonals and D l = (i,j)∈l σ i,j , is the occupancy in diagonal number l. This favors an exact match between the current and target configurations. By tuning the inverse temperature for the annealing from β = 0.01 to β = 20.48 in a geometric progression consisting of multiplying by 2 every 10L 2 steps, we achieve a stripe configuration which has a single Fourier component at the target k t . We ensure that the proposed configuration changes are long-ranged and respect the conservation laws discussed in the previous section, thus staying in the sector where every column and every row has exactly L/2 bosons. An example of a configuration created by this procedure is shown in Fig. 6a. To check the effectiveness of the annealing procedure, we record the final energy and find it to be O(L), implying that each D l is within an O(1) value of its target value, which is expected as D l is an integer and the target value is in general a real number, and not necessarily close to an integer value. We find that this procedure generates a sharp peak for the desired k t amid a weak background which decays with increasing system size. This can be seen in Fig. 6b, where we plot the Fourier component as a function of the distance from the target peak location in Fourier space for various sizes. Due to the periodic boundary conditions in Fourier space, given by (k x , k y ) ↔ (1 − k x , 1 − k y ), we consider the shortest distance to the expected peak over the naive distance for open boundary conditions. Note that such a configuration still contains a large number of flippable plaquettes due to the smooth nature of the target pattern.

B. Comparing with continuum field theory
In the context of the continuum field theory discussed in Sec.IIIC, the field g at time zero would now simply be equal to A sin( k · r) with k = k t . It is essential to numerically verify the validity of the assumption used in developing our continuum theory, namely that spatial correlations do not play a role in the initial dynamics.
Using this as the initial source field in Eq. (5), we find that ∂ t g reduces to This immediately indicates that we have lost linearity (which would imply no dependence on amplitude other than a global proportionality), and that the theory is no longer exactly separable in Fourier space. Note here that the terms proportional to A 2 within the brackets oppose the decay generated by the simple sub-diffusion, leading to a possibility of further slowing down the dynamics. The above equation also suggests that for A 1, we should expect to recover vanilla sub-diffusive dynamics. For k · r ∈ (0, π), the sign is controlled by the expression within the brackets, and we can see that a growth of the function can be achieved for sin( k · r) < ((8A 2 − 1)/(9A 2 )) 1/2 , provided 8A 2 > 1 (that is large enough amplitude). This condition is satisfied for k · r in the vicinity of 0 and π. This effect is rather nonintuitive as it implies that the local density tends away (towards ±1) from its equilibrium value (0) under nonequilibrium dynamics for a tunable range of r, leading to a local reduction of entropy as the number of states available locally reduces if we require their average to be closer to the most extreme values which it can take. To ascertain the extent to which our microscopic dynamics is consistent with the continuum field theory developed under the mean field assumption, we can now compare the averaged value of n x,y against a numerical evolution of Eq. (5). Looking specifically for the feature described above, we plot B x (t) = n x,x (t) / n x,x (0) in Fig. 7a as a function of x for the minimal period. Our initial condition already sets n x,x (0) = A sin(2kx). We see a qualitative match to the prediction from the continuum equation (in the sense of the sign of the difference to unity of the ratio is well captured), but a quantitative disagreement (in the amplitude of this difference) potentially due to the build up of correlation effects beyond mean-field.
Before we move to the study of the microscopic dynamics, it is worth considering the time scales over which the initial pattern melts under the continuum dynamics, as done in Sec. II. We show in Fig. 7b that the dynamics is still consistent with a single scaling variable, given by k 4 t, although the relaxation deviates from the single exponential seen in the simpler cases considered in Sec. II. This implies that the additional non-linear terms in the continuum dynamics do not alter the scaling of space-time, and this can be intuitively understood by observing that all non-linear terms have the same order in derivatives.

C. Exact numerical evolution: Beginning of melting
We now study the dynamical behavior starting from a configuration generated by the method discussed above. To carry out an analysis of the evolution of the coarsegrained structure of the initial configuration, we choose to study the decay of the dominant Fourier component via the already defined ratio R(t) = n k (t) n k (0) . To get an intuition about the effect of lattice spacing, we first consider the decay of R(t) over a range of decreasing values of k, starting from k = 1/2. As shown in Fig. 8a, for A = 0.2, there are clear deviations away from a single exponential for large values of k, whereas intermediate values of k appear to agree partially, and for small values of k, we once again deviate from the expected k −4 scaling. To study the initial stages of the decay of the wave pattern, we can look at R(t) in the regime where it is close to unity. Here we find a collapse consistent with a time scaling with k −4 for a large range of k, includ- ing large k values (as shown in Fig. 8b). This is quite surprising, because the regime of short-time scales and large-wave vectors is not the one where we expect the hydrodynamic prediction of Sec. III to hold. We have no simple explanation for this observation.
For larger values of A, we observe a completely different behavior. We find that R(t) ≈ 1 for a non-trivial amount of time after the initialization of the dynamics, an illustration of this behavior of R(t) is given in Fig. 9 for A = 0.5. We define the beginning of the melting process by the first time t 0 at which R(t) crosses 0.99. This threshold is chosen arbitrarily and a different threshold does not change the result qualitatively (as shown using a threshold of 0.9 in the inset of Fig. 9b). We study this for a few values of k, and for 20 realizations of the initial condition for each k. In addition to this, for each realization, we run sufficient number of realization of the stochastic dynamics to ensure that we get a good esti-mate for t 0 . We have taken systems of linear size in the range of 150 − 250, as the property of self-averaging allows us to narrow the spread in the values of t 0 . We find a strong dependence of the averaged t 0 on the inverse wave-vector 1/k. A linear fit on a log-log scale of t 0 vs 1/k reveals various power law regimes (see Fig. 9b) that depend on the amplitude A, the most extreme of which is achieved for A = 1, where the dependence is t 0 ∼ (1/k) 18.1(3) . Such a strong dependence on the initial pattern suggests that the mechanism for melting is initiated by a coordinated rearrangement of bosons which is favorable for dynamics. Note that for all values of A, there exists a window of 1/2 > k > k 0 (A), where we find good agreement with (1/k) 4 , as predicted from the continuum theory, with k 0 (A) reducing with reducing A.

D. Exact numerical evolution: Post-melt scaling
The continuum theory developed in Sec. III, which can be valid only for small wave-vectors and long times, cannot describe the results in Fig. 9b and the period for which we observed R(t) ≈ 1, and this suggests that it is not the appropriate theory to understand this "prethermal" behavior.
We can check however if the continuum prediction is upheld after the melting process with a simple scaling collapse. We first define a time t = 0, which is taken to be the well after the beginning of the melting process, as the time at which R(t) drops below a threshold of 0.01 (chosen just for the convenience of analysis, once again we checked that this value only plays a quantitative role). We can now attempt a parameter free scaling collapse of R(t ) and would expect to get a linear profile on a log-linear plot for a range of k by scaling the t → k 4 t . We present this analysis in Fig. 10 for A = 1.0 and find a reasonable agreement with the expectation developed above, even if not entirely satisfactory. To further check the applicability of the continuum theory, we also consider ring exchange dynamics over 2 × 1 plaquettes, which we expect to the same Eq. (3), but have much short pre-thermal times (as naturally be expected from a longer-range type of exchange). This allows to probe a larger range of wave-vectors k values, as presented in Fig. 10 where the adherence to a k 4 t is clearly improved.

V. PATTERN OF DYNAMICAL ACTIVITY IN PRETHERMAL AND MELTING PROCESS
The previous section shows that a mean field treatment is unable to recreate the slow dynamics seen numerically, thus hinting at the presence of correlations which are neglected at a mean field level. To understand the large-scale mechanisms involved in the prolongation of the prethermal lifetime and the onset of melting, we consider the build up and correlations of flippable plaquettes. This serves as a proxy for identifying dynamically active regions and their evolution. We study real-space flippability correlations through the normalized connected correlator C(x, y, t) = P (0, 0, t)P (x, y, t) − P (0, 0, t) P (x, y, t) where P (x, y, t) = 1 if the plaquette whose left bottom site is (x, y) is flippable and zero otherwise, and P (t) is the spatially averaged density of flippable plaquettes. Large (low) values of C(x, y) at a given time will thus indicate strong (weak) correlations of flippability with the initial point. Recall that we have chosen starting configurations from annealing to a potential which is only a function of the tilted coordinate x + y. This means that we have a freedom of choosing the origin for our correlation function at any x − y for fixed x + y. We take advantage of this symmetry by averaging over all equivalent positions of the origin. For the x + y position of the origin we use our potential function as defined in Eq. (6) and set our origin to be the point satisfying k · r = π/4. We make this choice as it lies at the threshold between highly active and inactive regions, defined by n = 1/2 (medium boson density, high density of flippable plaquettes) and n = 1 (high boson density, low density of flippable plaquettes).

A. Density wave patterns
We begin by considering the initial conditions studied in the previous sections. To ensure that we are able to observe large scale (slowly varying spatial) features in the correlation function, we choose a system size of 80 × 80 and k = 0.05 at A = 0.6 (a particular realization of this is shown in Fig. 11). To gather high quality statistics we average over 80 realizations for the initial condition and for each realization we run 128 realizations of the stochastic dynamics. The time t 0 , which denotes the beginning of the melting process, is ≈ 10 5 or 2 17 for k = 0.05, as seen in Fig. 9. We study C(x, y, t) up to time scales of 2 21 and find that strong anti-correlations develop at short range for the density of flippable plaquettes at early times, and sustain until the intermediate stages of melting. This is shown in Fig. 11, where snapshots of C(x, y, t) at t = 2 17 → 2 21 are presented.
This profile shows the development of active regions (bright) surrounded by inactive regions (dark), and suggests that a mechanism of excluded dynamical regions may be the source of dynamical slowing down. Note also that Fig. 11 suggests a bias of the dynamics towards the x and y axes, which may be expected from the continuum theory as well, due to the lack of radial symmetry. We see that even before the melting time, weak correlations are already built up along the x and y axes, but in a manner which is anti-correlated and modulated with the approximate stripe width. We find that this pattern begins to appear at times as short as t = 2 4 = 16, but we present data at the last possible time value before melting which we have recorded as the pattern is visible with significantly more clarity.
Another important aspect of the correlation pattern to note here is the dependence on the x − y coordinate ("perpendicular" to the initial wave-pattern), which is not built in in the initial condition as the energy used in the annealing process used to generate the configuration has no x − y dependence. This behavior cannot also be expected from the continuum theory since we start of with a single wave like configuration with k x = k y .
We find that an important condition for the existence of such correlation patterns is that they are specific values of A where the dynamics is not described by simple subdiffusion. For example, we studied the case of k = 0.05 at A = 0.2, which shows conventional dynamics (as far as we can see from Fig. 9), and did not find any non-trivial correlations down to a precision of 10 −4 in C(x, y, t). This is as expected from a mean-field treatment, as it forbids correlations by definition.
The evolution of the correlation landscape for the duration of the melting process presented in Fig. 11 suggests that it proceeds through a merging of dynamically active regions. This will be discussed in more detail for step-like initial conditions, where a complete melting process can be observed with higher clarity and on longer times.

B. Maximally active square wave initial conditions
To gain a better understanding of the melting process, we begin with a more artificial initial condition chosen to have regions of flippable plaquettes with maximal and minimal density. We consider alternating patterns such as the one shown in Fig. 12, with W (necessarily odd) consecutive diagonals filled, followed by W + 1 in a perfect staggered ("Néel") pattern, followed by W empty diagonals, and finally by another W + 1 in a Néel pattern to close out one period. These constraints are chosen to ensure that a filled region is bounded by empty diagonals and vice versa. The lattice size is then N (2W +2(W +1)), where N is the number of periods. Although these configurations show melting times (quantified by the decay of the appropriate Fourier component) which are faster than those studied above, they show a scaling fore pre-melting which is slower than the conventional k −4 , as seen in Fig.11. We find that the flippability correlator C(x, y, t) for these unconventional initial conditions shows an extension of the anti-correlated pattern in the x − y direction with a periodic modulation across distances large compared to lattice spacing. This is seen clearly for W = 13 with N = 2 at time of t = 2 23 = 8 × 10 6 in Fig. 12.
An important characteristic of the pattern described above is the value of the wave-vector associated to this pattern. We find that this is determined by the approximate stripe width, i.e. the wavelength 1/k. Our data suggest that this occurs due to a development of the correlation along the x and y axes, which is limited by the size of active regions (high density of flippable plaquettes). Growth of correlation within an active region is mediated by alternating patterns of correlation and anticorrelation, which can be thought of as being generated by reflections off of the boundary between active and inactive regions, as shown in Fig. 12. Repeated reflections create the observed periodic pattern, thus linking the periodicity and the width of the initial stripes. This is clearly seen in the approximate periodicities of the pattern for W = 13 and 27 in Fig. 12 where the profiles look similar even though the stripe width is changed by a factor of 2. Strong signatures of the periodic pattern for stripe widths (or alternatively wavelengths) of 13 and 27 at approximately the same time indicates the emergence of this pattern at a time scale which only weakly depends on the wave-vector k, and that it develops at a time scale parametrically much smaller than the melting time. Note that an extrapolation based on Fig.11b indicates the beginning of melting to occur earliest at To develop an estimate of the relevance of this patterning to the process of melting, we first remove the normalization of P (t) which we have absorbed into the definition and instead look at C(x, y, t) scaled by P di (t) P dj (t) , where P di (t) is the average number of flippable plaquettes in diagonal number i at time t. This makes the correlation function effectively lie between −1 and 1, with either limit describing a saturation to the largest (smallest) values. In particular, a value of −1 implies that the density of flippable plaquettes is zero, leading to a complete arrest of the dynamics. This version of C(x, y, t) is plotted in Fig. 12e and we see clearly that close to the frontier of dynamical activity, we do get a complete anti-correlation around the reference region. This suggests that the dynamical activity is likely strongly controlled by the unconventional pattern in the x − y direction. Now we turn to an investigation of the role of the periodic patterns observed above in the melting process. We can study this by choosing a small enough W so that the melting process is captured within the time scale of 10 9 which we can simulate. We find that for W = 7, the ratio R(t) begins to deviate from unity at a time t 0 ∼ 2 21 = 2 × 10 6 , and reaches a value of 0.01 by t b ∼ 2 24 = 1.6 × 10 7 (data plotted in Fig. 13). In this interval, we find at time t = 2 20 a formation of stable periodic structures in the x − y direction, and see a connection of the dynamically active regions for t = 2 22 . Finally by t = 2 24 , the correlation pattern has relaxed into a single wave in the x + y direction, restoring the symmetry in the x − y direction expected from the continuum prediction. The role of correlations is expected to be insignificant past this time point.

VI. SUMMARY AND OUTLOOK
We studied a system of hard-core bosons on a square lattice evolving under classical stochastic dynamics using ring exchanges. We find that boson density waves motivated by the patterns present in the Hilbert space fragmentation of this model approach equilibrium on an extremely long time scale, which diverges with the inverse wave-vector as k −α for small k, where α depends on the amplitude A of the initial density wave configuration. For decreasing values of A, we find increasingly larger windows in the large k regime where we observe a scaling of k −4 . Our numerical studies and mean-field treatments of simplified models show the emergence of sub-diffusion on the other-hand in the hydrodynamic low k limit. For the hard-core case, we derive a dynamical equation using a mean-field like assumption which is able to capture partially the unconventional feature of an increase in slope around the nodes of the pattern driven by the non-linear terms in the dynamical equation but which also predicts a time scaling of k −4 . This leads us to the conclusion that the microscopic dynamics are controlled by strong correlations which are built up in the prethermal regime, which cannot be accounted for by this mean-field treatment.
To understand the mechanism driving the long prethermal regime, we study the correlation function of flippable plaquette density for different instances in time. We find that a strong anti-correlation is built up in the direction orthogonal to the modulation direction for the initial condition, even extending to the entire length of the diagonals for intermediate system sizes. We thus observe the formation of isolated dynamical puddles and find that the mechanism leading to the melting is characterized by a merger of these puddles, which relax into a profile which is consistent with what we expect from the continuum theory. Thus, we have found an example where correlations control crucial aspects of the dynamical behavior, and a simple mean-field like hydrodynamic description is insufficient. An important question remains: how generic are the long melting processes that we observe, in particular to which additional classes of initial conditions and types of dynamics do they apply? From the discussion of perfect stripe configurations presented in Sec. II, it is easy to see that similar results can also be expected for stripes which have a different tilt, as long as k x , k y = 0. Regarding the possible set of dynamics which can be expected to engineer a similar phenomenology, we observe that a restricted spread of excitations similar to fractonic systems [29] may be at play here. For the specific model we study here this can be seen by considering the example discussed in Sec. II, where we see that the plaquette excitation is not able to move freely in the direction where there is a lack of similar excitations. More complex versions of this example can be generated by other spatially-larger dynamical terms conserving dipole moment, and can be expected to have a similar constraint when diffusing into spatial regions which lack dynamically active cells. This mechanism is nevertheless not as cleanly formulated and generic as the ones found in fractonic models [29,53,54]. In particular it is still possible to define special particle configurations in our model, in which an active plaquette can move locally without restriction. However, our numerical results and the intuition developed in Sec. II suggest that for modulated structures which vary over large length scales, the phenomenology which we have discussed should hold for generic dipole moment conserving local dynamics.
We have shown that the decay of the overlap with the initial state with time shows a similar behavior upon comparing the exact quantum dynamics and our classical automaton for a small system size. This leads us to speculate that the behavior studied in the automaton language might inherit properties which can also be relevant for the exact quantum evolution for long times, as already seen for the study of sub-diffusive behavior in 1D [24]. First of all, we expect that our hydrodynamic description might be equally valid also for the quantum dynamics, in case the quantum system exhibits a hydrodynamic long-time behavior. This would imply sub-diffusion also in the quantum case. Further, it might be possible that the onset of hydrodynamic behavior also in the quantum model might be delayed to very long times. This would have the consequence that a large intermediate time window exists with unconventionally slow dynamics. As this model is related to the quantum link model [15] via a relaxation of charge conservation [50], similar physics may be expected to emerge there in a more restricted setting. The category of initial conditions which we have studied may be seen as high energy states of a quantum XY model with ring exchange interactions [55]. We expect that an initial state with a fully staggered charge density pattern, located at low-energies for such a model with large ring-exchange, would relax extremely quickly as it has the maximal possible density of flippable plaquettes by the ring-exchange term.
Given that the stripe configurations we consider are qualitatively similar to domain walls in some ferromagnetic systems (such as found e.g. in Ising models [56]), it is important to note that the stability of domain walls in such cases is ensured by energetic reasons and is only expected at low temperature. In contrast, we do not have any such potential term which reinforces local correlations, and our long lived stripe structures are generated purely by the restricted dynamics available to the system. The crossover between Hilbert space fragmentation, which dictates preservation of the memory of the initial state to arbitrarily long times, and the long prethermal plateaus we see here, is another promising direction in which investigations can be carried out to better our understanding of the processes involved. Although we do not find any exact conserved quantities which could explain the slow relaxation seen in our dipole conserving model, another interesting follow-up would be to try to build a theoretical framework in terms of statistically localized integrals of motion [57].
Finally, the rapid increase of the lifetime of the prethermal behavior with inverse wave-vector suggests that rare events, which create configurations allowing the applicability of the hydrodynamic theory, play a key role. An improved understanding of the potentially large deviations which lead to the above mentioned phenomenon would help greatly in developing a coarse-grained description of the dynamics at the threshold between the prethermal and equilibrium regimes. tic dynamics. The authors of Ref. [51] estimate the correlation function in momentum space (C( k, t)) for a time much larger than the scale of the microscopic dynamics (using both automaton dynamics and an effective analytical ansatz) and find sub-diffusive features with "hidden" modulated symmetries corresponding to certain patterns in the Brillouin zone (see Fig.2b of Ref. [51]).
In this appendix, we would like to check if this analysis can help identify the slow dynamics we observe for large scale modulated patterns. The energy spectrum as a function of k has already been identified in a meanfield framework by Paramekanti et al. in Ref. [45] for the form of ring-exchange used in our work, to be given by E k ∝ | sin(k x /2) sin(k y /2)|. This trivially implies that the lines k x = 0 and k y = 0 host zero modes, and that this should be visible in C( k, t) in the long time limit.
Note that this profile does not explain the slow dynamics which we have explored in this work, as that would at the minimum require the presence of slow modes along tilted directions such as k x = k y . To check if this expectation is borne out by ring-exchange dynamics, we perform numerical evolutions of random configurations of a 200 × 200 lattice for t = 2 12 and 2 16 , shown in Fig. 14. We can clearly see that the profile is consistent with nonzero values being present only along the lines k x = 0 and k y = 0 at long times. This aspect of the dynamics is consistent with the mean-field treatment developed in the body of the manuscript. This is to be expected since random configurations correspond to the A → 0 limit of our study of modulated patterns, which was shown to be well-described by the "vanilla" mean-field dynamics.