Spectral statistics in constrained many-body quantum chaotic systems
SSpectral statistics in constrained many-body quantum chaotic systems
Sanjay Moudgalya, Abhinav Prem, David A. Huse, and Amos Chan Department of Physics, Princeton University, Princeton, NJ 08544, USA Princeton Center for Theoretical Science, Princeton University, Princeton, NJ 08544, USA (Dated: September 25, 2020)We study the spectral statistics of spatially-extended many-body quantum systems with on-siteAbelian symmetries or local constraints, focusing primarily on those with conserved dipole andhigher moments. In the limit of large local Hilbert space dimension, we find that the spectralform factor K ( t ) of Floquet random circuits can be mapped exactly to a classical Markov cir-cuit, and, at late times, is related to the partition function of a frustration-free Rokhsar-Kivelson(RK) type Hamiltonian. Through this mapping, we show that the inverse of the spectral gap ofthe RK-Hamiltonian lower bounds the Thouless time t Th of the underlying circuit. For systemswith conserved higher moments, we derive a field theory for the corresponding RK-Hamiltonian byproposing a generalized height field representation for the Hilbert space of the effective spin chain.Using the field theory formulation, we obtain the dispersion of the low-lying excitations of the RK-Hamiltonian in the continuum limit, which allows us to extract t Th . In particular, we analyticallyargue that in a system of length L that conserves the m th multipole moment, t Th scales subdiffu-sively as L m +1) . Our work therefore provides a general approach for studying spectral statistics inconstrained many-body chaotic systems. I. INTRODUCTION
Recent years have seen a surge of interest in under-standing the foundations of quantum statistical mechan-ics. A convergence of experimental progress in the en-gineering and manipulation of ultracold atomic gases,which provide excellent examples of isolated quantumsystems [1], and profound theoretical insight has broughtto the forefront of contemporary research the nature ofclosed quantum many-body systems evolving under uni-tary dynamics. Research in this direction has uneartheda plethora of novel non-equilibrium phenomena, suchas many-body localization [2–5], quantum many-bodyscarring [6–12], and Hilbert space fragmentation [13–16], which provide examples of non-integrable interact-ing systems which fail to obey the Eigenstate Thermal-ization Hypothesis (ETH) [17–21]. Concurrently, ideasfrom the dynamics of black holes have led to new per-spectives on characterizing quantum chaos and diagnos-tics thereof [22], including the decay of Out-of-Time-Order Correlators [23–27] and operator spreading [28–37]. These diagnostics complement familiar signaturesof quantum chaos derived from the eigenvalue spectrumof Hamiltonian or Floquet systems, such as level repul-sion [38–40] and the Spectral Form Factor (SFF) [41].These ideas rely on the widely-held belief that dynamicsof generic many-body quantum chaotic systems beyonda timescale t Th , dubbed the “Thouless time,” follow pre-dictions from Random Matrix Theory (RMT) [42] i.e.,their late-time behavior resembles that of a random ma-trix chosen from an ensemble consistent with the system’ssymmetries [27, 41, 43].Despite the significant difficulty in analytically study-ing dynamics in generic many-body quantum systems,substantial progress has been made in delineating thedynamics of chaos in random quantum circuits via thetwo-point spectral form factor K ( t ), defined in terms of the spectral properties of the evolution operator (cid:99) W as K ( t ) := (cid:42) N (cid:88) m,n =1 e i ( θ m − θ n ) t (cid:43) = (cid:28)(cid:12)(cid:12)(cid:12) Tr[ (cid:99) W ( t )] (cid:12)(cid:12)(cid:12) (cid:29) , (1)where { θ m } is the set of eigenphases of (cid:99) W , (cid:99) W ( t ) ≡ (cid:99) W t denotes the t th power of (cid:99) W , N is the Hilbert space dimen-sion, and (cid:104)·(cid:105) denotes the average over an ensemble of sta-tistically similar systems. The SFF is the Fourier trans-form of the two-level correlation function, with time t asthe variable conjugate to ω ∼ θ m − θ n , the (quasi)-energydifference. For (cid:99) W with Poisson-distributed eigen-levels, K ( t ) = N for all t , while for (cid:99) W chosen as random matri- FIG. 1. (Color online) Illustration of the SFF versus time on alog-log scale. The dashed blue line is the CUE RMT behavior.The black solid line is the generic behavior for many-bodyquantum chaotic systems, characterized by the Thouless time t Th and the Heisenberg time t Hei , the latter of which scalesexponentially with system size L . The early time behaviourof K ( t ) depends on the details of the system, particularlyon whether the system is defined by a Floquet operator or aHamiltonian. a r X i v : . [ c ond - m a t . s t a t - m ec h ] S e p FIG. 2. Summary of results: a) We study a class of constrained Floquet random quantum circuits whose SFF in the large- q limit can be related to the partition function of an RK-Hamiltonian. b) Graph representation of the RK-Hamiltonian, withnodes representing an appropriate set of basis states and links representing possible moves between them. c) Generalized heightfield representation for each spin configuration, with symmetries enforced via boundary constraints. The RK ground state isthen written as an equal-weight superposition of height fields (equivalently, as random walks in the space direction). ces from the Circular Unitary Ensemble (CUE), K ( t ) = | t | (the “ramp”) for times well below the Heisenberg time t Heis = N , after which it plateaus to K ( t ) = N [41]. TheSFF serves as a barometer for quantum chaos, as the ap-pearance of RMT predictions in the SFF provide a crisptime-scale for characterizing the many-body chaos in afinite system, which we take as defining the Thoulesstime t Th1 ; specifically, we define t Th as the time scaleat which the SFF approaches the CUE RMT behavior(see Fig. 1). In contrast to nearest-neighbour level spac-ing distributions, the SFF encodes spectral correlationsat all time (equivalently, energy) scales and is also rela-tively simple to analyse, as it involves only two sums overthe (quasi)-energy eigenvalues.One approach for analytically computing the SFF ex-ploits a self-duality present in certain models [44, 45].However, the applicability of this approach is limited onlyto self-dual circuits and does not extend to generic in-teracting circuits, possibly with conserved quantities. Asecond, more generic approach [46–50] studies Floquetrandom quantum circuits (FRQCs) in the limit of largelocal Hilbert space dimension, which are amenable to ex-act analytic calculations of the SFF and hence enableone to study the implications of conserved quantities ondynamics. For circuits with a globally conserved U(1)charge, t Th was shown to scale diffusively as ∼ L , vali-dating the idea that RMT behaviour is established onlyonce all conserved quantities have diffused through thesystem. In contrast, for certain systems without any con-served charge, t Th ∼ log L [47]. In these latter systems, In systems with diffusive or subdiffusive dynamics, this is ex-pected to be a good estimate of the actual Thouless time, sinceit is always slow compared to the operator spreading time. the SFF is not sensitive to the slower “ballistic” dynamicsdue to locality and causality which means the operatorspreading differs from RMT dynamics up to longer timesof order L .Besides systems with conserved charges, there has beengrowing interest in the dynamics of constrained non-integrable quantum systems with more general symme-tries, driven partly by the discovery of their anomalousdynamics [8, 10, 51–58], which resembles that of classicalKinetically-Constrained-Models [59]. Of particular inter-est are systems which conserve both charge and dipolemoment (or center-of-mass) [13, 15, 16, 60–64], sym-metries which naturally appear in systems subjected tostrong electric fields [14, 16, 65] and in fracton models [66,67]. Dipole moment conserving systems exhibit variousnovel dynamical phenomena including operator localiza-tion [60] and Hilbert space fragmentation [13, 14, 61],coexistence of integrable and non-integrable subspacesleading to a restricted form of ETH [16, 62], presenceof topological edge modes in highly excited states [15],and subdiffusive transport [64, 68, 69], which has beenexperimentally observed [65].In this paper, we develop a general approach for study-ing features of the spectral statistics in constrained quan-tum chaotic many-body systems, focusing on those withconserved dipole and higher moments. Using the SFF K ( t ) as a diagnostic for many-body chaos, we extract thescaling of this t Th with system size L for one-dimensional(1D) FRQCs with a local Hilbert space comprising q ‘color’ degrees of freedom (DOFs) coupled with auxiliaryspins through which the constraints are imposed. In thelarge- q ( q → ∞ ) limit, we express K ( t ) in terms of a clas-sical Markov circuit which inherits the constraints of theunderlying FRQC, as illustrated in Fig. 2(b). Utilizingan established correspondence between classical Markovprocesses and Rokhsar-Kivelson (RK) type Hamiltoni-ans, we can equivalently relate K ( t ) at late times to thepartition function of a positive-definite, frustration-freeRK-Hamiltonian acting solely on the spin DOFs. Con-sequently, a lower bound on t Th can be extracted fromthe spectral gap of the RK-Hamiltonian; this mappinghence allows us to borrow techniques from equilibriumphysics to establish dynamical properties of the underly-ing circuit [71–73].For circuits with conserved higher moments, we finda continuum representation for the ground state (GS)of the RK-Hamiltonian in terms of generalized “height”fields (see Fig. 2(c)), from which we identify a continuumparent Hamiltonian for the corresponding GS and estab-lish a lower bound on t Th . We find a sub-diffusive scalingof t Th ∼ L m +1) for circuits of length L with a con-served m th moment i.e., the timescale at which randommatrix behavior ensues in such systems is parametricallylonger than in systems with only a conserved U(1) charge,which spreads diffusively ( m = 0: t Th ∼ L ). Note thatwhile the height field representation is specific to higher-moment conserving systems, the mapping from K ( t ) toan emergent RK-Hamiltonian in the large- q limit holdsgenerally.This paper is organized as follows. We start by defin-ing the class of FRQCs under consideration in Sec. II.In Sec. III, we study these circuits in the limit of largelocal Hilbert space and establish a mapping between theSFF of the FRQC and the dynamics of a classical Markovchain, which we further show is equivalent at late times tothe partition function of an emergent RK-Hamiltonian.Through these mappings, we find that the Thouless time t Th of the underlying circuit is lower bounded by thespectral gap of this emergent Hamiltonian. Focusing onsystems with conserved higher moments in Sec. IV, weverify that t Th ∼ L in charge conserving systems andprovide numerical evidence for subdiffusive scaling of t Th in systems that additionally conserve the global dipolemoment. In Sec. V, we take the continuum limit of theemergent RK-Hamiltonian for systems which conserve allmoments up to the m th highest moment. We extract abound on t Th from the dispersion relation of the con-tinuum Hamiltonian and show that t Th ∼ L m +1) . Weconclude in Sec. VI with a discussion of open questionsand future directions. II. CONSTRAINED FLOQUET RANDOMQUANTUM CIRCUITS
Our primary object of interest in this paper is theThouless time t Th of constrained many-body quantum As discussed later in Sec. III, we take “RK-Hamiltonian” to meana quantum Hamiltonian that is proportional to the transitionmatrix of a Markov process that satisfies detailed balance. Thistaxonomy stems from the quantum dimer context [70], wherethe terms “RK-Hamiltonian” and “quantum dimer model at theRK-point” are often used interchangeably. chaotic systems, where we define t Th as the time-scaleafter which the behaviour of the SFF K ( t ) closely ap-proaches RMT predictions. To probe K ( t ), we considerone-dimensional L -site spatially-random FRQCs with lo-cal Hilbert space at each site of the chain given by H loc = C q ⊗ C s +1 , where C q and C s +1 are the localHilbert spaces of the color and spin DOFs respectively.The color DOFs facilitate Haar averaging and allow us toretain analytical control in the q → ∞ limit [36, 46, 47];the spins, on the other hand, allow us to encode on-siteAbelian symmetries, such as U(1) charge conservation(previously considered in Refs. [36, 48]), or impose localconstraints on the dynamics.More precisely, we consider Floquet circuits defined bya time-evolution operator (cid:99) W over a single period com-posed of unitary gates acting on a finite number (cid:96) ≥ (cid:96) min of contiguous sites, where (cid:96) min is the “minimal” gate sizefor non-trivial local dynamics under the symmetry orconstraints of interest. Without loss of generality, localdynamics on all sets of (cid:96) contiguous sites within a singletime-period can be ensured by choosing (cid:99) W to be com-posed of (cid:96) layers of operators { (cid:99) W a } , where (cid:99) W a is com-posed of r = (cid:98) L/(cid:96) (cid:99) spatially random local unitary gates { (cid:98) U [ j,j + (cid:96) − } , and has the form: (cid:99) W = (cid:96) (cid:89) a =1 (cid:99) W a , (cid:99) W a = r (cid:79) n =1 (cid:98) U [ a +( n − (cid:96),a + n(cid:96) − , (2)where a is the layer index. As shown in Fig. 2(a), (cid:98) U [ a +( n − (cid:96),a + n(cid:96) − labels the n th local gate in the a th layer and acts non-trivially only on sites j ∈ { a + n(cid:96) − (cid:96), . . . , a + n(cid:96) − } , where 1 ≤ n ≤ r and the site index j isdefined mod L for periodic boundary conditions (PBC).Each of the local gates has the following block-diagonalstructure: (cid:98) U [ j,j + (cid:96) − = D (cid:77) α =1 u ( j, α ) , (3)where α denotes each set (block) of (cid:96) -site spin-configurations within this gate which are connectedthrough local moves permissible under symmetries orconstraints, and D denotes the total number of suchblocks. Block α contains d α spin configurations withinthis gate, with (cid:80) D α =1 d α = (2 s + 1) (cid:96) . Note that we do notimpose any constraints on the color DOFs, only on thespins. Each u ( j, α ) is thus a d α q (cid:96) × d α q (cid:96) unitary drawn in-dependently from the Haar ensemble acting on the statesin block α , while it gives zero when it acts on all other In certain cases, considering Floquet operators with m < (cid:96) layersper period is sufficient to ensure that non-trivial dynamics occursin all sets of (cid:96) contiguous sites. However, this choice of m < (cid:96) layers leads to identical late-time dynamical features as that ofoperators with (cid:96) layers. states. In particular, for systems without any symme-tries or dynamical constraints, the local gates (cid:98) U [ j,j + (cid:96) − are (2 s + 1) (cid:96) q (cid:96) × (2 s + 1) (cid:96) q (cid:96) independent Haar randomunitaries. The block diagonal structure of the local gatesencodes the symmetries or constraints of interest. Specif-ically for systems which have global symmetry sectorslabelled by a set of quantum numbers S = { s , s , . . . } ,each local (cid:96) -site gate (cid:98) U is block diagonal, with each blockcontaining all spin states with the same S . As a technicalaside, we note that since we are keeping the gate size (cid:96) fixed while treating all transitions involving those (cid:96) siteson equal footing, this also includes all allowed processesinvolving spin transitions on any subset of those (cid:96) sites.As an example, let us consider an FRQC with s = 1 / (cid:98) Q = (cid:80) Lx =1 (cid:98) S zx of the spins, where (cid:98) S zx is the Pauli- Z matrix acting onsite x [48]. To allow non-trivial dynamics, we choose (cid:96) tobe equal to (cid:96) min = 2, such that each local gate (cid:98) U is a4 q × q block-diagonal matrix. Each local gate is com-posed of D = 3 blocks: two q × q blocks act on thetensor product subspaces associated with the spin con-figurations |↑↑(cid:105) and |↓↓(cid:105) , and a single 2 q × q block actson the subspace associated with the spin configurations |↑↓(cid:105) , |↓↑(cid:105) . Each of these three blocks locally preserves theU(1) charge over 2-sites and is an independently-drawnHaar random unitary.In this paper, we mainly focus on circuits which con-serve not only the total charge, but also all higher mo-ments up to the m th moment (cid:98) Q m = (cid:80) Lx =1 x m (cid:98) S zx ≡ (cid:80) Lx =1 (cid:98) Q m ( x ). Such circuits neatly fall into the largerclass of FRQCs defined earlier via Eqs. (2) and (3). Forinstance, for systems with both charge and dipole mo-ment conservation, we can consider s = 1 DOFs and (cid:96) = (cid:96) min = 3 site gates, where each local gate is a 27 q × q block diagonal matrix, with each block corresponding tothose spin configurations which are connected under local(3-site) dipole moment preserving dynamics (see Ref. [60]for details). It is straightforward to generalize the abovecircuits to higher spins s , larger gate sizes (cid:96) , and highermoment conservation laws. While our primary focus inthis paper will be systems with higher conserved mo-ments, the FRQCs defined in Eqs. (2) and (3) define amuch broader class, including those with arbitrary on-site Abelian symmetries as well as circuits which obeydynamical constraints, such as those present in the PXPmodel [56].We characterise the spectral features of the above classof FRQCs using the SFF defined in Eq. (1). For a cir-cuit (cid:99) W invariant under a set of global symmetries, cor-responding to a set of operators { (cid:98) S , (cid:98) S , . . . } , we have[ (cid:99) W , (cid:98) S i ] = 0 ∀ i . Therefore, (cid:99) W = (cid:77) S (cid:99) W ( S ) (4)is block-diagonal and quasienergy levels of (cid:99) W from blocks (cid:99) W ( S ) , corresponding to distinct quantum number sec- FIG. 3. (Color online) Example of a Floquet operator ex-hibiting Hilbert space fragmentation in the Z -basis. Symme-try sectors are denoted by S , and Krylov subspaces K withinsymmetry sectors are denoted by K ( S ) i . tors S = { s , s , . . . } , do not repel [21]. In addition, cer-tain systems, such as those with higher moment symme-tries, further exhibit the phenomenon of Hilbert spacefragmentation, wherein the dynamics does not connectall products states in the (cid:98) S z (equivalently, charge) basiseven within the same symmetry sector [13, 14, 61, 62].As a consequence, within each symmetry sector S theremay exist up to exponentially many disjoint “Krylov sub-spaces,” labelled by K ( S ) i i.e., (cid:99) W ( S ) = D ( S ) (cid:77) i =1 (cid:99) W ( K ( S ) i ) , (5)where D ( S ) denotes the number of disjoint Krylov sub-spaces generated from product states with the samequantum numbers S . This fragmented structure of (cid:99) W isschematically depicted in Fig. 3. Hence, given the possi-bility of Hilbert space fragmentation in quantum many-body systems, we define the SFF restricted to a givenKrylov subspace K : K ( t ; K ) ≡ (cid:28)(cid:12)(cid:12)(cid:12) Tr K [ (cid:99) W t ] (cid:12)(cid:12)(cid:12) (cid:29) , (6)where the subscript K denotes the restriction of (cid:99) W to aKrylov subspace and (cid:104)·(cid:105) denotes averaging over the Haarrandom unitaries in the FRQC. Note that this definitionencompasses systems with global symmetries but no frag-mentation, since in that case, each Krylov subspace K ( S ) fully spans its global symmetry sector S and all D ( S ) = 1. III. MAPPING TO CLASSICAL MARKOVCHAIN AND EMERGENT RK-HAMILTONIAN
Computing the SFF Eq. (6) for many-body quantumsystems is analytically difficult in general. Refs. [46–48]developed a diagrammatic approach for evaluating theensemble averaging in (6), by effectively “integrating out”the color DOFs for random quantum circuits with chargeconservation. In Appendix A, we generalize this tech-nique to the general class of constrained FRQCs dis-cussed in the previous section and find that, to leadingorder in the large- q limit, K ∞ ( t ; K ) ≡ lim q →∞ K ( t ; K ) = | t | Tr K (cid:104) (cid:99) M t (cid:105) , (7)where the factor of | t | stems from | t | leading order dia-grams as q → ∞ [46]. The Markov matrix (cid:99) M is a clas-sical bi-stochastic circuit acting only on effective spin- s DOFs and is composed of local (cid:96) -site gates. Furthermore, (cid:99) M inherits the circuit geometry, symmetries, and Krylovsubspaces of (cid:99) W and can be expressed as (cid:99) M = (cid:96) (cid:89) a =1 (cid:99) M a , (cid:99) M a = r (cid:79) n =1 (cid:98) m [ a +( n − (cid:96),a + n(cid:96) − (8)where, in analogy with Eq. (2), a is the layer index and (cid:98) m [ a +( n − (cid:96),a + n(cid:96) − labels the n th local gate in the a th layer, acting on sites j ∈ { a + n(cid:96) − (cid:96), . . . , a + n(cid:96) − } .As before, 1 ≤ n ≤ r and j is defined mod L for PBC.However, unlike the underlying gates (cid:98) U , the (cid:96) -site gates { (cid:98) m [ j,j + (cid:96) − } are non-random : (cid:98) m [ j,j + (cid:96) − = D (cid:77) α =1 m ( d α ) , m ( d α ) = 1 d α . . . . . . ... ... . . . d α × d α , (9)where the d α ’s are the sizes of the blocks of (cid:96) -site spinconfigurations that are dynamically connected and D isthe number of dynamically connected blocks in (cid:98) m . Thefact that Eq. (9) retains the block-diagonal form withequal matrix elements (within each block) is consistentwith the fact that the dynamical constraints are imposedvia the spin DOFs and that the local Haar random gatesare invariant upon a change of basis.Since (cid:99) M inherits the symmetries and Krylov subspaces(Eqs. (4) and (5)) of (cid:99) W , and (cid:99) M also has a block diag-onal structure (see Fig. 3), where each block is itself anirreducible bi-stochastic matrix; hence, each block has aunique largest magnitude eigenvalue 1. Restricting ourattention to the subspace K of interest, let us denote theeigenvalues of the corresponding block by { Λ ( K ) j } , withΛ ( K )1 = 1 and ordered such that | Λ ( K ) j | ≥ | Λ ( K ) j +1 | ∀ j . Wecan then write K ∞ ( t ; K ) = | t | Tr K (cid:104) (cid:99) M t (cid:105) = | t | (cid:88) j> (cid:16) Λ ( K ) j (cid:17) t . (10) An N × N non-negative matrix M is called bi-stochastic if (cid:80) Ni =1 M i,j = 1 ∀ j and (cid:80) Nj =1 M i,j = 1 ∀ i . Note that the eigenvalues of (cid:99) M can be complex since is not asymmetric matrix due to the brick-wall structure of the circuit At sufficiently long times t (cid:29)
1, we can then expandEq. (10) as K ∞ ( t ; K ) = | t | (cid:16) d ( K ) exp (cid:16) − t ∆ ( K ) (cid:99) M (cid:17) + · · · (cid:17) , (11)where | Λ ( K )2 | = exp (cid:16) − ∆ ( K ) (cid:99) M (cid:17) is the magnitude of thesecond largest eigenvalue of (cid:99) M restricted to the subspace K and d ( K ) is a constant factor encoding its degeneracyalong with any complex phases.Using Eq. (11), we see that the SFF K ∞ ( t ; K ) ap-proaches the linear in | t | RMT behavior after a time t (cid:39) ( K ) (cid:99) M ( L ) ≡ t ( K )Th . (12)Thus, we have shown that extracting the scaling behav-ior of t Th with system size L is equivalent, in the large- q limit, to the problem of obtaining the scaling of the“gap” ∆ ( K ) (cid:99) M ( L ) of the Markov circuit (cid:99) M (within the sub-space K ) with L . Henceforth, we will restrict our atten-tion to a single (exponentially large) subspace K of (cid:99) M and suppress the sub/superscript K for ease of notation.To determine the gap ∆ (cid:99) M , we will now establish a re-lation between (cid:99) M (within a subspace K ) and a quantumHamiltonian. We proceed by first observing that (cid:99) M cor-responds to the classical stochastic time evolution of aprobability density (cid:126)p ( t ), defined over all product statesin the usual Z -basis for the spin Hilbert space: p α ( t + 1) = (cid:88) β (cid:99) M αβ p β ( t ) , (13)where (cid:99) M αβ represents the matrix elements of (cid:99) M and thebistochasticity of (cid:99) M ensures that the total probabilityis conserved under time-evolution. In particular, underthe action of each local gate (cid:98) m [ j,j + (cid:96) − (see Eq. (9)), theprobability density (cid:126)p evolves (with equal probability) toall product states that can be reached via the allowed lo-cal moves i.e., moves that are allowed by the constraintsimposed on the spin DOFs. Now consider starting with aprobability density where all the weight is concentratedon a single product state within some subspace K . Un-der the stochastic evolution, the probability density willeventually reach a unique equilibrium state, specified bythe uniform distribution over all product states in K i.e.,by the eigenvector of (cid:99) M corresponding to the eigenvalueΛ = 1. Thus, obtaining ∆ (cid:99) M is related to obtaining the Technically, we require that t H (cid:29) t (cid:29)
1, but the so-calledHeisenberg time t H proportional to the inverse level spacing isinfinite in the large- q limit. The gap of a bistochastic matrix is traditionally defined to be1 − exp (cid:16) − ∆ ( K ) (cid:99) M (cid:17) , which we have approximated as ∆ ( K ) (cid:99) M here. inverse of the mixing time for this process, which is ingeneral not analytically tractable.To derive the gap ∆ (cid:99) M , we are interested in the stochas-tic evolution under (cid:99) M at time-scales of O (1 / ∆ (cid:99) M ). If∆ (cid:99) M → (cid:99) M is gap-less), we expect that the dynamics under (cid:99) M at latetimes is well approximated by a continuous-time pro-cess. In particular, since (cid:99) M corresponds to a stochasticprocess with local moves occurring independently withequal probability, its late-time behavior should be wellapproximated by that of a continuous time process com-posed of the same local moves occurring at equal rates,with the additional requirement of detailed balance. Inother words, the evolution of the probability density (cid:126)p should be governed by a Master equation of the form: dp α ( t ) dt = (cid:88) β (cid:54) = α ( T αβ p β ( t ) − T βα p α ( t )) , (14)where p α ( t ) is the probability of a classical system oc-cupying state α and T αβ is the transition rate fromstate β to state α , which we specify below. Defining T αα ≡ − (cid:80) β (cid:54) = α T βα , we can rewrite Eq. (14) as a matrixequation in terms of the transition matrix T , d(cid:126)p ( t ) dt = T (cid:126)p ( t ) , (15)which ensures that local moves that occur with equalprobability in the discrete-time process Eq. (13) occur atequal rates in the continuous-time process Eq. (16). Thelate time behavior of the stochastic process governed by (cid:99) M is then given by Eq. (15) with T = − Γ H = ⇒ ˙ (cid:126)p ( t ) = − Γ H(cid:126)p ( t ) , (16)where Γ is an overall positive constant that sets the rateat which local moves occur and H is defined as H = (cid:88) j Π [ j,j + (cid:96) − , Π [ j,j + (cid:96) − ≡ − (cid:98) m [ j,j + (cid:96) − . (17)The matrix Π [ j,j + (cid:96) − is a projector and has the formΠ [ j,j + (cid:96) − = D (cid:77) α =1 π ( d α ) ,π ( d α ) = 1 d α ( d α − − . . . − d α − . . . ... ... . . . d α × d α , (18) Note that this approximation does not hold for gapped systemsas they relax to their equilibrium distribution on time-scales of O (1), that are much smaller than the time-scale at which thecontinuous time description is valid. As we show in Sec. IV, Γ is a non-universal constant determinedby the detailed microscopic properties of the underlying FRQCe.g., it depends on the number of layers (cid:96) . However, obtaining itsprecise value is not important for our purposes. thereby ensuring that the transitions taking place in thecontinuous-time process are identical to those specifiedby the gates (cid:98) m [ j,j + (cid:96) − .In fact, H can be interpreted as a quantum Hamil-tonian in the product state basis (in the Z -basis) forthe spin Hilbert space and has the same symmetriesand Krylov subspaces as the stochastic circuit (cid:99) M . Moreimportantly, H belongs to the class of so-called RK-Hamiltonians , where an RK-Hamiltonian is defined asa quantum Hamiltonian that is proportional to the tran-sition matrix T of a discrete classical stochastic processwhich satisfies detailed balance [72]. Consequently, theground state wave function of an RK-Hamiltonian canbe interpreted as a classical equilibrium distribution, itslow-lying excited states correspond to classical relaxationmodes, and its gap coincides with the relaxation time ofthe corresponding transition matrix. Such Hamiltonianswere first studied in the context of quantum dimer mod-els [70] and have subsequently been explored extensivelyin various settings [71, 72, 74].To emphasize the relation between (cid:99) M and an emer-gent RK-Hamiltonian, we henceforth adopt the notation H → H RK . In the picture developed above, Eq. (16)then has the clear interpretation of an imaginary-timeSchr¨odinger evolution under H RK Eq. (17). In effect, thecorrespondence between (cid:99) M and H RK amounts to a rela-tion between Tr K [ (cid:99) M t ] and the partition function of H RK restricted to K at an inverse temperature β = Γ t , namely:Tr K [ (cid:99) M t ] ≈ Tr K [ e − Γ H RK t ]. We can therefore approximatethe SFF Eq. (10) at late times as K ∞ ( t ; K ) t (cid:29) ≈ | t | Tr K (cid:2) e − Γ H RK t (cid:3) , (19)such that the gap of (cid:99) M is related to ∆ RK ( L ), the gap ofthe Hamiltonian H RK (restricted to the subspace K ):∆ (cid:99) M ( L ) ≈ Γ∆ RK ( L ) . (20)Indeed, as we discuss in Sec. IV, we find numerical evi-dence that supports Eq. (20) in multipole conserving cir-cuits. Further evidence for the correspondence between H RK and (cid:99) M is obtained by studying the system-size de-pendence of the overlap between the “first-excited eigen-states” of | ψ RK (cid:105) and (cid:12)(cid:12) ψ (cid:99) M (cid:11) of H RK and (cid:99) M respectively.As shown in the inset of Fig. 4(a), we find that this over-lap approaches 1, suggesting that | ψ RK (cid:105) is an asymptot-ically exact eigenstate of (cid:99) M in the thermodynamic limit.We also numerically observe that the overlap does notapproach 1 in the cases when (cid:99) M is gapped, further sug-gesting that correspondence between (cid:99) M and H RK is onlyvalid when (cid:99) M is gapless.The preceding discussion shows that obtaining the gapof H RK is sufficient for obtaining the scaling of the Thou-less time t Th . Taken together, Eqs. (12) and (20) consti-tute one of the central results of this paper, whereby adynamical property of the FRQC is determined by thelow-energy, equilibrium behavior of an emergent quan-tum Hamiltonian. Since we made no reference to the mi-croscopic structure of the underlying circuit, this rela-tionship between t Th and ∆ RK holds generally for theclass of circuits specified in Sec. II.We now briefly discuss some properties of the Hamil-tonians H RK . Denoting the basis set of π ( d α ) in Eq. (18)by B α , we obtain π ( d α ) = 1 d α (cid:88) C , C (cid:48) ∈B α ( |C(cid:105) − |C (cid:48) (cid:105) ) ( (cid:104)C| − (cid:104)C (cid:48) | ) , (21)where C and C (cid:48) represent (cid:96) -site configurations chosenfrom the basis set B α . We can then re-express H RK as H RK = (cid:88) (cid:104)C , C (cid:48) (cid:105) w C , C (cid:48) (cid:98) Q C , C (cid:48) , (cid:98) Q C , C (cid:48) ≡ ( |C(cid:105) − |C (cid:48) (cid:105) ) ( (cid:104)C| − (cid:104)C (cid:48) | ) , (22)where (cid:104)C , C (cid:48) (cid:105) denotes product states C and C (cid:48) that areconnected under the action of H RK and the weights w C , C (cid:48) ≥ H RK is a positive semidef-inite Hamiltonian and the zero-energy ground state wavefunction within a subspace K is given by (cid:12)(cid:12)(cid:12) G ( K )RK (cid:69) = 1 √Z (cid:88) C∈K |C(cid:105) , (23)where C runs over all the product states in the subspace K and Z is a normalization factor. We remark that, up to anoverall normalization factor, | G RK (cid:105) can be interpreted asthe equilibrium probability distribution of the stochasticprocess described by Eqs. (13) or (16).Before proceeding to focus on multipole conservingcircuits, we briefly discuss some important aspects ofthe RK-Hamiltonian H RK which illustrate the potentialbenefits of the mapping developed in this section. First, H RK is a frustration-free Hamiltonian, i.e., the groundstate | G RK (cid:105) is the ground state of each of the termsΠ [ j,j + (cid:96) − as can be seen using Eqs. (22) and (23), since (cid:98) Q C , C (cid:48) (cid:12)(cid:12)(cid:12) G ( K ) RK (cid:69) = 0 for any C and C (cid:48) . This fact enables theuse of well-known methods for bounding the spectral gapof frustration-free Hamiltonians [75–77], which in turn al-low us to place a constraint on the scaling exponent of t Th in FRQC with constraints defined in Eq. (2) in thelarge q limit: t Th ∼ L α , (cid:40) α = 0 or α ≥ , PBC α = 0 or α ≥ / , OBC (24)We note that for circuits without any conserved quanti-ties, t Th has been shown to scale as log L for certain Flo-quet models [32, 47]. However, as evidenced by Eq. (24),for circuits of the form Eq. (2) (including the circuit dis-cussed in Ref. [46]), such scaling of t Th is suppressed inthe q → ∞ limit and we instead find that t Th can scaleas an O (1) number in FRQCs in this limit. Secondly, by virtue of the connection to classical Mas-ter equations, shown in Eq. (16), H RK is an example ofa stoquastic Hamiltonian, which can be efficiently stud-ied using Quantum Monte Carlo techniques [72, 78]. Wethus expect that the same techniques can be exploitedto efficiently study the late-time features of the SFF ina variety of settings at large- q . Moreover, in the contextof spectral graph theory, any Hamiltonian of the formEq. (22) restricted to a subspace K exactly correspondsto the Laplacian [79] of an undirected graph G , formedby the set of vertices {C} within K and by edges withweights w C , C (cid:48) between the vertices C and C (cid:48) . The gap∆ RK then corresponds to the gap of the Laplacian of thegraph G , which is closely related to the connectivity of G . In particular, the existence of bottlenecks in G , as de-tected by the Cheeger constant, results in a smaller gapof the Laplacian (and hence, in a larger t Th ). This estab-lishes a clear connection between the nature of transportin the presence of constraints and the connectivity of theHilbert space under those constraints.Finally, we note that earlier work has also discusseda relation between the Thouless time t Th of a charge-conserving FRQC and and the spectral gap of a U(1) in-variant classical bistochastic circuit: Ref. [48] constitutesa particular case of the results obtained in this paper,being derived in the large- q limit, while Ref. [80] invokesthe random phase approximation in a long-range inter-acting model at finite- q . While both of these works wererestricted to specific realizations of U(1) invariant sys-tems, the relations between K ∞ , (cid:99) M , and H RK obtainedin this section apply far more generally to the large classof circuits with arbitrary symmetries or constraints dis-cussed in Sec. II. IV. EXAMPLES FROM MULTIPOLECONSERVING CIRCUITS
In this section, we move our attention to FRQCs withconserved higher moments and provide explicit examplesof the mapping established in Sec. III. Specifically, weconsider circuits which conserve all moments of chargeup to the m th highest moment, where the m th multipolemoment is defined as (cid:98) Q m = L (cid:80) x =1 x m S zx , OBCexp (cid:18) πi L (cid:80) x =1 ( x/L ) m S zx (cid:19) , PBC (25)where m = 0 (1) corresponds to the charge (dipole) con-serving case. Where necessary, we will use the labels { Q m } to denote the m th multipole moment quantumnumber.We start by reviewing the charge conserving FRQC(see Sec. II), which was previously discussed in Ref. [48].Following the general discussion in Secs. II and III, thestochastic circuit (cid:99) M for a spin-1/2 U(1) charge conserv- . . . . . . . . ( L ) − . − . − . − . − . − . . . l og ( ∆ ) (a) H RK c M L . . . | h ψ R K | ψ c M i | m = 0 , s = , ‘ = 3 m = 1 , s = 1 , ‘ = 4 0 . . . . . . . . . ( L ) − . − . − . − . − . − . . . l og ( ∆ R K ) (b) m = 0 , s = , ‘ = 2 m = 0 , s = , ‘ = 3 m = 1 , s = 1 , ‘ = 4 m = 1 , s = 1 , ‘ = 5 L − . L − . L − . L − . FIG. 4. (Color online) Scaling of the gaps ∆ (cid:99) M ( L ) and ∆ RK ( L ) of (cid:99) M and H RK respectively for a system of spin- s and (cid:96) -sizedgates with conserved charges { Q j } , j ≤ m . (a) Gaps of H RK and (cid:99) M for two systems: one with charge conservation and one withdipole conservation. Note that for large system sizes, the gaps are related by a constant factor (see Eq. (20)). The inset showsthe overlap of normalized “first-excited eigenstates” of H RK ( | ψ RK (cid:105) ) and (cid:99) M ( (cid:12)(cid:12) ψ (cid:99) M (cid:11) ). Note that the overlap approaches 1 withincreasing system size, suggesting that the | ψ RK (cid:105) is an asymptotically exact first-excited state of (cid:99) M . (b) Scaling of the gapsof H RK with system size for systems with charge (green, blue) and dipole moment (orange, red) conservation. Note that thegaps scale diffusively ( ∼ L − β , β ≈
2) in the presence of only charge conservation whereas they scale subdiffusively ( ∼ L − β with β >
2) in the presence of dipole moment conservation. All data presented corresponds to the largest symmetry sector/Krylovsubspace containing the state | · · · (cid:105) for OBC. ing circuit with gate size (cid:96) = 2 (for PBC) is given by (cid:99) M = (cid:79) j odd (cid:98) m [ j,j +1] (cid:79) j even (cid:98) m [ j,j +1] , (cid:98) m [ j,j +1] = / / / / , (26)where the local 2-site gate (cid:98) m [ j,j +1] is written in the(ordered) basis {|↑↑(cid:105) , |↑↓(cid:105) , |↓↑(cid:105) , |↓↓(cid:105)} . Using Eqs. (18)and (26), we find that Π [ j,j +1] maps onto a ferromagneticspin-1/2 Heisenberg termΠ [ j,j +1] = / − / − / / = 12 ( |↑↓(cid:105) − |↓↑(cid:105) ) ( (cid:104)↑↓| − (cid:104)↓↑| ) j,j +1 = 14 (1 − (cid:126)σ j · (cid:126)σ j +1 ) , (27)where Π [ j,j +1] is written in the (ordered) basis {|↑↑(cid:105) , |↑↓(cid:105) , |↓↑(cid:105) , |↓↓(cid:105)} and (cid:126)σ = ( σ x , σ y , σ z ), with σ i theusual Pauli matrices. For the charge conserving FRQC(with (cid:96) = 2), we hence find that H RK is the Bethe-Ansatzintegrable ferromagnetic Heisenberg model, whose inte-grability was exploited in Ref. [48] to study the late-timebehavior of the SFF for this circuit. According to Eq. (23), the unique ground state withinany charge sector Q is the equal amplitude superpo-sition of all product states within that symmetry sec-tor. Indeed, such a state belongs to the SU(2) multi-plet of the spin-polarized ferromagnetic state with to-tal spin Q = L . Moreover, the low-energy excitationsabove the ferromagnetic state in the Heisenberg model ofEq. (27) are exactly known to be spin waves with disper-sion (cid:15) ( k ) = 2 sin ( k/ Q = L − H RK in any Q (cid:54) = L sector is then given by∆ RK ( L ) = (cid:15) (cid:18) k = 2 πL (cid:19) ≈ π L , (28)which is the energy corresponding to the lowest non-zeromomentum spin-wave. Using Eqs. (12) and (20), we findthat the Thouless time in any quantum number sectorin an FRQC with U(1) charge conservation scales diffu-sively with system size i.e., t Th ∼ L . For charge con-serving systems with higher spins or larger gate sizes (cid:96) , H RK is no longer integrable in general, but, as shown inFig. 4(b), we numerically observe the same diffusive scal-ing ∆ (cid:99) M ( L ) ∼ ∆ RK ∼ L − for the systems we studied.In fact, we find that the gaps are identical for spin-1/2and spin-1 systems with gate size (cid:96) = 2 even though therest of the spectrum is different, strongly suggesting auniversal origin of the scaling.As evidenced through the above example, the corre-spondence between the FRQC and the RK-Hamiltonianunveils a curious feature of the large- q limit. While theoriginal FRQC only has U(1) symmetry, after Haar av-eraging and taking q → ∞ , K ∞ ( t ; K )—related to H RK through Eq. (19)—exhibits an enlarged SU(2) invari-ance in the spin DOFs. Indeed, we expect this enlargedsymmetry to be a generic feature in the large- q limit,since RK-points typically exhibit enhanced symmetries,although not necessarily SU(2) [81–83]. On the otherhand, to our knowledge, the emergent integrability inthe above example is not generic and is specific to thespin-1/2 system with 2-site gates.We now turn our attention to systems which conservethe dipole moment (cid:98) Q in addition to the charge (cid:98) Q ,for which the nature of low-energy excitations above theground state of the corresponding RK-Hamiltonian H RK is not immediately apparent. An additional feature insuch systems is the fragmentation of the Hilbert space ofthe FRQC [13, 14], which leads to the formation of ex-ponentially many Krylov subspaces (see Eq. (5)). Hilbertspace fragmentation is typically classified into two types:strong or weak, where the size of the largest Krylov sub-space is respectively a zero or non-zero fraction of thetotal Hilbert space dimension within a given quantumnumber sector in the thermodynamic limit. Refs. [13, 14]numerically observed that spin-1 and spin-1/2 dipole con-serving systems with the minimal gate sizes (cid:96) = (cid:96) min = 3and (cid:96) = (cid:96) min = 4 respectively show strong fragmenta-tion whereas the inclusion of moves requiring larger gatesizes leads to weak fragmentation. Furthermore, Ref. [63]found that the nature of fragmentation can vary evenfor a given gate size depending on the quantum num-ber sector. Generically, however, experimentally relevantmultipole conserving systems are expected to show weakfragmentation.In strongly fragmented systems, the ratio between thedimension of the largest Krylov subspace within a sym-metry sector and the size of that symmetry sector ex-ponentially decays to zero in the thermodynamic limit.As a consequence, typical initial states do not thermal-ize [13, 14], although certain initial states do thermal-ize with respect to smaller Krylov subspaces [16, 62]. Incontrast, for weakly fragmented systems there always ex-ists a dominant Krylov subspace K within a given quan-tum number sector, such that its size asymptotically ap-proaches that of the symmetry sector in the thermo-dynamic limit. Due to this, typical eigenstates withina quantum number sector carry non-zero weight in thedominant Krylov subspace of that symmetry sector andlook thermal. As a consequence, frozen configurations,despite being exponential in number, are expected tohave a negligible effect on t Th in a weakly fragmentedsystem. Since our interest in this work is the behavior of generic multipole conserving systems, we focus only onthermalizing weakly fragmented systems here. Hence, wewill study the SFF, the scaling of the Thouless time t Th and the gaps ∆ (cid:99) M ( L ) and ∆ RK ( L ) all restricted to thedominant Krylov subspace K within a specified quantumnumber sector.In Fig. 4, we show the scaling of the gaps ∆ (cid:99) M ( L ) and∆ RK ( L ) for spin-1/2 and spin-1 dipole conserving sys-tems for several gate sizes (cid:96) > (cid:96) min . Fig. 4(a) shows thatthe numerics are in good agreement with Eq. (20) i.e.,they support the correspondence between the stochasticcircuit (cid:99) M and the emergent RK-Hamiltonian H RK de-veloped in Sec. III. In principle, we can also extract themicroscopic constant Γ for a specific circuit by comparingthe gaps for the corresponding (cid:99) M and H RK . Furthermore,as evident from the numerics shown in Fig. 4(b), we findthat ∆ RK ( L ) scales as ∼ L − β with β >
2; thus, the Thou-less time scales subdiffusively t Th ∼ L β ( β >
2) for systemsizes accessible to exact diagonalization. Importantly,this subdiffusive scaling appears to be a generic feature ofweakly fragmented, dipole conserving systems and doesnot show a strong dependence on the microscopic de-tails ( (cid:96) or s ) of the circuit. This mirrors the behavior ofsystems with only charge conservation (Fig. 4(a)), whichshow diffusive scaling of t Th independent of microscopicdetails.Due to the longer gate sizes required for systemsconserving even higher multipole moments Q m ≥ (e.g., (cid:96) min = 8 for spin-1/2 quadrupole conserving systems,with longer gates likely needed for weak fragmentation),we are unable to eliminate finite-size effects in such cases.Nevertheless, the numerics suggest a universality in thescaling of the gaps of charge and dipole conserving RK-Hamiltonians H RK i.e., ∆ RK ( L ) ∼ L − for charge con-serving systems and ∆ RK ( L ) ∼ L − for dipole conserv-ing systems, regardless of the ultraviolet details of therespective Hamiltonians. The appearance of this univer-sality suggests the existence of a universal field theoreticdescription which effectively captures the low-energy be-havior of generic multipole conserving RK-Hamiltonians,such as the scaling of their gap. The derivation of theseuniversal effective field theories will be the subject of thenext section. V. CONTINUUM LIMIT FOR MULTIPOLECONSERVING SYSTEMS
As discussed in the previous section, the ground stateof an RK-Hamiltonian is well-known as being the equal-weight superposition of all states in the correspond-ing Hilbert space. However, understanding the scalingof the gap ∆ RK requires knowledge of low-lying statesabove the GS, which are generically not known ex-actly. Nevertheless, motivated by our numerical obser-vation of a universal scaling of ∆ RK with L for genericcharge and (weakly fragmented) dipole conserving RK-Hamiltonians, we derive continuum field theories for mul-tipole conserving systems through a coarse-graining pro-cedure, detailed in Appendix D. We find that the re-sultant continuum field theories accurately capture the0 n − − s n , φ ( m ) n s n φ (0) n φ (1) n FIG. 5. (Color online) The generalized height variables of anexample spin state with Q = 0 and Q = 0 for L = 30. Notethat the end points φ (2)0 and φ (2) N +1 are fixed to be zero. ground state and low-energy excitations of the corre-sponding RK-Hamiltonians, therefore providing an an-alytic route to understanding the scaling of t Th in theunderlying FRQC.Throughout this section, we will only consider OBC.We denote the number of spins as N and the systemsize as L = N ∆ x , where ∆ x is the lattice spacing. Thecontinuum limit then corresponds to taking the limits N → ∞ and ∆ x → L fixed. For systems which conserve all moments of chargeup to the m th moment (or, m th moment conserving sys-tems), we focus our attention on the quantum numbersector S = { Q = 0 , Q = 0 , · · · , Q m = 0 } . As dis-cussed in Sec. IV, for weakly fragmented systems thereexists a dominant Krylov subspace within each symmetrysector, such that the size of that subspace asymptoticallyapproaches the size of the full symmetry sector as thegate-size (cid:96) increases. Since taking the continuum limitinvolves coarse-graining and thus effectively taking thegate-size (cid:96) → ∞ , we can neglect the effect of fragmenta-tion in systems with dipole and higher moment conser-vation, and expect that our analysis holds as long as thesectors we study do not exhibit strong fragmentation. A. Generalized height fields
In order to take the continuum limit, we first needto introduce “generalized” height fields, in analogy withfamiliar height fields in the quantum dimer context [83].When taking the continuum limit of the ground statewave function and H RK , a crucial issue is the restrictionto the symmetry sector S ; this restriction imposes a global constraint on the spin DOFs { s n } , thereby resulting in anon-local action for system. It is to circumvent preciselythis issue, while preserving locality, that we work in termsof generalized height variables { φ ( m ) } for systems with m th moment conservation. In terms of these variables,the conservation of higher moments { Q m } are expressedas local boundary constraints on the height fields andtheir derivatives, as opposed to a global constraint onthe spin DOFs.We first illustrate this construction in terms of heightvariables for systems with charge conservation. Theheight DOFs { φ (0) n + } are defined on the links of the one-dimensional chain as φ (0) n + − φ (0) n − = s n or φ (0) n + = φ (0) + n (cid:88) j =1 s n , (29)which immediately suggests that the total charge Q (see Eq. (25)) is given by the flux of the height variablethrough the system: Q = φ (0) N + − φ (0) . (30)Since Eqs. (29) and (30) are invariant under an over-all constant shift of the height variables ( φ n +1 / → φ n +1 / + c ), we can choose φ (0) = 0 without loss ofgenerality. Thus, restricting to a given charge sectorcorresponds to imposing constraints on the boundaryheight variables once we forgo the spin language forthe height representation. Furthermore, as a consequenceof Eq. (29), any charge conserving process involving (cid:96) spins { s n , · · · , s n + (cid:96) − } involves ( (cid:96) −
1) height variables { φ (0) n + , · · · φ (0) n + (cid:96) − } . So, the mapping from the spinDOFs to the height variables preserves the locality ofthe Hamiltonian.We now generalize the height representation to sys-tems with m th moment conservation. For a system withall moments up to the m th moment conserved, we recur-sively define the m th “generalized” height variable φ ( m ) ,that lives on links (sites) if m is even (odd), as φ ( m ) n + − φ ( m ) n − = φ ( m − n if m ∈ Z φ ( m ) n +1 − φ ( m ) n = φ ( m − n + if m ∈ Z + 1 , (31)with { φ (0) n + } defined in Eq. (29). An example of a chargeconfiguration in a dipole conserving ( m = 1) system, ex-pressed in terms of height variables, is depicted in Fig. 5.As mentioned earlier, the key advantage of forgoing thespin language is that the total m th multipole momentscan be expressed in terms of boundary height variablesand their “derivatives”, an observation that will proveessential when taking the continuum limit. For instance,the total dipole moment can be expressed as Q = N (cid:88) j =1 js j = N Q − N − (cid:88) j =1 j (cid:88) k =1 s k = N (cid:16) φ (1) N +1 − φ (1) N (cid:17) − (cid:16) φ (1) N − φ (1)0 (cid:17) , (32)1where we have used Eqs. (30) and (31). Similarly to thecharge conserving case, locality is also preserved whengoing to the generalized height representation. This canbe seen from Eq. (31), since any m th multipole conservingprocess involving (cid:96) spins { s n , · · · , s n + (cid:96) − } involves ( (cid:96) − m + 1) generalized height variables { φ ( m ) n + } .We now take the continuum limit by defining general-ized height fields through an appropriate rescaling of theheight variables by the lattice spacing ∆ x i.e., ρ ( x ) = s n ∆ x , φ ( m ) ( x ) = φ ( m ) n (∆ x ) m , (33)where x = n ∆ x . Here, ρ ( x ) can be interpreted as thecharge density and, as we will see, the height fields φ ( m ) ( x ) are related to the multipolar densities. In termsof the height fields, Eq. (29) in the continuum becomes ∂ x φ (0) ( x ) = ρ ( x ) . (34)Note that Eq. (34) closely resembles a Gauss’ Law. Sim-ilarly to Eq. (30), the total charge Q in the continuumis expressed in terms of the height fields as Q = (cid:90) L d x ρ ( x ) = φ (0) ( L ) − φ (0) (0) . (35)The preceding discussion illustrates how the globalcharge constraint on ρ ( x ) is re-expressed as a local bound-ary constraint on φ (0) ( x ), clarifying why the latter isphysically more appropriate as the field variable forcharge conserving systems. Similarly, using Eq. (33), wecan express Eq. (31) in terms of the generalized heightfields as ∂ x φ ( m ) ( x ) = φ ( m − ( x ) or ∂ m +1 x φ ( m ) ( x ) = ρ ( x ) . (36)In the continuum, the conservation of the m th and alllower moments then amounts to fixing the left andright boundary constraints on φ ( m ) ( x ) and its derivatives ∂ nx φ ( m ) ( x ) for all n ≤ m . For instance, the total dipolemoment and charge can be expressed in terms of bound-ary constraints on φ (1) ( x ) and ∂ x φ (1) ( x ) as Q = ∂ x φ (1) ( L ) − ∂ x φ (1) (0) ,Q = (cid:90) L d x x ρ ( x )= x φ (0) ( x ) (cid:12)(cid:12)(cid:12) L − (cid:90) L d x φ (0) ( x )= L ∂ x φ (1) ( L ) − ( φ (1) ( L ) − φ (1) (0)) , (37)where we have invoked Eqs. (34) and (36). In fact, it isstraightforward to show that the n th multipole momentcan be expressed in terms of the m th height field as ( m ≥ n ) Q n = (cid:90) L d x x n ρ ( x )= (cid:90) L d x x n ∂ m +1 x φ ( m ) ( x )= n ! n (cid:88) j =0 ( − j j ! (cid:16) L j ∂ m − n + jx φ ( m ) ( L ) − δ j, ∂ m − nx φ ( m ) (0) (cid:17) . (38)The multipole moments Q n in Eq. (38) are invariant un-der a polynomial shift of the height fields φ ( m ) , underwhich Eq. (36) is still satisfied: φ ( m ) ( x ) → φ ( m ) ( x ) + P ( m ) ( x ) , (39)where P ( m ) ( x ) is an arbitrary polynomial of degree ≤ m ;Eq. (39) can thus be used to set ∂ nx φ ( m ) (0) = 0 for all n ≤ m without loss of generality, so that Q n = n ! n (cid:88) j =0 ( − j L j j ! ∂ m − n + jx φ ( m ) ( L ) . (40)In the sector of primary interest, Q n = 0 for all n ≤ m ,these boundary constraints further simplify to ∂ nx φ ( m ) ( L ) = 0 ∀ n ≤ m . (41)For OBC, we need to further supplement these bound-ary constraints, which fix the symmetry sectors, with physical boundary conditions, which ensure that no mul-tipole currents flow through the boundaries. As discussedin Ref. [68], the fundamental hydrodynamic quantities formultipole conserving systems are the charge density ρ ( x )and the multipole current J ( m ) , from which one can inferthe conventional charge current; however, it is the mul-tipole current that is fundamental and is related to thecharge density as J ( m ) ( x ) ∼ ∂ m +1 x ρ ( x ) , (42)for systems which conserve all moments up to the m th highest moment, which is the generalization of Fick’slaw to multipole conserving systems [64]. The physicalrequirement that no multipole current flows through theboundaries, phrased in terms of the height fields, can bestated as ∂ m +1) x φ ( m ) (0) = ∂ m +1) x φ ( m ) ( L ) = 0 . (43) B. Ground state
Before we obtain the continuum limit of the Hamilto-nian H RK , we express the ground state wave function ofan m th moment conserving system, discussed in Sec. IV,in terms of the height field φ ( m ) ( x ). Recall that the GS2Eq. (23) is the equal-weight superposition of all allowedbasis states i.e., all possible height variable configura-tions that satisfy the boundary constraints, which fix thequantum number sectors.Treating the spin- s DOFs as “random variables” thatassume integer or half-integer values in [ − s, s ], undercoarse-graining the distribution of the spins flows to aGaussian as a direct consequence of the central limit the-orem [84]. The variance of the resulting coarse-grainedDOFs then scales as σ = ∆ x/κ where ∆ x is the lat-tice spacing and κ is a parameter chosen such that mi-croscopic correlation functions are accurately reproducedat long distances; effectively, one can think of κ as thecoarse-graining length scale. After coarse-graining, thewave functional Φ ( m )0 [ ρ ( x )] corresponding to a chargedensity profile ρ ( x ) is thus simply given by a Gaussian,albeit subject to global constraints specified by the con-served quantities (see Appendix B for details):Φ ( m )0 [ ρ ( x )] = 1 √Z e − κ (cid:82) L d x ( ρ ( x )) × G [ ρ ( x )] (44)where Z is a normalization constant and G [ ρ ( x )] en-forces the global symmetry constraints; namely, it fixesthe quantum number sector of interest. For instance, fora dipole conserving system in the { Q , Q } sector, G [ ρ ( x )] = δ (cid:18)(cid:90) ρ ( x ) − Q (cid:19) δ (cid:18)(cid:90) xρ ( x ) − Q (cid:19) . (45)To circumvent the global constraint in Eq. (44), itis convenient to work in terms of the m th height fieldsfor m th multipole conserving systems—as discussed inSec. V A, in this language, the quantum number sectorsare instead expressed as local boundary constraints. Moreexplicitly, Eq. (36) allows us to express the wave func-tional Eq. (44) in terms of the height field φ ( m ) ( x ) asΦ ( m )0 [ φ ( m ) ( x )] = 1 √Z e − κ (cid:82) L d x ( ∂ m +1 x φ ( m ) ( x )) B [ φ ( m ) ( x )] , (46)where the global constraints encoded in G [ ρ ( x )] are re-placed with local boundary constraints B [ φ ( m ) ( x )]. Theseconstraints are imposed by δ -functions that fix theboundary constraints on the height fields, correspondingto the quantum number sector of interest (see Eq. (38)).For the sector with Q n = 0 ∀ n ≤ m B [ φ ( m ) ( x )] = m (cid:89) n =0 δ (cid:16) ∂ nx φ ( m ) ( L ) (cid:17) , (47)which follows from Eq. (41). Note that we also need toimpose the physical boundary conditions Eq. (43) on thegeneralized height fields.Recall that in the discrete setting, the GS is an equalweight superposition of allowed configurations, while tak-ing the continuum limit introduces Gaussian weights intothe GS due to coarse-graining. Concurrently, the corre-sponding continuum Hamiltonian will no longer be of the form Eq. (22) but instead belongs to the class of “SMFdecomposable” Hamiltonians, that are related to clas-sical master equations and include the RK-HamiltoniansEq. (22) as a subclass [72]. This correspondence will proveuseful in deriving the continuum expression for the RK-Hamiltonian, which we discuss in Sec. V C (see also Ap-pendix C).To close this discussion, we note that expressions ofthe form Eq. (46) have previously been derived for thecontinuum limit of ground states of RK-Hamiltoniansusing various methods, albeit never in the context ofmultipole conserving systems. For instance, the expo-nent in Eq. (46) can be interpreted as the free energyfunctional corresponding to a configuration of the heightfield φ ( m ) ( x ), as is typically done in the context of RKpoints in dimer models [71, 83]. Alternately, the expres-sion Eq. (46) can also be derived using the path-integralformulation of Brownian motion: here, one interprets x as a time coordinate, ρ ( x ) in Eq. (36) as white-noise, and φ ( m ) ( x ) as a trajectory under the “Langevin dynamics”described by Eq. (36) [71, 84, 85]. C. Hamiltonian and dispersion relation
Having obtained the continuum expression for theground state wave functional, we now identify thecorresponding expression for the coarse-grained RK-Hamiltonian H RK . As discussed in the previous sec-tion, the coarse-grained wave functional Eq. (46) is theground state of a multipole conserving RK-Hamiltonian,which belongs to the generalized class of frustration-free positive-definite RK-like Hamiltonians discussed inRef. [72]. In Appendix D, we discuss two distinct ap-proaches for deriving the continuum limit of H RK : thefirst approach involves an appropriate choice of regula-tors, which allows us to explicitly obtain the continuumparent Hamiltonian corresponding to Eq. (46). The sec-ond, more commonly employed approach [71, 82, 83, 85,86] exploits the relationship between H RK and classicalmaster equations discussed in Sec. III (see Eq. (16)). Insummary, this approach proceeds by identifying the clas-sical process corresponding to H RK which equilibratesto a Gaussian distribution of height fields, as given byEq. (46). As we show in Appendix D 2, this classicalprocess describes the Langevin dynamics of the general-ized height fields under damping. The continuum expres-sion for H RK can then be derived via the Fokker-Planckequation for the probability functionals of the generalizedheight fields.Both approaches lead to the same continuum expres-sion for H RK , which is the parent Hamiltonian for the Real, symmetric, and irreducible matrices which admit aStochastic Matrix Form (SMF) decomposition were found to bein 1-to-1 correspondence with classical stochastic systems de-scribed by a master equation in Ref. [72] H ( m ) = γ (cid:90) L d x Q † m ( x ) Q m ( x ) , (48)where γ is an overall dimensionful constant. The creationand annihilation operators Q † m ( x ) and Q m ( x ) are definedas Q † m ( x ) = 1 √ (cid:18) δδφ ( m ) + ( − m κ ∂ m +1) x φ ( m ) (cid:19) Q m ( x ) = 1 √ (cid:18) − δδφ ( m ) + ( − m κ ∂ m +1) x φ ( m ) (cid:19) , (49)and satisfy the commutation relations[ Q m ( x ) , Q † m ( y )] = ( − m +1 κ ∂ m +1) x δ ( x − y ) . (50)We can directly verify that the wave functionalΦ ( m )0 [ φ ( m ) ] Eq. (46) is a “frustration-free” ground stateof the Hamiltonian H ( m ) Eq. (48) by noting that (seeEqs. (D1) and (D11)) δδφ ( m ) Φ ( m )0 [ φ ( m ) ] = ( − m κ (cid:16) ∂ m +1) x φ ( m ) (cid:17) Φ ( m )0 [ φ ( m ) ] , (51)resulting in Q m ( x )Φ ( m )0 [ φ ( m ) ] = 0 ∀ x. (52)Up to a constant (infinite) energy shift, the continuumHamiltonian Eq. (48) can be brought to more standardform [82, 85] H ( m ) = γ (cid:90) L d x (cid:20) (cid:16) Π ( m ) (cid:17) + κ (cid:16) ∂ m +1) x φ ( m ) (cid:17) (cid:21) , (53)where Π ( m ) ( x ) = iδ/δφ ( m ) ( x ) is the canonical momen-tum which satisfies [ φ ( m ) ( x ) , Π ( m ) ( y )] = iδ ( x − y ).We observe that the Hamiltonian Eq. (53) is invariantunder a polynomial shift of the form φ ( m ) ( x ) → φ ( m ) ( x ) + P (2 m +1) ( x ) , (54)where P (2 m +1) ( x ) is a polynomial in x of degree ≤ (2 m + 1). That is, it has additional symmetries beyondjust the m th multipole moment, as is typical of contin-uum RK-Hamiltonians [81]. However, using Eq. (38), itis straightforward to see that the transformation Eq. (54)changes the quantum number sector of the system. Thisshows that the continuum Hamiltonian Eq. (53) is thesame across all quantum number sectors, further imply-ing that the ground state sector is extensively degenerate.Now that we have established the form of the contin-uum Hamiltonian, we can study its lowest energy excitedstates to derive the dispersion relation and the gap. Us-ing Eqs. (48) and (50), the excited states Φ k [ φ ( m ) ] canbe written asΦ k [ φ ( m ) ( x )] = (cid:90) L d x f ( m ) ( kx ) Q † m ( x )Φ[ φ ( m ) ( x )] , (55) where the mode function f ( m ) ( kx ) is determined by theboundary constraints on the height field φ ( m ) ( x ), where k is the momentum of the mode. For large system sizes,we expect that deep within the bulk f ( m ) ( kx ) ∼ e ikx [86]from which we obtain the dispersion relation H ( m ) Φ k [ φ ( m ) ] = γκk m +1) Φ k [ φ ( m ) ] . (56)For a finite system of size L , we thus expect the gap ∆ ( m ) of H ( m ) to scale as ∆ ( m ) ∼ L m +1) . (57)For charge-conserving systems ( m = 0), we see that thecontinuum height field approach correctly reproduces thescaling of the spin-wave dispersion relation of the Heisen-berg model discussed in Sec. IV. More generally, we canfurther lower-bound the scaling of the Thouless time t ( m )Th for a system of size L conserving the m th multipole mo-ment as follows: t ( m )Th (cid:38) L m +1) . (58)Due to the polynomial shift symmetry (Eq. (54)) of thecontinuum Hamiltonian, we expect that this scaling ofthe Thouless time is independent of the quantum num-ber sector. Eq. (58) is one of the main results of this pa-per as it encodes the subdiffusive scaling of the Thoulesstime in systems with higher moment conservation laws.These results, obtained analytically through the general-ized height representation developed herein, are validatedby the numerical analysis performed on dipole conservingFRQCs (see Sec. IV).The applicability of our continuum analysis of H RK extends beyond the context of random quantum circuitsand is directly pertinent to the study of classical cellu-lar automata with conserved higher moments. Such au-tomata were studied in Refs. [63, 64] and are equivalentto the circuit (cid:99) M . To further test the validity of the contin-uum Hamiltonian obtained in Eq. (53), we can computethe two-point spin correlations using Eq. (36): (cid:104) ρ ( x, t ) ρ (0 , (cid:105) H ∝ κt ) / m +1) F (cid:18) x m +1) κt (cid:19) , (59)where ρ ( x, t ) = e − iH ( m ) t ρ ( x ) e iH ( m ) t , (cid:104)·(cid:105) H = (cid:82) D φ ( m ) ( · ) exp( − i (cid:82) dx dt H ( m ) [ φ ( m ) ]), and F isa hypergeometric scaling function. Eq. (59) is in agree-ment with scalings obtained from numerical calculationsand hydrodynamic considerations in Ref. [64]. VI. CONCLUDING REMARKS
In this paper, we have studied the spectral statistics,as encoded in the SFF K ( t ), for spatially-extended con-strained many-body quantum chaotic systems, focusingon FRQCs with conserved higher moments, such as the4dipole moment. As one of the key results of this paper,we have established a series of relations between K ( t ) inthe q → ∞ limit, a classical stochastic circuit (cid:99) M , andan emergent RK-Hamiltonian, such that the inverse gapof this RK-Hamiltonian lower bounds the Thouless time t Th of the underlying FRQC. As we have shown here, therelation between t Th and ∆ RK proves particularly effica-cious, since it relates a dynamical property of the FRQCto the low-energy physics of a sign-problem-free quantumHamiltonian.We emphasise that these relations are valid for generic local FRQCs with on-site Abelian symmetries or dynami-cal constraints, not only those with conserved higher mo-ments of charge. For example, we can consider circuit im-plementations of other fragmented models [87] or studyan FRQC inspired by the Rydberg blockade [56, 88], alsoknown as the PXP model [8]. The latter is implementedby taking e.g., (cid:96) = 3 site local gates with the only non-trivial dynamics contained within a 2 × |↓↓↓(cid:105) and |↓↑↓(cid:105) states. The resulting Floquet oper-ator (cid:99) W has no conserved quantities besides the (quasi)-energy, but fragments into dynamically disconnected sub-spaces; the largest of these corresponds to the constrainedHilbert space most often discussed in the context of quan-tum many-body scar dynamics [8].We have verified these general results on circuits withhigher conserved moments, which generically exhibitHilbert space fragmentation. Working in the q → ∞ limit, we derived the corresponding stochastic circuit (cid:99) M and emergent RK-Hamiltonian H RK for both charge and(weakly fragmented) dipole conserving systems. Our nu-merical study of these systems suggests a universalityin the scaling of t Th with system size, specifically, wepredict diffusive scaling t Th ∼ L for charge conserv-ing systems and subdiffusive behavior ∼ L for dipoleconserving systems, regardless of the microscopic detailsof the underlying circuit. Further evidence for this scal-ing is given by continuum field theoretic descriptions ofthe emergent RK-Hamiltonians for multipole conservingFRQCs in terms of generalized height fields. By analyti-cally computing the dispersion relation for the resultantfield theories, we find that t Th ∼ L m +1) in a circuitsthat conserve the m th multipole moment, consistent withnumerical results for charge and dipole conserving sys-tems.Our work opens many exciting avenues for future re-search: here, we have only focused on the class of mul-tipole conserving circuits which exhibit weak fragmen-tation. Dynamics in strongly fragmented systems, wheretypical initial states are ETH-violating, is highly con-strained; nevertheless, such systems exhibit large Krylov subspaces which eventually thermalize [16]. The scalingof t Th within such subspaces remains to be understoodand may lead to distinct continuum field theories thanthose we have introduced for weakly fragmented sys-tems. Another interesting avenue to explore is extendingour formalism to incorporate non-Abelian symmetries,for which the nature of transport and thermalization iscurrently being debated [89–91].We note that the large- q diagrammatics, and there-fore the mapping to a classical bistochastic circuit andRK-Hamiltonian, have so far only been developed for thetwo-point SFF K ( t ). Other observables, such as the sec-ond Renyi entropy and out-of-time-order correlator, canalso be mapped to stochastic classical dynamics uponensemble averaging, and will be discussed in forthcomingwork. Pushing these ideas further presents an importantbut technically-demanding theoretical challenge. Morestraightforward is extending our results to circuit geome-tries besides the brick-wall structure considered here aswell as to other RMT symmetry classes.More pressing, however, is building a systematic under-standing of FRQCs at finite- q , to delineate those featureswhich are an artefact of the q → ∞ limit from thosewhich are more generic properties of constrained ran-dom circuits. Numerically investigating finite- q circuitsremains prohibitive, particularly in the context of highermoment conserving circuits which already require large( (cid:96) ≥
4) local gates. Analytically, one could attempt tokeep track of diagrams at next to leading order in thelarge- q expansion to better quantify deviations of theSFF from the strict q → ∞ limit. We leave the devel-opment of such analytical techniques to future work. ACKNOWLEDGMENTS
We are particularly grateful to Shivaji Sondhi for en-lightening discussions. We also acknowledge useful con-versations with Nathan Benjamin, John Chalker, AndreaDe Luca, and Alan Morningstar. A. P. was supported inpart with funding from the Defense Advanced ResearchProjects Agency (DARPA) via the DRINQS program.The views, opinions and/or findings expressed are thoseof the authors and should not be interpreted as repre-senting the official views or policies of the Departmentof Defense or the U.S. Government. A. C. is supportedin part by the Croucher foundation. A. P. and A. C.are supported by fellowships at the PCTS at PrincetonUniversity. D. A. H. is supported in part by DOE grantDE-SC0016244. [1] I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. , 885 (2008).[2] I. V. Gornyi, A. D. Mirlin, and D. G. Polyakov, Phys. Rev. Lett. , 206603 (2005).[3] D. M. Basko, I. L. Aleiner, and B. L. Altshuler, Annalsof Physics , 1126 (2006). [4] A. Pal and D. A. Huse, Phys. Rev. B , 174411 (2010).[5] R. Nandkishore and D. A. Huse, Annual Review of Con-densed Matter Physics , 15 (2015).[6] N. Shiraishi and T. Mori, Phys. Rev. Lett. , 030601(2017).[7] S. Moudgalya, S. Rachel, B. A. Bernevig, and N. Reg-nault, Phys. Rev. B , 235155 (2018).[8] C. J. Turner, A. A. Michailidis, D. A. Abanin, M. Serbyn,and Z. Papi, Nature Physics , 745 (2018).[9] S. Moudgalya, N. Regnault, and B. A. Bernevig, Phys.Rev. B , 235156 (2018).[10] C. J. Turner, A. A. Michailidis, D. A. Abanin, M. Serbyn,and Z. Papi´c, Phys. Rev. B , 155134 (2018).[11] W. W. Ho, S. Choi, H. Pichler, and M. D. Lukin, Phys.Rev. Lett. , 040603 (2019).[12] V. Khemani, C. R. Laumann, and A. Chandran, PhysicalReview B , 161101 (2019).[13] P. Sala, T. Rakovszky, R. Verresen, M. Knap, andF. Pollmann, Phys. Rev. X , 011047 (2020).[14] V. Khemani, M. Hermele, and R. Nandkishore, Phys.Rev. B , 174204 (2020).[15] T. Rakovszky, P. Sala, R. Verresen, M. Knap, andF. Pollmann, Phys. Rev. B , 125126 (2020).[16] S. Moudgalya, A. Prem, R. Nandkishore, N. Reg-nault, and B. A. Bernevig, arXiv e-prints (2019),arXiv:1910.14048 [cond-mat.str-el].[17] J. M. Deutsch, Phys. Rev. A , 2046 (1991).[18] M. Srednicki, Phys. Rev. E , 888 (1994).[19] M. Rigol, V. Dunjko, and M. Olshanii, Nature , 854(2008).[20] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalat-tore, Rev. Mod. Phys. , 863 (2011).[21] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol,Advances in Physics , 239 (2016).[22] J. Kudler-Flam, L. Nie, and S. Ryu, JHEP , 175(2020).[23] A. Kitaev, in KITP strings seminar and Entanglement,Vol. 12 (2015).[24] J. Maldacena, S. H. Shenker, and D. Stanford, JHEP , 106 (2016).[25] J. Maldacena and D. Stanford, Phys. Rev. D , 106002(2016).[26] J. S. Cotler, G. Gur-Ari, M. Hanada, J. Polchinski,P. Saad, S. H. Shenker, D. Stanford, A. Streicher, andM. Tezuka, JHEP , 118 (2017).[27] J. Cotler, N. Hunter-Jones, J. Liu, and B. Yoshida, JHEP , 48 (2017).[28] A. Nahum, J. Ruhman, S. Vijay, and J. Haah, Phys.Rev. X , 031016 (2017).[29] C. W. von Keyserlingk, T. Rakovszky, F. Pollmann, andS. L. Sondhi, Phys. Rev. X , 021013 (2018).[30] A. Hamma, S. Santra, and P. Zanardi, Phys. Rev. Lett. , 040502 (2012).[31] A. Hamma, S. Santra, and P. Zanardi, Phys. Rev. A ,052324 (2012).[32] P. Kos, M. Ljubotina, and T. Prosen, Phys. Rev. X ,021062 (2018).[33] H. Gharibyan, M. Hanada, S. H. Shenker, andM. Tezuka, JHEP , 124 (2018).[34] S. Moudgalya, T. Devakul, C. W. von Keyserlingk, andS. L. Sondhi, Phys. Rev. B , 094312 (2019).[35] C. Jonay, D. A. Huse, and A. Nahum, arXiv e-prints(2018), arXiv:1803.00089 [cond-mat.stat-mech].[36] V. Khemani, A. Vishwanath, and D. A. Huse, Phys. Rev. X , 031057 (2018).[37] T. Rakovszky, F. Pollmann, and C. W. von Keyserlingk,Phys. Rev. X , 031058 (2018).[38] O. Bohigas, M. J. Giannoni, and C. Schmit, Phys. Rev.Lett. , 1 (1984).[39] G. Montambaux, D. Poilblanc, J. Bellissard, and C. Sire,Phys. Rev. Lett. , 497 (1993).[40] D. Poilblanc, T. Ziman, J. Bellissard, F. Mila, andG. Montambaux, Europhysics Letters (EPL) , 537(1993).[41] F. Haake, “Quantum signatures of chaos,” inQuantum Coherence in Mesoscopic Systems, editedby B. Kramer (Springer US, Boston, MA, 1991) pp.583–595.[42] D. J. Thouless, Phys. Rev. Lett. , 1167 (1977).[43] F. J. Dyson, Journal of Mathematical Physics , 1191(1962).[44] B. Bertini, P. Kos, and T. Prosen, Phys. Rev. Lett. ,264101 (2018).[45] A. Flack, B. Bertini, and T. Prosen, “Statistics of thespectral form factor in the self-dual kicked ising model,”(2020), arXiv:2009.03199 [nlin.CD].[46] A. Chan, A. De Luca, and J. T. Chalker, Phys. Rev. X , 041019 (2018).[47] A. Chan, A. De Luca, and J. T. Chalker, Phys. Rev.Lett. , 060601 (2018).[48] A. J. Friedman, A. Chan, A. De Luca, and J. T. Chalker,Phys. Rev. Lett. , 210603 (2019).[49] A. Chan, A. De Luca, and J. T. Chalker, Phys. Rev.Lett. , 220601 (2019).[50] S. J. Garratt and J. T. Chalker, arXiv e-prints (2020),arXiv:2008.01697 [cond-mat.stat-mech].[51] C. Chamon, Phys. Rev. Lett. , 040402 (2005).[52] I. H. Kim and J. Haah, Phys. Rev. Lett. , 027202(2016).[53] K. Siva and B. Yoshida, Phys. Rev. A , 032324 (2017).[54] A. Prem, J. Haah, and R. Nandkishore, Phys. Rev. B , 155133 (2017).[55] A. Chandran, M. D. Schulz, and F. J. Burnell, Phys.Rev. B , 235122 (2016).[56] H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Om-ran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres,M. Greiner, V. Vuleti´c, and M. D. Lukin, Nature ,579 (2017).[57] J. Feldmeier, F. Pollmann, and M. Knap, Phys. Rev.Lett. , 040601 (2019).[58] J. Iaconis, S. Vijay, and R. Nandkishore, Phys. Rev. B , 214301 (2019).[59] F. Ritort and P. Sollich, Advances in Physics , 219(2003).[60] S. Pai, M. Pretko, and R. M. Nandkishore, Physical Re-view X , 021003 (2019).[61] S. Pai and M. Pretko, Phys. Rev. Lett. , 136401(2019).[62] S. Moudgalya, B. A. Bernevig, and N. Regnault, arXive-prints (2019), arXiv:1906.05292 [cond-mat.str-el].[63] A. Morningstar, V. Khemani, and D. A. Huse, Phys.Rev. B , 214205 (2020).[64] J. Feldmeier, P. Sala, G. de Tomasi, F. Pollmann, andM. Knap, arXiv e-prints (2020), arXiv:2004.00635 [cond-mat.str-el].[65] E. Guardado-Sanchez, A. Morningstar, B. M. Spar, P. T.Brown, D. A. Huse, and W. S. Bakr, Phys. Rev. X ,011042 (2020). [66] R. M. Nandkishore and M. Hermele, Annual Review ofCondensed Matter Physics , 295 (2019).[67] M. Pretko, X. Chen, and Y. You, International Journalof Modern Physics A , 2030003 (2020).[68] A. Gromov, A. Lucas, and R. M. Nandkishore, Phys.Rev. Research , 033124 (2020).[69] P. Zhang, Phys. Rev. Research , 033129 (2020).[70] D. S. Rokhsar and S. A. Kivelson, Phys. Rev. Lett. ,2376 (1988).[71] C. L. Henley, Journal of Statistical Physics , 483(1997).[72] C. Castelnovo, C. Chamon, C. Mudry, and P. Pujol,Annals of Physics , 316 (2005).[73] R. D. Somma, C. D. Batista, and G. Ortiz, Phys. Rev.Lett. , 030603 (2007).[74] C. L. Henley, Journal of Physics: Condensed Matter ,S891 (2004).[75] S. Knabe, Journal of Statistical Physics , 627 (1988).[76] D. Gosset and E. Mozgunov, Journal of MathematicalPhysics , 091901 (2016).[77] M. Lemm and E. Mozgunov, Journal of MathematicalPhysics , 051901 (2019).[78] S. Bravyi, Quantum Info. Comput. , 11221140 (2015).[79] F. R. Chung and F. C. Graham, Spectral graph theory,92 (American Mathematical Soc., 1997).[80] D. Roy and T. Prosen, arXiv e-prints (2020),arXiv:2005.10489 [cond-mat.stat-mech].[81] E. Fradkin, D. A. Huse, R. Moessner, V. Oganesyan, andS. L. Sondhi, Phys. Rev. B , 224415 (2004).[82] E. Ardonne, P. Fendley, and E. Fradkin, Annals ofPhysics , 493 (2004).[83] R. Moessner and K. S. Raman, inIntroduction to Frustrated Magnetism (Springer, 2011)pp. 437–479.[84] S. N. Majumdar, “Brownian functionals in physics andcomputer science,” in The Legacy of Albert Einstein, pp.93–129.[85] X. Chen, E. Fradkin, and W. Witczak-Krempa, Phys.Rev. B , 180402 (2017).[86] X. Chen, E. Fradkin, and W. Witczak-Krempa, Journalof Physics A: Mathematical and Theoretical , 464002(2017).[87] Z.-C. Yang, F. Liu, A. V. Gorshkov, and T. Iadecola,Phys. Rev. Lett. , 207602 (2020).[88] I. Lesanovsky, Physical Review Letters , 025301(2011).[89] I. V. Protopopov, R. K. Panda, T. Parolini, A. Scardic-chio, E. Demler, and D. A. Abanin, Phys. Rev. X ,011025 (2020).[90] Z.-C. Yang, S. Nicholls, and M. Cheng, arXiv e-prints(2020), arXiv:2007.02950 [cond-mat.dis-nn].[91] P. Glorioso, L. V. Delacr´etaz, X. Chen, R. M. Nandkishore, and A. Lucas, arXiv e-prints (2020),arXiv:2007.13753 [cond-mat.stat-mech].[92] P. M. Chaikin, T. C. Lubensky, and T. A. Witten,Principles of condensed matter physics, Vol. 10 (Cam-bridge university press Cambridge, 1995).[93] H. Risken, “Fokker-planck equation,” inThe Fokker-Planck Equation: Methods of Solution and Applications(Springer Berlin Heidelberg, Berlin, Heidelberg, 1996)pp. 63–95. Appendix A: Mapping K ( t ) to a classical Markov circuit In this Appendix, we generalize the calculation of SFF K ( t ; K ) performed in Ref. [48] for an FRQC with conservedU(1) charge to FRQCs with arbitrary symmetries or local constraints, such as dipole moment conservation. Specifically,we show that in the q → ∞ limit, K ( t ; K ) is mapped to the trace of the t -th power of a bistochastic matrix.The ensemble average in K ( t ; K ), defined in Eq. (6) and illustrated in Fig. 6(a), can be evaluated as a sum ofdiagrams, each of which corresponds to a possible pairing between unitaries and their complex conjugates. In thelimit of large- q , diagrams with the same “cyclical” pairing at all sites dominate the sum (see illustration in Fig. 6(c)and Ref. [46] for details). There are t such diagrams, which we call “Gaussian” since they can be evaluated using theWick contractions between unitaries and their conjugates. For (cid:96) -site unitaries U and their conjugates U ∗ , shown inFig. 6(d), we have (cid:68) U (cid:126)j,(cid:126)β(cid:126)i,(cid:126)α U ∗ (cid:126)j (cid:48) ,(cid:126)β (cid:48) (cid:126)i (cid:48) ,(cid:126)α (cid:48) (cid:69) = 1 d (cid:126)α q (cid:96) δ (cid:126)i,(cid:126)i (cid:48) δ (cid:126)j,(cid:126)j (cid:48) δ (cid:126)α,(cid:126)α (cid:48) δ (cid:126)β,(cid:126)β (cid:48) , δ (cid:126)i,(cid:126)i (cid:48) = (cid:96) (cid:89) a =1 δ i a ,i (cid:48) a , (A1)where (cid:126)i = ( i , i , . . . , i (cid:96) ) and (cid:126)α = ( α , α , . . . , α (cid:96) ) are (cid:96) -component vectors denoting respectively the color and spindegrees of freedom on their support. Here, d (cid:126)α denotes the number of “out-going” (cid:96) -site spin configurations dynamicallyaccessible to the “in-coming” (cid:96) -site spin configuration (cid:126)α . That is, d (cid:126)α q (cid:96) is the size of the random matrix block in (cid:99) W to which the basis state ( (cid:126)i, (cid:126)α ) belongs.As an example, let us consider the simplest charge-conserving circuit [48]. We have two-site gates ( (cid:96) = 2) and (cid:126)α = ( α , α ), with α i = ↑ , ↓ . The dynamically connected two-site configurations can be labelled by a conservedcharge Q = (cid:80) (cid:96)x =1 S zx . The configurations ( ↑ , ↑ ) and ( ↓ , ↓ ) have Q ( (cid:126)α ) = +1 and Q ( (cid:126)α ) = − d (cid:126)α = 1. However, the FIG. 6. (a) A schematic illustration of the diagrammatic representation of K ( t ; K ) (adapted from Ref. [46]) for a FRQC with (cid:96) = 4 (which is appropriate for the minimal spin-1 dipole-conserving FRQC with weak fragmentation). The white and greysheets represent Tr K [ (cid:99) W ( t )] and Tr K [ (cid:99) W † ( t )] respectively. Space runs horizontally and time runs vertically. Curly lines on the topand the bottom of the sheet represents traces over the dof at each site. Gates with different support are denoted with differentcolors. (b) The diagrammatic representation of Tr K [ (cid:99) M t ] where, unlike (cid:99) W , (cid:99) M is a non-random circuit that acts only on the spindegrees of freedom of the FRQC (cid:99) W . Each (cid:96) -gate in (cid:99) M is a block-diagonal bistochastic matrix defined in Eq. (9). All charges arepreserved after the action of every (cid:96) -site gates in (cid:99) M . (c) The three leading diagrams in the diagrammatic expansion of K ( t ; S )at t = 3 in the limit of large- q [46]. Note that for each diagram, every site i takes the same configuration. The grey ribbonscorrespond to a ladder of (cid:96) number of contractions between unitaries and their conjugates, represented by colored dots. (d)The diagrammatic representation of U (cid:126)j,(cid:126)β(cid:126)i,(cid:126)α , where the Roman (Greek) indices correspond to color (spin) degrees of freedom. ↑ , ↓ ) and ( ↓ , ↑ ) both have Q ( α ) = 0, and are dynamically connected to each other. Thus we have d (cid:126)α = 2.Next, we can translate each of the t diagrams into algebraic terms by using Eq. (A1) and summing over the colordegrees of freedom. Consider the diagram where the pairing between unitaries and their conjugates takes the formof the leftmost diagram in Fig. 6(c) at every site. The sum over the colors precisely cancels out all factors of q inEq. (A1) [46]. It then remains to sum over the spin degrees of freedom. Observe that the choice of spin DOFs in (cid:99) W uniquely fixes those in (cid:99) W † due to the Wick contractions Eq. (A1). Consequently, the sum over spins can be computedby finding all possible ways of assigning spins in the diagrammatic representation of Tr K (cid:104)(cid:99) W ( t ) (cid:105) (Fig. 6(b), suchthat all charges { Q a } are preserved after the action of every (cid:96) -site gate. This sum over spins can be reproduced byTr K (cid:104) (cid:99) M t (cid:105) in Eq. (10), where (cid:99) M has the geometry of (cid:99) W in Eq. (2) and contains non-random block-diagonal (cid:96) -gates (cid:98) m , defined in Eq. (9). Note that (cid:98) m contains d (cid:126)α × d (cid:126)α matrix blocks (cid:98) m ( d (cid:126)α ) with entries 1 /d (cid:126)α to account for the factor d (cid:126)α in Eq. (A1).The overall factor of | t | in Eq. (7) arises from the following argument: any two of the t contracted diagrams inFig. 6(c) are related by a rotation of the arrowed loop on the right side of one of the diagrams i.e., the diagrams aretopologically equivalent. As a result, the corresponding algebraic terms for all t leading diagrams are identical. Thisconcludes the derivation of Eq. (7). Appendix B: Continuum Limit of the RK Ground State wave function
In this Appendix, we derive the continuum limit for the ground state wave function of multipole conserving RK-Hamiltonians H RK , thereby allowing us to analyze these systems using field-theoretic techniques. In the discretesetting, the ground state wave function Eq. (23) has the form | ψ (cid:105) = (cid:88) { s n }∈ Λ |{ s n }(cid:105) , (B1)where Λ denotes the set of allowed spin configurations that are allowed within the quantum number sector or Krylovsubspace of interest.As discussed in Sec. V A, the spins { s n } assume integer or half-integer values in [ − s, s ]. In the continuum limit i.e.,when the lattice spacing ∆ x →
0, under coarse-graining the distributions of s n ’s flows to Gaussian random variableswith zero mean and variance σ ∼ ∆ x as a consequence of the central limit theorem [84]. That is, the probabilitydensity P of a single spin s n is P ( s n ) ∼ exp (cid:18) − s n σ (cid:19) = ⇒ (cid:104) s n (cid:105) = 0 , (cid:104) s n s n (cid:48) (cid:105) = σ δ n,n (cid:48) , σ = ∆ xκ . (B2)Here, κ effectively serves as the coarse-graining length scale and is determined by requiring that the continuum theorycorrectly reproduces long-distance correlation functions in the microscopic model, which can in principle be computednumerically. Next, by expressing the s n ’s in terms of the m th height fields and by using Eqs. (29) and (31), we obtain s n = φ (0) n + − φ (0) n − = φ (1) n +1 − φ (1) n + φ (1) n − = · · · = m +12 (cid:88) j = − m +12 ( − j + m +12 (cid:18) m + 1 j + m +12 (cid:19) φ ( m ) n + j , (B3)while the joint probabilities satisfy P ( { s n } ) ∼ exp (cid:32) − σ (cid:88) n s n (cid:33) × G [ { s n } ]= ⇒ P ( { φ ( m ) n } ) ∼ exp − σ (cid:88) n m +12 (cid:88) j = − m +12 ( − j + m +12 (cid:18) m + 1 j + m +12 (cid:19) φ ( m ) n + j × B [ { φ ( m ) n } ] , (B4)where σ = ∆ x/κ . G [ { s n } ] and B [ { φ ( m ) n } ] are shorthand for global and boundary constraints written in terms of spinand height variables respectively; these constraints are required to fix the quantum number sector, as discussed in9Sec. V A. Note that the Jacobian of the transformation from { s n } → { φ ( m ) n } is 1. In order to obtain a continuumlimit in the case of m th height variables, we need to define height fields that scale with the lattice spacing as shownin Eq. (33). The probability density for a configuration of the m th height field in the continuum limit then reads P ( φ ( m ) ( x )) ∼ lim ∆ x → exp − κ (cid:88) n ∆ x m +12 (cid:80) j = − m +12 ( − j + m +12 (cid:0) m +1 j + m +12 (cid:1) φ ( m ) n + j (∆ x ) m (∆ x ) m +1 × B [ { φ ( m ) n } ]= exp (cid:32) − κ (cid:90) L dx ( ∂ m +1 x φ ( m ) ( x )) (cid:33) × B [ { φ ( m ) ( x ) } ] , (B5)from which we obtain Eq. (46) as the continuum limit of the ground state wave functional. The full continuum wavefunction can be written as (cid:12)(cid:12)(cid:12) Φ ( m )0 (cid:69) = 1 √Z (cid:90) B [ φ ( m ) ( x )] D φ ( m ) exp (cid:32) − κ (cid:90) L dx ( ∂ m +1 x φ ( m ) ( x )) (cid:33) (cid:12)(cid:12)(cid:12) φ ( m ) ( x ) (cid:69) . (B6) Appendix C: Brief Interlude on SMF Decomposable Hamiltonians
In this Appendix, we discuss a generalization of the RK wave functions in Eq. (23) to states of the form | ψ SMF (cid:105) = 1 √Z (cid:88) C∈ Λ p C |C(cid:105) , (C1)where the Hilbert space is spanned by basis states composed of an abstract set of configurations Λ = {C} , p C is the“probability” associated with the configuration C , and √Z is the normalization constant for the wave function. Werecover the standard RK wave functions of the form Eq. (23) by setting p C = 1 for all configurations C .As discussed in Refs. [72, 74], wave functions of the form Eq. (C1) are closely related to steady states of Markovprocesses satisfying detailed balance. Moreover, these generalized RK wave functions correspond to ground states ofHamiltonians which admit a “Stochastic Matrix Form” decomposition and are of the following form: H SMF = (cid:80) (cid:104)C , C (cid:48) (cid:105) w C , C (cid:48) (cid:98) Q C , C (cid:48) , (cid:98) Q C , C (cid:48) = (cid:104) p C |C(cid:105) (cid:104)C| + p C(cid:48) |C (cid:48) (cid:105) (cid:104)C (cid:48) | − p C p C(cid:48) ( |C(cid:105) (cid:104)C (cid:48) | + |C (cid:48) (cid:105) (cid:104)C| ) (cid:105) = (cid:16) p C |C(cid:105) − p C(cid:48) |C (cid:48) (cid:105) (cid:17) (cid:16) p C (cid:104)C| − p C(cid:48) (cid:104)C (cid:48) | (cid:17) , (C2)where (cid:104)C , C (cid:48) (cid:105) represents a pair of configurations from the set Λ that are connected under the action of the Hamiltonian.In Eq. (C2), w C , C (cid:48) ’s are positive real numbers to ensure that H SMF ≥
0. It is straightforward to verify by directcomputation that | ψ SMF (cid:105) is a frustration-free ground state of H SMF since (cid:98) Q C , C (cid:48) ( p C |C(cid:105) + p C (cid:48) |C (cid:48) (cid:105) ) = 0 = ⇒ H SMF | ψ SMF (cid:105) = 0 . (C3)Typically, p C and Z are represented as [72] p C = e − βE C , Z = (cid:88) C∈ Λ e − βE C . (C4)where E C is the “energy” associated with the configurations, β is the “inverse-temperature” and hence the Z resemblesthe partition function. We note that while the standard RK-Hamiltonian H RK Eq. (22) is proportional to the transitionmatrix of a Markov process, H SMF is related to a transition matrix T of the Markov process through a similaritytransformation, under which the ground state wave function remains of the form Eq. (C1). Appendix D: Continuum Limit of the RK Hamiltonian
In this Appendix, we discuss the derivation of the continuum limit of the RK-Hamiltonian for multipole conservingsystems. We outline two distinct approaches, both of which lead to the same Hamiltonian.0
1. Regulator approach
Let us begin by observing that the continuum wave function Eq. (B6) has the form Eq. (C1) i.e., it belongs tothe class of ground states of SMF decomposable Hamiltonians. We thus anticipate that the continuum limit of theRK-Hamiltonian is SMF decomposable and can be brought to the form Eq. (C2).To make this correspondence precise, we identify the set of configurations Λ with the set of height field configurations { φ ( m ) ( x ) } , while the inverse-temperature β , energy { E ( φ ( m ) ( x )) } , and partition function Z are respectively given by(see Eq. (C4)) β = κ, E ( φ ( m ) ( x )) = (cid:90) L d x (cid:16) ∂ m +1 x φ ( m ) ( x ) (cid:17) , Z = (cid:90) B [ φ ( m ) ( x )] D φ ( m ) exp (cid:16) − βE ( φ ( m ) ( x )) (cid:17) . (D1) B [ φ ( m ) ( x )] is shorthand for the boundary constraints Eq. (38) that fix the quantum number sector. Note that inthe case of systems involving (cid:96) -site gates, two configurations of spins (height variables) C = { s n } ( C = { φ ( m ) n } ) and C (cid:48) = { s (cid:48) n } ( C (cid:48) = { φ ( m ) n (cid:48) } ) are connected under the action of the Hamiltonian only if the (cid:96) local spins ( (cid:96) − s j + s j +1 = s (cid:48) j + s (cid:48) j +1 for some j and s n = s (cid:48) n ∀ n, n (cid:54) = j, j + 1 ⇐⇒ φ (0) j + = φ (0) j + (cid:48) + s j − s (cid:48) j for some j and φ (0) n + = φ (0) n + (cid:48) ∀ n, n (cid:54) = j . (D2)We see that the Hamiltonian remains local when switching to the height representation, since the action of theHamiltonian changes the height variable only at a single point. While this is strictly true only for multipole conservingprocesses of the minimal gate-size, after coarse-graining we expect that the long-wavelength physics is captured by aHamiltonian that only permits local fluctuations in the height field configurations.Consequently, after coarse-graining we can express the RK-Hamiltonian Eq. (22) for an m th multipole momentconserving system in the SMF-form Eq. (C2): H ( m ) = (cid:88) { φ ( m ) n } , { φ ( m ) n (cid:48) } (cid:88) j (cid:89) n (cid:54) = j δ φ ( m ) n ,φ ( m ) n (cid:48) (cid:98) Q { φ ( m ) n } , { φ ( m ) n (cid:48) } w { φ ( m ) n } , { φ ( m ) n (cid:48) } . (D3)Since in taking the continuum limit, we replace discrete spin variables by continuous Gaussian random variables withvariance σ ∼ ∆ x , we replace sums in Eq. (D3) by integrals as follows (cid:88) { φ ( m ) n } , { φ ( m ) n (cid:48) } → (cid:90) (cid:89) n dφ ( m ) n dφ ( m ) n (cid:48) P ( { φ ( m ) n } ) P ( { φ ( m ) n (cid:48) } ) , (D4)where P denotes the Gaussian probability density Eq. (B4). Discrete delta functions are replaced by continuous ones.We then obtain the following expression for the continuum Hamiltonian H ( m ) H ( m ) = (cid:88) j (cid:90) (cid:89) n (cid:54) = j dφ ( m ) n dφ ( m ) n (cid:48) δ ( φ ( m ) n − φ ( m ) n (cid:48) ) (cid:90) dφ ( m ) j dφ ( m ) j (cid:48) P ( { φ ( m ) n } ) P ( { φ ( m ) n (cid:48) } ) w { φ ( m ) n } , { φ ( m ) n (cid:48) } (cid:98) Q { φ ( m ) n } , { φ ( m ) n (cid:48) } . (D5)It is clear that with an appropriate choice of the regulator w { φ ( m ) n } , { φ ( m ) n (cid:48) } , we can write Eq. (D5) in the form H ( m ) = (cid:90) (cid:89) n dφ ( m ) n P ( { φ ( m ) n } ) (cid:88) j (cid:90) dλ j exp (cid:32) − λ j σ (cid:33) (cid:98) Q { φ ( m ) n } , { φ ( m ) n (cid:48) } , λ j ≡ φ (cid:48) j − φ j , (D6)which corresponds to the fluctuations of the height fields being Gaussian with variance σ . Note that since { φ ( m ) n } and { φ ( m ) n (cid:48) } only differ in φ ( m ) j (cid:48) and φ ( m ) j ( φ ( m ) n = φ ( m ) n (cid:48) otherwise), we can Taylor expand (cid:98) Q { φ ( m ) n } , { φ ( m ) n (cid:48) } in λ j asfollows: (cid:98) Q { φ ( m ) n } , { φ ( m ) n (cid:48) } = (cid:32) P ( { φ ( m ) n } ) (cid:12)(cid:12)(cid:12) { φ ( m ) n } (cid:69) − P ( { φ ( m ) n (cid:48) } ) (cid:12)(cid:12)(cid:12) { φ ( m ) n (cid:48) } (cid:69)(cid:33) (cid:32) P ( { φ ( m ) n } ) (cid:68) { φ ( m ) n } (cid:12)(cid:12)(cid:12) − P ( { φ ( m ) n (cid:48) } ) (cid:68) { φ ( m ) n (cid:48) } (cid:12)(cid:12)(cid:12)(cid:33) = λ j δδφ ( m ) j (cid:32) P ( { φ ( m ) n } ) (cid:12)(cid:12)(cid:12) { φ ( m ) n } (cid:69)(cid:33) (cid:32) P ( { φ ( m ) n } ) (cid:68) { φ ( m ) n } (cid:12)(cid:12)(cid:12)(cid:33) (cid:32) δδφ ( m ) j (cid:33) † + O (cid:0) λ j (cid:1) . (D7)1Using Eq. (D6) and (D7), we obtain H ( m ) = σ (cid:88) j (cid:90) (cid:89) n dφ ( m ) n P ( { φ ( m ) n } ) δδφ ( m ) j (cid:32) P ( { φ ( m ) n } ) (cid:12)(cid:12)(cid:12) { φ ( m ) n } (cid:69)(cid:33) (cid:32) P ( { φ ( m ) n } ) (cid:68) { φ ( m ) n } (cid:12)(cid:12)(cid:12)(cid:33) (cid:32) δδφ ( m ) j (cid:33) † + O (cid:0) σ (cid:1) (D8)Using the fact that σ ∼ ∆ x , up to an overall dimensionful factor, the continuum Hamiltonian is H ( m ) = γ (cid:90) dx (cid:90) D φ ( m ) P ( φ ( m ) ( x )) δδφ ( m ) ( x ) (cid:18) P ( φ ( m ) ( x )) (cid:12)(cid:12)(cid:12) φ ( m ) ( x ) (cid:69)(cid:19) (cid:18) P ( φ ( m ) ( x )) (cid:68) φ ( m ) ( x ) (cid:12)(cid:12)(cid:12)(cid:19) (cid:18) δδφ ( m ) ( x ) (cid:19) † . (D9)Next, from Eqs. (C4) and (B5) we compute P ( φ ( m ) ( x )) δδφ ( m ) ( x ) (cid:18) P ( φ ( m ) ( x )) (cid:12)(cid:12)(cid:12) φ ( m ) ( x ) (cid:69)(cid:19) = (cid:18) δδφ ( m ) ( x ) + β δE ( φ ( m ) ( x )) δφ ( m ) ( x ) (cid:19) (cid:12)(cid:12)(cid:12) φ ( m ) ( x ) (cid:69) , (D10)and δE ( φ ( m ) ( x )) δφ ( m ) ( x ) = − δδφ ( m ) ( x ) (cid:82) L d y (cid:0) ∂ m +1 y φ ( m ) ( y ) (cid:1) = − (cid:82) L d y ∂ m +1 y φ ( m ) ( y ) ∂ m +1 y δ ( y − x )= 2 × ( − m ∂ m +1) x φ ( m ) ( x ) , (D11)where we have integrated by parts. Finally, using Eq. (D10) we obtain the functional form of the Hamiltonian H ( m ) H ( m ) = γ (cid:90) dx (cid:18) δδφ ( m ) ( x ) + ( − m κ ∂ m +1) x φ ( m ) ( x ) (cid:19) (cid:18) − δδφ ( m ) ( x ) + ( − m κ ∂ m +1) x φ ( m ) ( x ) (cid:19) = γ (cid:32)(cid:90) dx (cid:32) Π ( m )2 κ ∂ m +1) x φ ( m ) ( x )) (cid:33) + ( − m +1 κ (cid:90) dx ∂ m +1) x δ ( x ) (cid:33) , (D12)matching the expressions in Eqs. (48), (49) and (53) in the main text, up to an overall constant energy shift.
2. Fokker-Planck approach
An alternate, more standard approach for deriving the continuum Hamiltonian proceeds by invoking an analogywith Langevin dynamics [71]. We begin with the observation that that the ground state wave functional of thecontinuum RK Hamiltonian—when interpreted as the equilibrium probability distribution of the Markov processin the continuum—given by the absolute square of the amplitudes in Eqs. (B6) and (D1), resembles a Boltzmanndistribution W ( φ ( m ) ( x )) = 1 Z exp( − κE ( φ ( m ) ( x ))) , (D13)where κ is the “inverse-temperature” and E ( φ ( m ) ( x )) is the energy of the system. Given that the underlying Markovprocess relaxes to the Gibbs distribution, we expect the long-wavelength behavior at late-times to be captured byLangevin dynamics of the generalized height fields φ ( m ) ( x, t ) [71, 92] dφ ( m ) ( x, t ) dt = − γκ δE ( φ ( m ) ( x, t )) δφ ( m ) ( x ) + ζ ( x, t ) . (D14)Here, γ is a constant that sets the overall rate of relaxation and ζ ( x, t ) is δ -correlated “white-noise” i.e., it satisfies aGaussian distribution with (cid:104) ζ ( x, t ) (cid:105) = 0 , (cid:104) ζ ( x, t ) ζ ( x (cid:48) , t (cid:48) ) (cid:105) = γδ ( x − x (cid:48) ) δ ( t − t (cid:48) ) . (D15)To see that the distribution Eq. (D13) is indeed the equilibrium distribution of Eq. (D14), we derive the time-evolution equation for the joint probability densities, also known as the Fokker-Planck equation. We follow methodselucidated in Refs. [71, 74, 92, 93]. Note that although we consider OBC in the main text, we will use PBC in thefollowing for convenience. Taking the the Fourier transform of Eq. (D14), we obtain1 √ L (cid:88) k (cid:54) =0 e ikx dφ ( m ) k dt = − γκ (cid:88) k (cid:54) =0 dE ( φ ( m ) k ) dφ ( m ) − k δφ ( m ) − k δφ ( m ) ( x ) + 1 √ L (cid:88) k (cid:54) =0 e ikx ζ k ( t ) = ⇒ dφ ( m ) k dt = γκk m +1) φ ( m ) k + ζ k ( t ) , (D16)2where we have defined φ ( m ) ( x ) = 1 √ L (cid:88) k (cid:54) =0 e ikx φ ( m ) k , φ ( m ) k = 1 √ L (cid:90) L dx e − ikx φ ( m ) ( x ) , E ( φ ( m ) k ) ≡ (cid:88) k (cid:54) =0 k m +1) φ ( m ) k φ ( m ) − k , (D17)where k is a discrete momentum variable, and we have suppressed the dependence on t and without loss of generalityset φ ( m ) k =0 = 0. Further, using Eqs. (D15) and (D16), we obtain (cid:104) ζ k ( t ) (cid:105) = 0 , (cid:104) ζ k ( t ) ζ k (cid:48) ( t ) (cid:105) = γδ k (cid:48) , − k δ ( t − t (cid:48) ) . (D18)Following standard methods in Langevin theory [92, 93], we write the expression for the probability distribution attime t + ∆ t as W ( φ ( m ) k , t + ∆ t ) = (cid:90) D φ ( m ) (cid:48) P ( φ ( m ) k , t + ∆ t | φ ( m ) k (cid:48) , t ) W ( φ ( m ) k (cid:48) , t ) , D φ ( m ) (cid:48) ≡ (cid:89) k (cid:54) =0 dφ ( m ) k (cid:48) (D19)where W ( φ ( m ) k , t ) represents the joint probability distribution of the height variables { φ ( m ) k } at time t , and P ( φ ( m ) k , t +∆ t | φ ( m ) k (cid:48) , t ) denotes the probability of transition of the height variables from { φ ( m ) k (cid:48) } to { φ ( m ) k } in time ∆ t . Defining∆ φ ( m ) k ≡ φ ( m ) k − φ ( m ) k (cid:48) , (D20)we rewrite Eq. (D19) as W ( φ ( m ) k , t + ∆ t ) = (cid:90) D (cid:16) ∆ φ ( m ) (cid:17) P ( φ ( m ) k , t + ∆ t | φ ( m ) k − ∆ φ ( m ) k , t ) W ( φ ( m ) k − ∆ φ ( m ) k , t ) , D (cid:16) ∆ φ ( m ) (cid:17) ≡ (cid:89) k (cid:54) =0 d (∆ φ ( m ) k ) . (D21)We then Taylor-expand the integrand in Eq. (D21) as P ( φ ( m ) k − ∆ φ ( m ) k + ∆ φ ( m ) k , t + ∆ t | φ ( m ) k − ∆ φ ( m ) k , t ) P ( φ ( m ) k − ∆ φ ( m ) k , t )= (cid:32) − (cid:80) k (cid:54) =0 ∆ φ ( m ) k ddφ ( m ) k + (cid:80) k,k (cid:48) (cid:54) =0 ∆ φ ( m ) k ∆ φ ( m ) k (cid:48) d dφ ( m ) k dφ ( m ) k (cid:48) + · · · (cid:33) P ( φ ( m ) k + ∆ φ ( m ) k , t + ∆ t | φ ( m ) k , t ) P ( φ ( m ) k , t ) . (D22)Defining the quantities (cid:104) ∆ φ ( m ) k (cid:105) ≡ (cid:90) D (∆ φ ( m ) )∆ φ ( m ) k P ( φ ( m ) k + ∆ φ ( m ) k , t + ∆ t | φ ( m ) k , t ) (cid:104) ∆ φ ( m ) k ∆ φ ( m ) k (cid:48) (cid:105) ≡ (cid:90) D (∆ φ ( m ) )∆ φ ( m ) k ∆ φ ( m ) k (cid:48) P ( φ ( m ) k + ∆ φ ( m ) k , t + ∆ t | φ ( m ) k , t ) , (D23)we can write the Kramers-Moyal expansion [93] of Eq. (D21) as W ( φ ( m ) k , t + ∆ t ) = − (cid:88) k (cid:54) =0 ddφ ( m ) k (cid:104) ∆ φ ( m ) k (cid:105) + 12 (cid:88) k,k (cid:48) (cid:54) =0 d dφ ( m ) k dφ ( m ) k (cid:48) (cid:104) ∆ φ ( m ) k ∆ φ ( m ) k (cid:48) (cid:105) + · · · P ( φ ( m ) k , t ) , (D24)where the derivatives also act on P . Using Eq. (D16), we find that (cid:104) ∆ φ ( m ) k (cid:105) = γκk m +1) φ ( m ) k ∆ t + (cid:90) t +∆ tt dt (cid:48) (cid:104) ζ k ( t (cid:48) ) (cid:105) = γk m +1) φ ( m ) k ∆ t , (cid:104) ∆ φ ( m ) k ∆ φ ( m ) k (cid:48) (cid:105) = γ ∆ t δ k (cid:48) , − k + O ((∆ t ) ) , (D25)where we have used Eq. (D18). Using Eqs. (D24) and (D25), we obtain dW ( φ ( m ) k , t ) dt = γ (cid:88) k (cid:54) =0 ddφ ( m ) k (cid:32) ddφ ( m ) − k + 2 κk m +1) φ ( m ) k (cid:33) W ( φ ( m ) k , t ) . (D26)3Note that Eq. (D26) is a master equation of the form of Eq. (15). We can then directly verify that the equilibriumprobability distribution is given by W ( φ ( m ) k ) = 1 Z exp − κ (cid:88) k (cid:54) =0 k m +1) φ ( m ) k φ ( m ) − k , (D27)which is the Fourier transform of the Boltzmann distribution Eq. (D13).To obtain an symmetric (SMF decomposable) quantum Hamiltonian corresponding to the master equation Eq. (D26)with the wavefunction of Eq. (B6) as the ground state, we perform a similarity transformation on the transition matrixin Eq. (D26); equivalently, we write Eq. (D26) in the form [71, 72, 74] d Φ ( m ) ( φ ( m ) k , t ) dt = − H ( m ) Φ ( m ) ( φ ( m ) k , t ) , Φ ( m ) ( φ ( m ) k , t ) = W ( φ ( m ) k , t ) (cid:113) W ( φ ( m ) k ) , (D28)which, using Eqs. (D26) and the form of Eq. (D27), can be shown to be [71] d Φ ( m ) ( φ ( m ) k , t ) dt = − γ (cid:88) k (cid:54) =0 (cid:32) − ddφ ( m ) k + κk m +1) φ ( m ) k (cid:33) (cid:32) ddφ ( m ) − k + κk m +1) φ ( m ) − k (cid:33) Φ ( m ) ( φ ( m ) k , t ) . (D29)Thus, we find that the continuum Hamiltonian H ( m ) is given by H ( m ) = γκ (cid:88) k (cid:54) =0 k m +1) Q † m ( k ) Q m ( k ) , (D30)where the creation and annihilation operators are defined as Q † m ( k ) = (cid:114) κk m +1) (cid:32) − ddφ ( m ) k + κk m +1) φ ( m ) − k (cid:33) , Q m ( k ) = (cid:114) κk m +1) (cid:32) ddφ ( m ) − k + κk m +1) φ ( m ) k (cid:33) (D31)and satisfy the algebra[ Q m ( k ) , Q m ( k (cid:48) )] = 0 , (cid:2) Q † m ( k ) , Q † m ( k (cid:48) ) (cid:3) = 0 , (cid:2) Q m ( k ) , Q † m ( k (cid:48) ) (cid:3) = δ k,k (cid:48) . (D32)As a direct consequence of the construction, the equilibrium wavefunction Φ ( m )0 ( φ ( m ) k ) is given byΦ ( m )0 ( φ ( m ) k ) = 1 √ Z exp − κ (cid:88) k (cid:54) =0 k m +1) φ ( m ) k φ ( m ) − k , (D33)which is precisely the wavefunction of Eq. (B6) written in terms of Fourier variables. Using Eqs. (D30) and (D32), wedirectly obtain the dispersion relation of Eq. (56). Further, the real-space expressions for H ( m ) and the creation andannihilation operators shown in Eqs. (48)-(50) can be directly obtained by rescaling Q m ( k ) and Q † m ( k ) in Eq. (D31)by a factor of √ κk m +1)+1)