Spectral statistics in constrained many-body quantum chaotic systems

We study the spectral statistics of spatially-extended many-body quantum systems with on-site Abelian symmetries or local constraints, focusing primarily on those with conserved dipole and higher moments. In the limit of large local Hilbert space dimension, we find that the spectral form factor $K(t)$ of Floquet random circuits can be mapped exactly to a classical Markov circuit, 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 of the RK-Hamiltonian lower bounds the Thouless time $t_{\mathrm{Th}}$ of the underlying circuit. For systems with conserved higher moments, we derive a field theory for the corresponding RK-Hamiltonian by proposing 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_{\mathrm{Th}}$. In particular, we analytically argue that in a system of length $L$ that conserves the $m^{th}$ multipole moment, $t_{\mathrm{Th}}$ scales subdiffusively as $L^{2(m+1)}$. We also show that our formalism directly generalizes to higher dimensional circuits, and that in systems that conserve any component of the $m^{th}$ multipole moment, $t_{\mathrm{Th}}$ has the same scaling with the linear size of the system. Our work therefore provides a general approach for studying spectral statistics in constrained many-body chaotic systems.

Despite the significant difficulty in analytically studying dynamics in generic many-body quantum systems, substantial progress has been made in delineating the dynamics of chaos in random quantum circuits via the two-point spectral form factor K(t), defined in terms of the spectral properties of the evolution operator W as where {θ m } is the set of eigenphases of W , W (t) ≡ W t denotes the t th power of W , N is the Hilbert space dimension, and · denotes the average over an ensemble of statistically similar systems. The SFF is the Fourier transform of the two-level correlation function, with time t as the variable conjugate to ω ∼ θ m −θ n , the (quasi)-energy difference. For W with Poisson-distributed eigen-levels, K(t) = N for all t, while for W chosen as random matrices 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]. The SFF serves as a barometer for quantum chaos, as the appearance of RMT predictions in the SFF provide a crisp time-scale for characterizing the many-body chaos in a finite system, which we take as defining the Thouless The black solid line is the generic behavior for many-body quantum chaotic systems, characterized by the Thouless time t Th and the Heisenberg time tHei, the latter of which scales exponentially with system size L. The early time behaviour of K(t) depends on the details of the system, particularly on whether the system is defined by a Floquet operator or a Hamiltonian. time t Th 1 ; specifically, we define t Th as the time scale at which the SFF approaches the CUE RMT behavior (see Fig. 1). In contrast to nearest-neighbour level spacing distributions, the SFF encodes spectral correlations at all time (equivalently, energy) scales and is also relatively simple to analyse, as it involves only two sums over the (quasi)-energy eigenvalues.
One approach for analytically computing the SFF exploits a self-duality present in certain models [44,45]. However, the applicability of this approach is limited only to self-dual circuits and does not extend to generic interacting circuits, possibly with conserved quantities. A second, more generic approach [46][47][48][49][50] studies Floquet random quantum circuits (FRQCs) in the limit of large local Hilbert space dimension, which are amenable to exact analytic calculations of the SFF and hence enable one to study the implications of conserved quantities on dynamics. For circuits with a globally conserved U(1) charge, t Th was shown to scale diffusively as ∼ L 2 , validating the idea that RMT behaviour is established only once all conserved quantities have diffused through the system. In contrast, for certain systems without any conserved charge, t Th ∼ log L [47]. In these latter systems, the SFF is not sensitive to the slower "ballistic" dynamics due to locality and causality which means the operator spreading differs from RMT dynamics up to longer times of order L.
In this paper, we develop a general approach for studying features of the spectral statistics in constrained quantum chaotic many-body systems, focusing on those with conserved dipole and higher moments. Using the SFF K(t) as a diagnostic for many-body chaos, we extract the scaling 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 auxiliary 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, with nodes representing an appropriate set of basis states and links representing possible moves between them. c) Generalized height field representation for each spin configuration, with symmetries enforced via boundary constraints. The RK ground state is then written as an equal-weight superposition of height fields (equivalently, as random walks in the space direction).
spins through which the constraints are imposed. In the large-q (q → ∞) limit, we express K(t) in terms of a classical Markov circuit which inherits the constraints of the underlying FRQC, as illustrated in Fig. 2(b). Utilizing an established correspondence between classical Markov processes and Rokhsar-Kivelson (RK) type Hamiltonians, we can equivalently relate K(t) at late times to the partition function of a positive-definite, frustration-free RK-Hamiltonian 2 acting solely on the spin DOFs. Consequently, a lower bound on t Th can be extracted from the spectral gap of the RK-Hamiltonian; this mapping hence allows us to borrow techniques from equilibrium physics to establish dynamical properties of the underlying circuit [71][72][73].
For circuits with conserved higher moments in one dimension, we find a 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 continuum parent Hamiltonian for the corresponding GS and establish a lower bound on t Th . We find a sub-diffusive scaling of t Th ∼ L 2(m+1) for circuits of length L with a conserved m th moment i.e., the timescale at which random matrix behavior ensues in such systems is parametrically longer than in systems with only a conserved U(1) charge, which spreads diffusively (m = 0: t Th ∼ L 2 ). We also find similar results in higher dimensions by constructing continuum representations of the ground state and the RK-Hamiltonian in terms of gener-2 As discussed later in Sec. III, we take "RK-Hamiltonian" to mean a quantum Hamiltonian that is proportional to the transition matrix of a Markov process that satisfies detailed balance. This taxonomy stems from the quantum dimer context [70], where the terms "RK-Hamiltonian" and "quantum dimer model at the RK-point" are often used interchangeably.
alized tensor fields, which predicts t Th ∼ L 2(m+1) for a system with linear-size L with conserved m th moments in all directions. Note that while the field theories are specific to higher-moment conserving systems, the mapping from K(t) to an emergent RK-Hamiltonian in the large-q limit holds generally. This paper is organized as follows. We start by defining the class of FRQCs under consideration in Sec. II. In Sec. III, we study these circuits in the limit of large local Hilbert space and establish a mapping between the SFF of the FRQC and the dynamics of a classical Markov chain, which we further show is equivalent at late times to the 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 the spectral gap of this emergent Hamiltonian. Focusing on systems with conserved higher moments in Sec. IV, we verify that t Th ∼ L 2 in charge conserving systems and provide numerical evidence for subdiffusive scaling of t Th in systems that additionally conserve the global dipole moment. In Sec. V, we take the continuum limit of the emergent RK-Hamiltonian for systems which conserve all moments up to the m th highest moment. We extract a bound on t Th from the dispersion relation of the continuum Hamiltonian and show that t Th ∼ L 2(m+1) . In Sec. VI, we generalise our results to multipole conserving systems in higher dimensions. We conclude in Sec. VII with a discussion of open questions and future directions.

II. CONSTRAINED FLOQUET RANDOM QUANTUM CIRCUITS
Our primary object of interest in this paper is the Thouless time t Th of constrained many-body quantum chaotic systems, where we define t Th as the time-scale after which the behaviour of the SFF K(t) closely approaches RMT predictions. To probe K(t), we consider one-dimensional L-site spatially-random FRQCs with local Hilbert space at each site of the chain given by H loc = C q ⊗ C 2s+1 , where C q and C 2s+1 are the local Hilbert spaces of the color and spin DOFs respectively. The color DOFs facilitate Haar averaging and allow us to retain analytical control in the q → ∞ limit [36,46,47]; the spins, on the other hand, allow us to encode on-site Abelian symmetries, such as U(1) charge conservation (previously considered in Refs. [36,48]), or impose local constraints on the dynamics.
More precisely, we consider Floquet circuits defined by a time-evolution operator W over a single period composed of unitary gates acting on a finite number ≥ min of contiguous sites, where min is the "minimal" gate size for non-trivial local dynamics under the symmetry or constraints of interest. Without loss of generality, local dynamics on all sets of contiguous sites within a single time-period can be ensured by choosing W to be composed of layers of operators { W a }, where W a is composed of r = L/ spatially random local unitary gates { U [j,j+ −1] }, and has the form: 3 where a is the layer index. As shown in Fig. 2(a), U [a+(n−1) ,a+n −1] labels the n th local gate in the a th layer and acts non-trivially only on sites j ∈ {a + n − , . . . , a + n − 1}, where 1 ≤ n ≤ r and the site index j is defined mod L for periodic boundary conditions (PBC). Each of the local gates has the following block-diagonal structure: where α denotes each set (block) of -site spinconfigurations within this gate which are connected through local moves permissible under symmetries or constraints, and D denotes the total number of such blocks. Block α contains d α spin configurations within this gate, with D α=1 d α = (2s + 1) . Note that we do not impose any constraints on the color DOFs, only on the spins. Each u(j, α) is thus a d α q ×d α q unitary drawn independently from the Haar ensemble acting on the states in block α, while it gives zero when it acts on all other states. In particular, for systems without any symmetries or dynamical constraints, the local gates U [j,j+ −1] are (2s + 1) q × (2s + 1) q independent Haar random unitaries. The block diagonal structure of the local gates encodes the symmetries or constraints of interest. Specifically for systems which have global symmetry sectors labelled by a set of quantum numbers S = {s 1 , s 2 , . . . }, each local -site gate U is block diagonal, with each block containing all spin states with the same S. As a technical aside, we note that since we are keeping the gate size fixed while treating all transitions involving those sites on equal footing, this also includes all allowed processes involving spin transitions on any subset of those sites.
As an example, let us consider an FRQC with s = 1/2 DOFs which preserves the total charge Q 0 = L x=1 S z x of the spins, where S z x is the Pauli-Z matrix acting on site x [48]. To allow non-trivial dynamics, we choose to be equal to min = 2, such that each local gate U is a 4q 2 × 4q 2 block-diagonal matrix. Each local gate is composed of D = 3 blocks: two q 2 × q 2 blocks act on the tensor product subspaces associated with the spin configurations |↑↑ and |↓↓ , and a single 2q 2 ×2q 2 block acts on the subspace associated with the spin configurations |↑↓ , |↓↑ . Each of these three blocks locally preserves the U(1) charge over 2-sites and is an independently-drawn Haar random unitary.
In this paper, we mainly focus on circuits which conserve not only the total charge, but also all higher moments up to the m th moment . Such circuits neatly fall into the larger class of FRQCs defined earlier via Eqs. (2) and (3). For instance, for systems with both charge and dipole moment conservation, we can consider s = 1 DOFs and = min = 3 site gates, where each local gate is a 27q 3 × 27q 3 block diagonal matrix, with each block corresponding to those spin configurations which are connected under local (3-site) dipole moment preserving dynamics (see Ref. [60] for details). It is straightforward to generalize the above circuits to higher spins s, larger gate sizes , and higher moment conservation laws. While our primary focus in this paper will be systems with higher conserved moments, the FRQCs defined in Eqs. (2) and (3) define a much broader class, including those with arbitrary onsite Abelian symmetries as well as circuits which obey dynamical constraints, such as those present in the PXP model [56].
We characterise the spectral features of the above class of FRQCs using the SFF defined in Eq. (1). For a circuit W invariant under a set of global symmetries, corresponding to a set of operators { S 1 , S 2 , . . . }, we have [ W , S i ] = 0 ∀ i. Therefore, is block-diagonal and quasienergy levels of W from blocks W (S) , corresponding to distinct quantum number sectors S = {s 1 , s 2 , . . . }, do not repel [21]. In addition, certain systems, such as those with higher moment symme- tries, further exhibit the phenomenon of Hilbert space fragmentation, wherein the dynamics does not connect all products states in the S z (equivalently, charge) basis even within the same symmetry sector [13,14,61,62]. As a consequence, within each symmetry sector S there may exist up to exponentially many disjoint "Krylov subspaces," labelled by K i.e., where D (S) denotes the number of disjoint Krylov subspaces generated from product states with the same quantum numbers S. This fragmented structure of W is schematically depicted in Fig. 3. Hence, given the possibility of Hilbert space fragmentation in quantum manybody systems, we define the SFF restricted to a given Krylov subspace K: where the subscript K denotes the restriction of W to a Krylov subspace and · denotes averaging over the Haar random unitaries in the FRQC. Note that this definition encompasses systems with global symmetries but no fragmentation, 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 MARKOV CHAIN AND EMERGENT RK-HAMILTONIAN
Computing the SFF Eq. (6) for many-body quantum systems is analytically difficult in general. Refs. [46][47][48] developed a diagrammatic approach for evaluating the ensemble averaging in (6), by effectively "integrating out" the color DOFs for random quantum circuits with charge conservation. In Appendix A, we generalize this technique to the general class of constrained FRQCs discussed in the previous section and find that, to leading order in the large-q limit, where the factor of |t| stems from |t| leading order diagrams as q → ∞ [46]. The Markov matrix M is a classical bi-stochastic 4 circuit acting only on effective spin-s DOFs and is composed of local -site gates. Furthermore, M inherits the circuit geometry, symmetries, and Krylov subspaces of W and can be expressed as where, in analogy with Eq. (2), a is the layer index and m [a+(n−1) ,a+n −1] labels the n th local gate in the a th layer, acting on sites j ∈ {a + n − , . . . , a + n − 1}. As before, 1 ≤ n ≤ r and j is defined mod L for PBC. However, unlike the underlying gates U , the -site gates { m [j,j+ −1] } are non-random: where the d α 's are the sizes of the blocks of -site spin configurations that are dynamically connected and D is the number of dynamically connected blocks in m. The fact that Eq. (9) retains the block-diagonal form with equal matrix elements (within each block) is consistent with the fact that the dynamical constraints are imposed via the spin DOFs and that the local Haar random gates are invariant upon a change of basis.
Since M inherits the symmetries and Krylov subspaces (Eqs. (4) and (5)) of W , and M also has a block diagonal structure (see Fig. 3), where each block is itself an irreducible bi-stochastic matrix; hence, each block has a unique largest magnitude eigenvalue 1. Restricting our attention to the subspace K of interest, let us denote the eigenvalues of the corresponding block by {Λ 5 We can then write 5 Note that the eigenvalues of M can be complex since is not a symmetric matrix due to the brick-wall structure of the circuit At sufficiently long times 6 t 1, we can then expand Eq. (10) as is the magnitude of the second largest eigenvalue of M restricted to the subspace K and d (K) is a constant factor encoding its degeneracy along with any complex phases.
Using Eq. (11), we see that the SFF K ∞ (t; K) approaches the linear in |t| RMT behavior after a time Thus, we have shown that extracting the scaling behavior of t Th with system size L is equivalent, in the largeq limit, to the problem of obtaining the scaling of the "gap" 7 ∆ (K) M (L) of the Markov circuit M (within the subspace K) with L. Henceforth, we will restrict our attention to a single (exponentially large) subspace K of M and suppress the sub/superscript K for ease of notation.
To determine the gap ∆ M , we will now establish a relation between M (within a subspace K) and a quantum Hamiltonian. We proceed by first observing that M corresponds to the classical stochastic time evolution of a probability density p(t), defined over all product states in the usual Z-basis for the spin Hilbert space: where M αβ represents the matrix elements of M and the bistochasticity of M ensures that the total probability is conserved under time-evolution. In particular, under the action of each local gate m [j,j+ −1] (see Eq. (9)), the probability density p evolves (with equal probability) to all product states that can be reached via the allowed local moves i.e., moves that are allowed by the constraints imposed on the spin DOFs. Now consider starting with a probability density where all the weight is concentrated on a single product state within some subspace K. Under the stochastic evolution, the probability density will eventually reach a unique equilibrium state, specified by the uniform distribution over all product states in K i.e., by the eigenvector of M corresponding to the eigenvalue Λ 1 = 1. Thus, obtaining ∆ M is related to obtaining the 6 Technically, we require that t H t 1, but the so-called Heisenberg time t H proportional to the inverse level spacing is infinite in the large-q limit. 7 The gap of a bistochastic matrix is traditionally defined to be , which we have approximated as ∆ inverse of the mixing time for this process, which is in general not analytically tractable.
To derive the gap ∆ M , we are interested in the stochastic evolution under M at time-scales of O(1/∆ M ). If ∆ M → 0 in the thermodynamic limit (i.e., if M is gapless), we expect that the dynamics under M at late times is well approximated by a continuous-time process. 8 In particular, since M corresponds to a stochastic process with local moves occurring independently with equal probability, its late-time behavior should be well approximated by that of a continuous time process composed of the same local moves occurring at equal rates, with the additional requirement of detailed balance. In other words, the evolution of the probability density p should be governed by a Master equation of the form: where p α (t) is the probability of a classical system occupying state α and T αβ is the transition rate from state β to state α, which we specify below. Defining T βα , we can rewrite Eq. (14) as a matrix equation in terms of the transition matrix T , which ensures that local moves that occur with equal probability in the discrete-time process Eq. (13) occur at equal rates in the continuous-time process Eq. (16). The late time behavior of the stochastic process governed by M is then given by Eq. (15) with where Γ is an overall positive constant 9 that sets the rate at which local moves occur and H is defined as The matrix Π [j,j+ −1] is a projector and has the form 8 Note that this approximation does not hold for gapped systems as they relax to their equilibrium distribution on time-scales of O(1), that are much smaller than the time-scale at which the continuous time description is valid. 9 As we show in Sec. IV, Γ is a non-universal constant determined by the detailed microscopic properties of the underlying FRQC e.g., it depends on the number of layers . However, obtaining its precise value is not important for our purposes.
thereby ensuring that the transitions taking place in the continuous-time process are identical to those specified by the gates m [j,j+ −1] . In fact, H can be interpreted as a quantum Hamiltonian in the product state basis (in the Z-basis) for the spin Hilbert space and has the same symmetries and Krylov subspaces as the stochastic circuit M . More importantly, H belongs to the class of so-called RK-Hamiltonians, where an RK-Hamiltonian is defined as a quantum Hamiltonian that is proportional to the transition matrix T of a discrete classical stochastic process which satisfies detailed balance [72]. Consequently, the ground state wave function of an RK-Hamiltonian can be interpreted as a classical equilibrium distribution, its low-lying excited states correspond to classical relaxation modes, and its gap coincides with the relaxation time of the corresponding transition matrix. Such Hamiltonians were first studied in the context of quantum dimer models [70] and have subsequently been explored extensively in various settings [71,72,74].
To emphasize the relation between M and an emergent 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-time Schrödinger evolution under H RK Eq. (17). In effect, the correspondence between M and H RK amounts to a relation between Tr K [ M t ] and the partition function of H RK restricted to K at an inverse temperature β = Γt, namely: . We can therefore approximate the SFF Eq. (10) at late times as such that the gap of M is related to ∆ RK (L), the gap of the Hamiltonian H RK (restricted to the subspace K): Indeed, as we discuss in Sec. IV, we find numerical evidence that supports Eq. (20) in multipole conserving circuits. Further evidence for the correspondence between H RK and M is obtained by studying the system-size dependence of the overlap between the "first-excited eigenstates" of |ψ RK and ψ M of H RK and M respectively. As shown in the inset of Fig. 4(a), we find that this overlap approaches 1, suggesting that |ψ RK is an asymptotically exact eigenstate of M in the thermodynamic limit.
We also numerically observe that the overlap does not approach 1 in the cases when M is gapped, further suggesting that correspondence between M and H RK is only valid when M is gapless.
The preceding discussion shows that obtaining the gap of H RK is sufficient for obtaining the scaling of the Thouless time t Th . Taken together, Eqs. (12) and (20) constitute one of the central results of this paper, whereby a dynamical property of the FRQC is determined by the low-energy, equilibrium behavior of an emergent quantum Hamiltonian. Since we made no reference to the microscopic structure of the underlying circuit, this relationship between t Th and ∆ RK holds generally for the class of circuits specified in Sec. II.
We now briefly discuss some properties of the Hamiltonians H RK . Denoting the basis set of π(d α ) in Eq. (18) by B α , we obtain where C and C represent -site configurations chosen from the basis set B α . We can then re-express H RK as where C, C denotes product states C and C that are connected under the action of H RK and the weights w C,C ≥ 0 are defined in accordance with Eq. (21). From this expression, it is clear that H RK is a positive semidefinite Hamiltonian and the zero-energy ground state wave function within a subspace K is given by where C runs over all the product states in the subspace K and Z is a normalization factor. We remark that, up to an overall normalization factor, |G RK can be interpreted as the equilibrium probability distribution of the stochastic process described by Eqs. (13) or (16). Before proceeding to focus on multipole conserving circuits, we briefly discuss some important aspects of the RK-Hamiltonian H RK which illustrate the potential benefits of the mapping developed in this section. First, H RK is a frustration-free Hamiltonian, i.e., the ground state |G RK is the ground state of each of the terms Π [j,j+ −1] as can be seen using Eqs. (22) and (23), since Q C,C G (K) RK = 0 for any C and C . This fact enables the use of well-known methods for bounding the spectral gap of frustration-free Hamiltonians [75][76][77], which in turn allow us to place a constraint on the scaling exponent of t Th in FRQC with constraints defined in Eq. (2) in the large q limit: We note that for circuits without any conserved quantities, t Th has been shown to scale as log L for certain Floquet models [32,47]. However, as evidenced by Eq. (24), for circuits of the form Eq. (2) (including the circuit discussed in Ref. [46]), such scaling of t Th is suppressed in the q → ∞ limit and we instead find that t Th can scale as an O(1) number in FRQCs in this limit.
Secondly, by virtue of the connection to classical Master equations, shown in Eq. (16), H RK is an example of a stoquastic Hamiltonian, which can be efficiently studied using Quantum Monte Carlo techniques [72,78]. We thus expect that the same techniques can be exploited to efficiently study the late-time features of the SFF in a variety of settings at large-q. Moreover, in the context of spectral graph theory, any Hamiltonian of the form Eq. (22) restricted to a subspace K exactly corresponds to the Laplacian [79] of an undirected graph G, formed by the set of vertices {C} within K and by edges with weights w C,C between the vertices C and C . The gap ∆ RK then corresponds to the gap of the Laplacian of the graph G, which is closely related to the connectivity of G. In particular, the existence of bottlenecks in G, as detected by the Cheeger constant, results in a smaller gap of the Laplacian (and hence, in a larger t Th ). This establishes a clear connection between the nature of transport in the presence of constraints and the connectivity of the Hilbert space under those constraints.
Finally, we note that earlier work has also discussed a relation between the Thouless time t Th of a chargeconserving FRQC and and the spectral gap of a U(1) invariant classical bistochastic circuit: Ref. [48] constitutes a particular case of the results obtained in this paper, being derived in the large-q limit, while Ref. [80] invokes the random phase approximation in a long-range interacting model at finite-q. While both of these works were restricted to specific realizations of U(1) invariant systems, the relations between K ∞ , M , and H RK obtained in this section apply far more generally to the large class of circuits with arbitrary symmetries or constraints discussed in Sec. II.

IV. EXAMPLES FROM MULTIPOLE CONSERVING CIRCUITS
In this section, we move our attention to FRQCs with conserved higher moments and provide explicit examples of the mapping established in Sec. III. Specifically, we consider circuits which conserve all moments of charge up to the m th highest moment, where the m th multipole moment is defined as where m = 0 (1) corresponds to the charge (dipole) conserving case. Where necessary, we will use the labels {Q m } to denote the m th multipole moment quantum number. 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, the stochastic circuit M for a spin-1/2 U(1) charge conserv-ing circuit with gate size = 2 (for PBC) is given by where the local 2-site gate m [j,j+1] is written in the (ordered) basis {|↑↑ , |↑↓ , |↓↑ , |↓↓ }. Using Eqs. (18) and (26), we find that Π [j,j+1] maps onto a ferromagnetic spin-1/2 Heisenberg term where Π [j,j+1] is written in the (ordered) basis {|↑↑ , |↑↓ , |↓↑ , |↓↓ } and σ = (σ x , σ y , σ z ), with σ i the usual Pauli matrices. For the charge conserving FRQC (with = 2), we hence find that H RK is the Bethe-Ansatz integrable ferromagnetic Heisenberg model, whose integrability was exploited in Ref. [48] to study the late-time behavior of the SFF for this circuit. According to Eq. (23), the unique ground state within any charge sector Q 0 is the equal amplitude superposition of all product states within that symmetry sector. Indeed, such a state belongs to the SU(2) multiplet of the spin-polarized ferromagnetic state with total spin Q 0 = L. Moreover, the low-energy excitations above the ferromagnetic state in the Heisenberg model of Eq. (27) are exactly known to be spin waves with dispersion (k) = 2 sin 2 (k/2). As a consequence of the SU(2) symmetry of the Heisenberg model, the lowest energy excited state within each symmetry sector belongs to the multiplet of spin-wave states with total spin Q 0 = L − 1; the gap of H RK in any Q 0 = L sector is then given by which is the energy corresponding to the lowest non-zero momentum spin-wave. Using Eqs. (12) and (20), we find that the Thouless time in any quantum number sector in an FRQC with U(1) charge conservation scales diffusively with system size i.e., t Th ∼ L 2 . For charge conserving systems with higher spins or larger gate sizes , H RK is no longer integrable in general, but, as shown in Fig. 4(b), we numerically observe the same diffusive scaling ∆ M (L) ∼ ∆ RK ∼ L −2 for the systems we studied.
In fact, we find that the gaps are identical for spin-1/2 and spin-1 systems with gate size = 2 even though the  20)). The inset shows the overlap of normalized "first-excited eigenstates" of HRK (|ψRK ) and M ( ψ M ). Note that the overlap approaches 1 with increasing system size, suggesting that the |ψRK is an asymptotically exact first-excited state of M . (b) Scaling of the gaps of HRK with system size for systems with charge (green, blue) and dipole moment (orange, red) conservation. Note that the gaps 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/Krylov subspace containing the state |0 0 · · · 0 0 for OBC. rest of the spectrum is different, strongly suggesting a universal origin of the scaling.
As evidenced through the above example, the correspondence between the FRQC and the RK-Hamiltonian unveils a curious feature of the large-q limit. While the original FRQC only has U(1) symmetry, after Haar averaging and taking q → ∞, K ∞ (t; K)-related to H RK through Eq. (19)-exhibits an enlarged SU(2) invariance in the spin DOFs. Indeed, we expect this enlarged symmetry to be a generic feature in the large-q limit, since RK-points typically exhibit enhanced symmetries, although not necessarily SU(2) [81][82][83]. On the other hand, to our knowledge, the emergent integrability in the above example is not generic and is specific to the spin-1/2 system with 2-site gates.
We now turn our attention to systems which conserve the dipole moment Q 1 in addition to the charge Q 0 , for which the nature of low-energy excitations above the ground state of the corresponding RK-Hamiltonian H RK is not immediately apparent. An additional feature in such systems is the fragmentation of the Hilbert space of the FRQC [13,14], which leads to the formation of exponentially many Krylov subspaces (see Eq. (5)). Hilbert space fragmentation is typically classified into two types: strong or weak, where the size of the largest Krylov subspace is respectively a zero or non-zero fraction of the total Hilbert space dimension within a given quantum number sector in the thermodynamic limit. Refs. [13,14] numerically observed that spin-1 and spin-1/2 dipole conserving systems with the minimal gate sizes = min = 3 and = min = 4 respectively show strong fragmentation whereas the inclusion of moves requiring larger gate sizes leads to weak fragmentation. Furthermore, Ref. [63] found that the nature of fragmentation can vary even for a given gate size depending on the quantum number sector. Generically, however, experimentally relevant multipole conserving systems are expected to show weak fragmentation.
In strongly fragmented systems, the ratio between the dimension of the largest Krylov subspace within a symmetry sector and the size of that symmetry sector exponentially decays to zero in the thermodynamic limit. As a consequence, typical initial states do not thermalize [13,14], although certain initial states do thermalize with respect to smaller Krylov subspaces [16,62]. In contrast, for weakly fragmented systems there always exists a dominant Krylov subspace K within a given quantum number sector, such that its size asymptotically approaches that of the symmetry sector in the thermodynamic limit. Due to this, typical eigenstates within a quantum number sector carry non-zero weight in the dominant Krylov subspace of that symmetry sector and look thermal. As a consequence, frozen configurations, despite being exponential in number, are expected to have a negligible effect on t Th in a weakly fragmented system. Since our interest in this work is the behavior of generic multipole conserving systems, we focus only on thermalizing weakly fragmented systems here. Hence, we will study the SFF, the scaling of the Thouless time t Th and the gaps ∆ M (L) and ∆ RK (L) all restricted to the dominant Krylov subspace K within a specified quantum number sector.
In Fig. 4, we show the scaling of the gaps ∆ M (L) and ∆ RK (L) for spin-1/2 and spin-1 dipole conserving systems for several gate sizes > min . Fig. 4(a) shows that the numerics are in good agreement with Eq. (20) i.e., they support the correspondence between the stochastic circuit M and the emergent RK-Hamiltonian H RK developed in Sec. III. In principle, we can also extract the microscopic constant Γ for a specific circuit by comparing the gaps for the corresponding M and H RK . Furthermore, as evident from the numerics shown in Fig. 4(b), we find that ∆ RK (L) scales as ∼ L −β with β > 2; thus, the Thouless time scales subdiffusively t Th ∼ L β (β > 2) for system sizes accessible to exact diagonalization. Importantly, this subdiffusive scaling appears to be a generic feature of weakly fragmented, dipole conserving systems and does not show a strong dependence on the microscopic details ( or s) of the circuit. This mirrors the behavior of systems with only charge conservation ( Fig. 4(a)), which show diffusive scaling of t Th independent of microscopic details.
Due to the longer gate sizes required for systems conserving even higher multipole moments Q m≥2 (e.g., 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 the scaling of the gaps of charge and dipole conserving RK-Hamiltonians H RK i.e., ∆ RK (L) ∼ L −2 for charge conserving systems and ∆ RK (L) ∼ L −4 for dipole conserving systems, regardless of the ultraviolet details of the respective Hamiltonians. The appearance of this universality suggests the existence of a universal field theoretic description which effectively captures the low-energy behavior of generic multipole conserving RK-Hamiltonians, such as the scaling of their gap. The derivation of these universal effective field theories will be the subject of the next section.

V. CONTINUUM LIMIT FOR MULTIPOLE CONSERVING SYSTEMS
As discussed in the previous section, the ground state of an RK-Hamiltonian is well-known as being the equalweight superposition of all states in the corresponding Hilbert space. However, understanding the scaling of the gap ∆ RK requires knowledge of low-lying states above the GS, which are generically not known exactly. Nevertheless, motivated by our numerical observation of a universal scaling of ∆ RK with L for generic charge and (weakly fragmented) dipole conserving RK-Hamiltonians, we derive continuum field theories for mul- N +1 are fixed to be zero.
tipole conserving systems through a coarse-graining procedure, detailed in Appendix D. We find that the resultant continuum field theories accurately capture the ground state and low-energy excitations of the corresponding RK-Hamiltonians, therefore providing an analytic route to understanding the scaling of t Th in the underlying FRQC. Throughout this section, we will only consider OBC. We denote the number of spins as N and the system size as L = N ∆x, where ∆x is the lattice spacing. The continuum limit then corresponds to taking the limits N → ∞ and ∆x → 0 simultaneously, while keeping L fixed. For systems which conserve all moments of charge up to the m th moment (or, m th moment conserving systems), we focus our attention on the quantum number sector S = {Q 0 = 0, Q 1 = 0, · · · , Q m = 0}. As discussed in Sec. IV, for weakly fragmented systems there exists a dominant Krylov subspace within each symmetry sector, such that the size of that subspace asymptotically approaches the size of the full symmetry sector as the gate-size increases. Since taking the continuum limit involves coarse-graining and thus effectively taking the gate-size → ∞, we can neglect the effect of fragmentation in systems with dipole and higher moment conservation, and expect that our analysis holds as long as the sectors we study do not exhibit strong fragmentation.

A. Generalized height fields
In order to take the continuum limit, we first need to introduce "generalized" height fields, in analogy with familiar height fields in the quantum dimer context [83]. When taking the continuum limit of the ground state wave function and H RK , a crucial issue is the restriction to the symmetry sector S; this restriction imposes a global constraint on the spin DOFs {s n }, thereby resulting in a non-local action for system. It is to circumvent precisely this issue, while preserving locality, that we work in terms of generalized height variables {φ (m) } for systems with m th moment conservation. In terms of these variables, the conservation of higher moments {Q m } are expressed as local boundary constraints on the height fields and their derivatives, as opposed to a global constraint on the spin DOFs.
We first illustrate this construction in terms of height variables for systems with charge conservation. The height DOFs {φ (0) n+ 1 2 } are defined on the links of the onedimensional chain as which immediately suggests that the total charge Q 0 (see Eq. (25)) is given by the flux of the height variable through the system: Since Eqs. (29) and (30)  = 0 without loss of generality. Thus, restricting to a given charge sector corresponds to imposing constraints on the boundary height variables once we forgo the spin language for the height representation. Furthermore, as a consequence of Eq. (29), any charge conserving process involving spins {s n , · · · , s n+ −1 } involves ( − 1) height variables {φ }. So, the mapping from the spin DOFs to the height variables preserves the locality of the Hamiltonian.
We now generalize the height representation to systems with m th moment conservation. For a system with all moments up to the m th moment conserved, we recursively define the m th "generalized" height variable φ (m) , that lives on links (sites) if m is even (odd), as with {φ (0) n+ 1 2 } defined in Eq. (29). An example of a charge configuration in a dipole conserving (m = 1) system, expressed in terms of height variables, is depicted in Fig. 5.
As mentioned earlier, the key advantage of forgoing the spin language is that the total m th multipole moments can be expressed in terms of boundary height variables and their "derivatives", an observation that will prove essential when taking the continuum limit. For instance, the total dipole moment can be expressed as where we have used Eqs. (30) and (31). Similarly to the charge conserving case, locality is also preserved when going to the generalized height representation. This can be seen from Eq. (31), since any m th multipole conserving process involving spins {s n , · · · , s n+ −1 } involves ( − m + 1) generalized height variables {φ (m) n+ 1 2 }. We now take the continuum limit by defining generalized height fields through an appropriate rescaling of the height variables by the lattice spacing ∆x i.e., where x = n∆x. Here, ρ(x) can be interpreted as the charge density and, as we will see, the height fields φ (m) (x) are related to the multipolar densities. In terms of the height fields, Eq. (29) in the continuum becomes Note that Eq. (34) closely resembles a Gauss' Law. Similarly to Eq. (30), the total charge Q 0 in the continuum is expressed in terms of the height fields as The preceding discussion illustrates how the global charge constraint on ρ(x) is re-expressed as a local boundary constraint on φ (0) (x), clarifying why the latter is physically more appropriate as the field variable for charge conserving systems. Similarly, using Eq. (33), we can express Eq. (31) in terms of the generalized height fields as In the continuum, the conservation of the m th and all lower moments then amounts to fixing the left and right boundary constraints on φ (m) (x) and its derivatives ∂ n x φ (m) (x) for all n ≤ m. For instance, the total dipole moment and charge can be expressed in terms of boundary constraints on φ (1) (x) and ∂ x φ (1) (x) as where we have invoked Eqs. (34) and (36). In fact, it is straightforward to show that the n th multipole moment can be expressed in terms of the m th height field as (m ≥ n) The multipole moments Q n in Eq. (38) are invariant under a polynomial shift of the height fields φ (m) , under which Eq. (36) is still satisfied: where P (m) (x) is an arbitrary polynomial of degree ≤ m; Eq. (39) can thus be used to set ∂ n x φ (m) (0) = 0 for all n ≤ m without loss of generality, so that In the sector of primary interest, Q n = 0 for all n ≤ m, these boundary constraints further simplify to For OBC, we need to further supplement these boundary constraints, which fix the symmetry sectors, with physical boundary conditions, which ensure that no multipole currents flow through the boundaries. As discussed in Ref. [68], the fundamental hydrodynamic quantities for multipole conserving systems are the charge density ρ(x) and the multipole current J (m) , from which one can infer the conventional charge current; however, it is the multipole current that is fundamental and is related to the charge density as for systems which conserve all moments up to the m th highest moment, which is the generalization of Fick's law to multipole conserving systems [64]. The physical requirement that no multipole current flows through the boundaries, phrased in terms of the height fields, can be stated as B. Ground state Before we obtain the continuum limit of the Hamiltonian H RK , we express the ground state wave function of an m th moment conserving system, discussed in Sec. IV, in terms of the height field φ (m) (x). Recall that the GS Eq. (23) is the equal-weight superposition of all allowed basis states i.e., all possible height variable configurations that satisfy the boundary constraints, which fix the quantum number sectors.
Treating the spin-s DOFs as "random variables" that assume integer or half-integer values in [−s, s], under coarse-graining the distribution of the spins flows to a Gaussian as a direct consequence of the central limit theorem [84]. The variance of the resulting coarse-grained DOFs then scales as σ 2 = ∆x/κ where ∆x is the lattice spacing and κ is a parameter chosen such that microscopic correlation functions are accurately reproduced at long distances; effectively, one can think of κ as the coarse-graining length scale. After coarse-graining, the wave functional Φ (m) 0 [ρ(x)] corresponding to a charge density profile ρ(x) is thus simply given by a Gaussian, albeit subject to global constraints specified by the conserved quantities (see Appendix B for details): where Z is a normalization constant and G[ρ(x)] enforces the global symmetry constraints; namely, it fixes the quantum number sector of interest. For instance, for a dipole conserving system in the {Q 0 , Q 1 } sector, To circumvent the global constraint in Eq. (44), it is convenient to work in terms of the m th height fields for m th multipole conserving systems-as discussed in Sec. V A, in this language, the quantum number sectors are instead expressed as local boundary constraints. More explicitly, Eq. (36) allows us to express the wave functional Eq. (44) in terms of the height field φ (m) (x) as which follows from Eq. (41). Note that we also need to impose the physical boundary conditions Eq. (43) on the generalized height fields.
Recall that in the discrete setting, the GS is an equal weight superposition of allowed configurations, while taking the continuum limit introduces Gaussian weights into the GS due to coarse-graining. Concurrently, the corresponding continuum Hamiltonian will no longer be of the form Eq. (22) but instead belongs to the class of "SMF decomposable" Hamiltonians, 10 that are related to classical master equations and include the RK-Hamiltonians Eq. (22) as a subclass [72]. This correspondence will prove useful in deriving the continuum expression for the RK-Hamiltonian, which we discuss in Sec. V C (see also Appendix C).
To close this discussion, we note that expressions of the form Eq. (46) have previously been derived for the continuum limit of ground states of RK-Hamiltonians using various methods, albeit never in the context of multipole conserving systems. For instance, the exponent in Eq. (46) can be interpreted as the free energy functional corresponding to a configuration of the height field φ (m) (x), as is typically done in the context of RK points in dimer models [71,83]. Alternately, the expression Eq. (46) can also be derived using the path-integral formulation 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 the ground state wave functional, we now identify the corresponding expression for the coarse-grained RK-Hamiltonian H RK . As discussed in the previous section, the coarse-grained wave functional Eq. (46) is the ground state of a multipole conserving RK-Hamiltonian, which belongs to the generalized class of frustrationfree positive-definite RK-like Hamiltonians discussed in Ref. [72]. In Appendix D, we discuss two distinct approaches for deriving the continuum limit of H RK : the first approach involves an appropriate choice of regulators, which allows us to explicitly obtain the continuum parent Hamiltonian corresponding to Eq. (46). The second, more commonly employed approach [71,82,83,85,86] exploits the relationship between H RK and classical master equations discussed in Sec. III (see Eq. (16)). In summary, this approach proceeds by identifying the classical process corresponding to H RK which equilibrates to a Gaussian distribution of height fields, as given by Eq. (46). As we show in Appendix D 2, this classical process describes the Langevin dynamics of the generalized height fields under damping. The continuum expression for H RK can then be derived via the Fokker-Planck equation for the probability functionals of the generalized height fields.
Both approaches lead to the same continuum expression for H RK , which is the parent Hamiltonian for the GS wave functional Eq. (46) and is given by where γ is an overall dimensionful constant. The creation and annihilation operators Q † m (x) and Q m (x) are defined as and satisfy the commutation relations We can directly verify that the wave functional Φ Up to a constant (infinite) energy shift, the continuum Hamiltonian Eq. (48) can be brought to more standard form [82,85] We observe that the Hamiltonian Eq. (53) is invariant under a polynomial shift of the form where P (2m+1) (x) is a polynomial in x of degree ≤ (2m + 1). That is, it has additional symmetries beyond just the m th multipole moment, as is typical of continuum RK-Hamiltonians [81]. However, using Eq. (38), it is straightforward to see that the transformation Eq. (54) changes the quantum number sector of the system. This shows that the continuum Hamiltonian Eq. (53) is the same across all quantum number sectors, further implying that the ground state sector is extensively degenerate. Now that we have established the form of the continuum Hamiltonian, we can study its lowest energy excited states to derive the dispersion relation and the gap. Using Eqs. (48) and (50), the excited states Φ k [φ (m) ] can be written as (55) where the mode function f (m) (kx) is determined by the boundary 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 For a finite system of size L, we thus expect the gap ∆ (m) of H (m) to scale as For charge-conserving systems (m = 0), we see that the continuum height field approach correctly reproduces the scaling of the spin-wave dispersion relation of the Heisenberg model discussed in Sec. IV. More generally, we can further lower-bound the scaling of the Thouless time t (m) Th for a system of size L conserving the m th multipole moment as follows: Due to the polynomial shift symmetry (Eq. (54)) of the continuum Hamiltonian, we expect that this scaling of the Thouless time is independent of the quantum number sector. Eq. (58) is one of the main results of this paper as it encodes the subdiffusive scaling of the Thouless time in systems with higher moment conservation laws. These results, obtained analytically through the generalized height representation developed herein, are validated by the numerical analysis performed on dipole conserving FRQCs (see Sec. IV). The applicability of our continuum analysis of H RK extends beyond the context of random quantum circuits and is directly pertinent to the study of classical cellular automata with conserved higher moments. Such automata were studied in Refs. [63,64] and are equivalent to the circuit M . To further test the validity of the continuum Hamiltonian obtained in Eq. (53), we can compute the two-point spin correlations using Eq. (36): where , and F is a hypergeometric scaling function. Eq. (59) is in agreement with scalings obtained from numerical calculations and hydrodynamic considerations in Ref. [64].

VI. HIGHER DIMENSIONAL CIRCUITS
We now briefly discuss extensions of our results to constrained FRQCs in dimensions d > 1, and in particular, systems on a hypercubic lattice that conserve all components of the m th multipole moment. First, we note that the discussions in Secs. II and III generalize directly mutatis mutandis to higher dimensions.
We start with a d-dimensional spatially-random FRQC W acting on a set of sites carrying color and spin DOFs, with the local Hilbert space given by H loc = C q ⊗ C 2s+1 . The circuit W takes the form of Eq. (2), comprising several layers { W a } composed of local unitary gates U [·] . The layers of { W a } are arranged in a "Trotterized" form: (i) For a given { W a }, U [·] commute with each other, and all sites are being acted upon by exactly one gate; and (ii) each group of neighboring sites will be acted by a U [·] in some { W a } in W once and only once. An example of a two-dimensional system with charge conservation will be provided below in Eq. (62).
As before, we impose symmetries or local constraints on the spin DOFs and take the large-q limit in the color DOFs; thus, the local gates U [·] have the block-diagonal forms shown in Eq. (3). Using techniques directly generalized from Appendix A, we find that in the q → ∞ limit, the SFF is expressed as Eq. (10), where M is a bistochastic matrix that retains the geometry of the original circuit W but with its unitary gates U In what follows, we will be interested in d-dimensional systems that conserve all components of the m th multipole moment. We also restrict ourselves to hypercubic lattices with OBC in all directions, with coordinates labelled by a d-dimensional vector x = (x 1 , · · · , x d ). The m th multipole moment operators are given by rank-(m + 1) symmetric tensors Q i1···im m , defined as where S z x is the Pauli-Z matrix acting on site x, the indices of the tensor {i j } (1 ≤ i j ≤ d) represent the d lattice directions, and the summation runs over all sites of the hypercubic lattice. Note that when d = 1, we recover Eq. (25). Quantum numbers associated with the operators { Q i1···im m } will be denoted by {Q i1···im m }. For example, the expressions for charge (m = 0) and dipole moment (m = 1) are We now illustrate the above with an example and calculate t Th for an FRQC composed of charge-conserving gates acting on spin-1/2 DOFs living on neighboring sites of a square lattice (with OBC). The circuit W in this case can be implemented in four layers: where each of the m [(x,y),(x+α,y+β)] is a 4 × 4 matrix that has the form shown in Eq. (26). Following Eqs. (17) and (27), the corresponding RK-Hamiltonian is the spin-1/2 ferromagnetic Heisenberg Hamiltonian in two dimensions: where i, j represents nearest-neighboring sites on the square lattice. Similar to the one-dimensional case, the Hamiltonian Eq. (64) has a ferromagnetic ground state and its lowest energy excitations can be solved exactly; these are known to be spin-waves with a dispersion relation (k x , k y ) = 2 sin 2 (k x /2) + 2 sin 2 (k y /2), where k x and k y represent the momenta of the spin wave in the x and y directions respectively. Furthermore, the Hamiltonian Eq. (64) is SU (2) symmetric so that the low-energy spectrum is the same within any of the Q 0 sectors. The gap within any S z sector thus scales as (if L x > L y ) Following Eq. (20), the Thouless time for a charge conserving system hence scales with the square of the longest linear-size of the system, consistent with expected results from diffusion. This discussion generalizes directly to charge conserving FRQCs acting on d-dimensional hypercubic lattices, where the emergent RK-Hamiltonian is the ferromagnetic Heisenberg Hamiltonian in d-dimensions with spin-wave excitations and the Thouless time scales as the square of the linear-size of the system. For dipole and higher multipole moment conserving systems, in general, or for charge-conserving FRQCs with higher spins or longer-range gates, the emergent RK-Hamiltonian is generically non-integrable. Similar to the one-dimensional case, we hence consider systems with weak fragmentation [14], take the continuum limit and resort to field theoretic arguments to obtain the gap scaling of the resulting Hamiltonians.
Recall that the ground state of an RK-Hamiltonian is an equal superposition of all configurations within a given quantum number sector, similar to the one-dimensional case (see Sec. V B and App. B); in the continuum limit, the ground state wavefunctional in d-dimensions is then where ρ( x) is the charge density, Z is a normalization factor, and G[ρ( x)] enforces the global symmetry constraints i.e., it fixes the quantum numbers of the sector of interest. For example, for a dipole conserving system in d-dimensions with quantum numbers {Q 0 , {Q j 1 }}, we have (67) We now need to derive a continuum parent RK-Hamiltonian for the wavefunctional Eq. (66). To circumvent the global constraints in Eq. (66), we need some analog of the generalized height fields that we had introduced for one-dimensional systems in Sec. V A. As emphasized in that section, the key role of the generalized height fields is to translate the global constraints in the wavefunctional into boundary constraints. As shown in Eq. (38), for m th multipole conserving 1D systems in the continuum, this was accomplished by demanding that the generalized height fields φ (m) (x) satisfy the generalized Gauss law Eq. (36). The natural analog of the m th generalized height fields φ (m) (x) in higher dimensions are given by symmetric rank-(m + 1) tensor fields {E j0···jm ( x)}, versions of which have previously been studied in the context of fracton models [87][88][89][90].
To recast the global symmetry constraints enforcing the conservation of the n th multipole moments (n < m) in terms of boundary constraints on the tensor fields, we impose the following generalized Gauss law on the rank-(m + 1) tensor fields: where we sum over repeated indices. We can then express the conserved quantities {Q i1···im m } in terms of boundary constraints on the tensor fields. For example, in charge conserving systems (m = 0), Eq. (68) reduces to the usual Gauss law for electric fields ∂ j0 E j0 ( x) = ρ( x), and the total charge Q 0 can be expressed as where dn j0 represents the "area" element on the boundary of the system, and we have used integration by parts along with Stokes' theorem. Similarly, in dipole conserving systems, Eq. (68) reduces to the generalized Gauss law for rank-2 symmetric tensor fields ∂ j0 ∂ j1 E j0j1 ( x) = ρ( x) and the total charge Q 0 and dipole moments {Q i 1 } can be expressed as [67] It is straightforward to show that a general expression for the n th multipole moment can also be derived in terms of boundary integrals of rank-(m + 1) symmetric tensor fields {E j0···jm ( x)} for any m ≥ n, although the general expressions are rather tedious to show here and are not particularly illuminating. Thus, for a system with m th multipole moment conservation in all directions, we work in terms of rank-(m+1) symmetric tensor fields, with the ground state wavefunctional Eq. (66) re-expressed as where B[{E i0···im ( x)}] represents a boundary constraint on the fields {E i0···im ( x)} that fixes the quantum number sectors corresponding to all the n th multipole moments for n ≤ m.
We now proceed to derive the expression for the parent RK-Hamiltonian corresponding to the wavefunctional Eq. (71). The derivation closely follows the onedimensional case discussed in Sec. V C. The crucial idea is that in the long-wavelength limit, the Markov process corresponding to the RK-Hamiltonian of an m th multipole conserving system is simply the independent Langevin dynamics of each component of the rank-(m+1) tensor field at each point. We can intuitively understand this on a two-dimensional square lattice, where we label the two directions byx andŷ. The generalized Gauss law of Eq. (68) is then discretized appropriately, and acts locally around each site of the lattice. In charge conserving systems, the rank-1 electric field E i has two components E x and E y , which can be thought of as DOFs on the links of the lattice along thex-andŷ-directions respectively (see Fig. 6a). As a consequence of the discrete Gauss law, any nearest-neighbor charge conserving process along a link in thex (resp.ŷ) direction only modifies the fields E x (resp. E y ) on that link, whereas the electric fields far away remain unchanged. In the continuum, such processes are modeled by the independent Langevin dynamics of E x and E y on each link. Similarly, in dipole conserving systems on a lattice, the rank-2 symmetric tensor field has three independent components: E xx , E yy , E xy = E yx . The components E xx and E yy are DOFs on the vertices of the square lattice whereas E xy = E yx are DOFs living on plaquettes of the square lattice (see Fig. 6). As a consequence of the discrete generalized Gauss law, various local dipole conserving processes that occur independently result in independent fluctuations of these tensor fields. Furthermore, after coarse graining, we expect that the fluctuations in each component of the local fields will be Gaussian and that the fluctuations of different components of the tensors E j0···jm will be uncorrelated.
Using the expression Eq. (71), (in Appendix E) we derive the following expression for the continuum Hamiltonian: where repeated indices are summed over, and (Q † ( x)) l0···lm and (Q( x)) l0···lm are respectively creation and annihilation operators for the fluctuations of the component E l0···lm ; their explicit expressions are given by Eq. (E12). Note that we obtain separate creation and annihilation operators for each component of the rank-(m + 1) tensor since their fluctuations are independent. Further, using the properties of these operators shown in Eq. (E12), the lowest excited state Φ where repeated indices are summed over, and we have supressed the arguments in Φ can also be shown to satisfy For a system with linear-size L j in the j th direction, we thus expect the gap to scale as (assuming L = max j (L j )) thereby showing that the Thouless time follows the scaling of Eq. (58) i.e., the Thouless time for multipole conserving circuits in higher dimensions follows the same subdiffusive scaling with the linear extent of the system as that of one-dimensional multipole conserving circuits. While we have primarily focused on systems that conserve all components of the m th multipole moment, this formalism directly generalizes to systems where only a few components of m th multipole moments are conserved. Such a setting is directly relevant to many physical systems, for instance in recent experiments that impose dipole moment conservation only along a single direction by subjecting the system to a strong electric field in that particular direction. Continuum wavefunctions of the form Eq. (71) for such systems can also be expressed in terms of tensor fields that obey anisotropic versions of the Gauss law Eq. (68) [89]. For example, in a two-dimensional system with charge conservation in the x-direction and dipole moment conservation in the y-direction, we obtain Following similar ideas as in the isotropic case, it is then straightforward to derive expressions for the continuum Hamiltonian similarly to Eq. (72), which corresponds to Langevin dynamics of each of the tensor components involved, and to then derive the scaling of the Thouless time. We find that the Thouless time for the entire system is dominated by the highest multipole moment conserved, i.e. t Th ∼ L 2(m+1) if some component of the m th multipole moment (but none higher) is conserved, consistent with intuition and experimental observations [65].

VII. CONCLUDING REMARKS
In this paper, we have studied the spectral statistics, as encoded in the SFF K(t), for spatially-extended constrained many-body quantum chaotic systems, focusing on FRQCs with conserved higher moments, such as the dipole moment. As one of the key results of this paper, we have established a series of relations between K(t) in the q → ∞ limit, a classical stochastic circuit M , and an emergent RK-Hamiltonian, such that the inverse gap of this RK-Hamiltonian lower bounds the Thouless time t Th of the underlying FRQC. As we have shown here, the relation between t Th and ∆ RK proves particularly efficacious, since it relates a dynamical property of the FRQC to the low-energy physics of a sign-problem-free quantum Hamiltonian.
We emphasise that these relations are valid for generic local FRQCs with on-site Abelian symmetries or dynamical constraints, not only those with conserved higher mo-ments of charge. For example, we can consider circuit implementations of other fragmented models [91] or study an FRQC inspired by the Rydberg blockade [56,92], also known as the PXP model [8]. The latter is implemented by taking e.g., = 3 site local gates with the only nontrivial dynamics contained within a 2×2 block connecting the |↓↓↓ and |↓↑↓ states. The resulting Floquet operator W has no conserved quantities besides the (quasi)energy, but fragments into dynamically disconnected subspaces; the largest of these corresponds to the constrained Hilbert space most often discussed in the context of quantum many-body scar dynamics [8].
We have verified these general results on circuits with higher conserved moments, which generically exhibit Hilbert space fragmentation. Working in the q → ∞ limit, we derived the corresponding stochastic circuit M and emergent RK-Hamiltonian H RK for both charge and (weakly fragmented) dipole conserving systems. Our numerical study of these systems suggests a universality in the scaling of t Th with system size, specifically, we predict diffusive scaling t Th ∼ L 2 for charge conserving systems and subdiffusive behavior ∼ L 4 for dipole conserving systems, regardless of the microscopic details of the underlying circuit. Further evidence for this scaling is given by continuum field theoretic descriptions of the emergent RK-Hamiltonians for multipole conserving FRQCs in terms of generalized height fields. By analytically computing the dispersion relation for the resultant field theories, we find that t Th ∼ L 2(m+1) in circuits that conserve the m th multipole moment, consistent with numerical results for charge and dipole conserving systems. We further generalize our formalism to higher dimensions, where we derive continuum field theories for emergent RK-Hamiltonians for systems that conserve dipole and higher multipole moments. We obtain the same scaling of the Thouless time with the largest linear size of the system i.e., t Th ∼ L 2(m+1) for circuits that conserve any component of the m th multipole moment (but none higher) in any number of dimensions, consistent with expectations from the one-dimensional result.
Our work opens many exciting avenues for future research: here, we have only focused on the class of multipole conserving circuits which exhibit weak fragmentation. Dynamics in strongly fragmented systems, where typical initial states are ETH-violating, is highly constrained; nevertheless, such systems exhibit large Krylov subspaces which eventually thermalize [16]. The scaling of t Th within such subspaces remains to be understood and may lead to distinct continuum field theories than those we have introduced for weakly fragmented systems. Another interesting avenue to explore is extending our formalism to incorporate non-Abelian symmetries, for which the nature of transport and thermalization is currently being debated [93][94][95].
We note that the large-q diagrammatics, and therefore the mapping to a classical bistochastic circuit and RK-Hamiltonian, have so far only been developed for the two-point SFF K(t). Other observables, such as the sec-ond Renyi entropy and out-of-time-order correlator, can also be mapped to stochastic classical dynamics upon ensemble averaging, and will be discussed in forthcoming work. Pushing these ideas further presents an important but technically-demanding theoretical challenge. More straightforward is extending our results to circuit geometries besides the brick-wall structure considered here as well as to other RMT symmetry classes.
More pressing, however, is building a systematic understanding of FRQCs at finite-q, to delineate those features which are an artefact of the q → ∞ limit from those which are more generic properties of constrained random circuits. Numerically investigating finite-q circuits remains prohibitive, particularly in the context of higher moment conserving circuits which already require large ( ≥ 4) local gates. Analytically, one could attempt to keep track of diagrams at next to leading order in the large-q expansion to better quantify deviations of the SFF 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 enlightening discussions. We also acknowledge useful con- In this Appendix, we generalize the calculation of SFF K(t; K) performed in Ref. [48] for an FRQC with conserved U(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. 7(a), can be evaluated as a sum of diagrams, each of which corresponds to a possible pairing between unitaries and their complex conjugates. In the limit of large-q, diagrams with the same "cyclical" pairing at all sites dominate the sum (see illustration in Fig. 7(c) and Ref. [46] for details). There are t such diagrams, which we call "Gaussian" since they can be evaluated using the Wick contractions between unitaries and their conjugates. For -site unitaries U and their conjugates U * , shown in Fig. 7(d), we have where i = (i 1 , i 2 , . . . , i ) and α = (α 1 , α 2 , . . . , α ) are -component vectors denoting respectively the color and spin degrees of freedom on their support. Here, d α denotes the number of "out-going" -site spin configurations dynamically accessible to the "in-coming" -site spin configuration α. That is, d α q is the size of the random matrix block in W to which the basis state ( i, α) belongs.
As an example, let us consider the simplest charge-conserving circuit [48]. We have two-site gates ( = 2) and α = (α 1 , α 2 ), with α i =↑, ↓. The dynamically connected two-site configurations can be labelled by a conserved charge Q 0 = x=1 S z x . The configurations (↑, ↑) and (↓, ↓) have Q 0 ( α) = +1 and Q 0 ( α) = −1 respectively and are not dynamically connected to any other two-site configurations. Thus, these correspond to d α = 1. However, the configurations (↑, ↓) and (↓, ↑) both have Q 0 (α) = 0, and are dynamically connected to each other. Thus we have d α = 2. Next, we can translate each of the t diagrams into algebraic terms by using Eq. (A1) and summing over the color degrees of freedom. Consider the diagram where the pairing between unitaries and their conjugates takes the form of the leftmost diagram in Fig. 7(c) at every site. The sum over the colors precisely cancels out all factors of q in Eq. (A1) [46]. It then remains to sum over the spin degrees of freedom. Observe that the choice of spin DOFs in W uniquely fixes those in W † due to the Wick contractions Eq. (A1). Consequently, the sum over spins can be computed by finding all possible ways of assigning spins in the diagrammatic representation of Tr K W (t) (Fig. 7(b), such that all charges {Q a } are preserved after the action of every -site gate. This sum over spins can be reproduced by Tr K M t in Eq. (10), where M has the geometry of W in Eq. (2) and contains non-random block-diagonal -gates m, defined in Eq. (9). Note that m contains d α × d α matrix blocks m(d α ) with entries 1/d α to account for the factor d α in Eq. (A1).
The overall factor of |t| in Eq. (7) arises from the following argument: any two of the t contracted diagrams in Fig. 7(c) are related by a rotation of the arrowed loop on the right side of one of the diagrams i.e., the diagrams are topologically equivalent. As a result, the corresponding algebraic terms for all t leading diagrams are identical. This concludes 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 discrete setting, the ground state wave function Eq. (23) has the form where Λ denotes the set of allowed spin configurations that are allowed within the quantum number sector or Krylov subspace 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 variables with zero mean and variance σ 2 ∼ ∆x as a consequence of the central limit theorem [84]. That is, the probability density P of a single spin s n is P (s n ) ∼ exp − s 2 n 2σ 2 =⇒ s n = 0, s n s n = σ 2 δ n,n , σ 2 = ∆x κ .
Here, κ effectively serves as the coarse-graining length scale and is determined by requiring that the continuum theory correctly reproduces long-distance correlation functions in the microscopic model, which can in principle be computed numerically. Next, by expressing the s n 's in terms of the m th height fields and by using Eqs. (29) and (31), we obtain while the joint probabilities satisfy n }] are shorthand for global and boundary constraints written in terms of spin and height variables respectively; these constraints are required to fix the quantum number sector, as discussed in Sec. V A. Note that the Jacobian of the transformation from {s n } → {φ (m) n } is 1. In order to obtain a continuum limit in the case of m th height variables, we need to define height fields that scale with the lattice spacing as shown in Eq. (33). The probability density for a configuration of the m th height field in the continuum limit then reads from which we obtain Eq. (46) as the continuum limit of the ground state wave functional. The full continuum wave function can be written as

Regulator approach
Let us begin by observing that the continuum wave function Eq. (B6) has the form Eq. (C1) i.e., it belongs to the class of ground states of SMF decomposable Hamiltonians. We thus anticipate that the continuum limit of the RK-Hamiltonian is SMF decomposable and can be brought to the form Eq. (C2).

Fokker-Planck approach
An alternate, more standard approach for deriving the continuum Hamiltonian proceeds by invoking an analogy with Langevin dynamics [71]. We begin with the observation that that the ground state wave functional of the continuum RK Hamiltonian-when interpreted as the equilibrium probability distribution of the Markov process in the continuum-given by the absolute square of the amplitudes in Eqs. (B6) and (D1), resembles a Boltzmann distribution where κ is the "inverse-temperature" and E(φ (m) (x)) is the energy of the system given by Given that the underlying Markov process relaxes to the Gibbs distribution, we expect the long-wavelength behavior at late-times to be captured by Langevin dynamics of the generalized height fields φ (m) (x, t) [71,96] dφ (m) (x, t) dt = − γκ 2 δE(φ (m) (x, t)) δφ (m) (x) + ζ(x, t) .
To see that the distribution Eq. (D13) is indeed the equilibrium distribution of Eq. (D15), we derive the timeevolution equation for the joint probability densities, also known as the Fokker-Planck equation. We follow methods elucidated in Refs. [71,74,96,97]. Note that although we consider OBC in the main text, we will use PBC in the following for convenience. We first define Fourier transforms of the height fields {φ (m) k } as follows where k is a discrete momentum variable, and we have suppressed the dependence on t and without loss of generality set φ (m) k=0 = 0. Taking the the Fourier transform of Eq. (D15), we then obtain Further, using Eqs. (D16) and (D18), we obtain ζ k (t) = 0, ζ k (t)ζ k (t ) = γ δ k ,−k δ(t − t ) . (D18) which obey the commutation relations [Q m (x), Q m (y)] = 0, Q † m (x), Q † m (y) = 0, Q m (x), Q † m (y) = (−1) m+1 κ ∂ 2(m+1) x δ(x − y).
Consequently, the ground state wavefunction is annihilated by Q m (x) and is given by Furthermore, following Eq. (D34), the excited wavefunction Φ In this Appendix, we briefly sketch the derivation of the continuum Hamiltonian for m th multipole conserving systems in d-dimensions, following the Fokker-Planck approach illustrated in Appendix D 2. As discussed in Sec. VI, the ground state wavefunctional of the continuum RK-Hamiltonian, when interpreted as the equilibrium probability distribution of the Markov process in the continuum, is given by the absolute square of the amplitudes in Eq. (66): where κ is the "inverse-temperature" and E(E j0···jm ( x)) is the "energy" of the system, given by where repeated indices are summed over. Similar to Eq. (D15) in the one dimensional case, we expect the longwavelength behavior at late-times to be captured by Langevin dynamics of the generalized tensor fields E j0···jm ( x, t): where ζ l0···lm ( x, t) is δ-correlated "white-noise" governing the dynamics of the component E l0···lm ( x, t) of the tensor field i.e., it satisfies a Gaussian distribution with where δ {l0···lm},{p0···pm} is 0 unless the sets {l 0 , · · · , l m } and {p 0 , · · · , p m } are equal -this ensures that the fluctuations of different components of the tensor field {E j0···jm } are not correlated, as discussed in Sec. VI. Considering a system in d-dimensions with linear-size L in each direction, we define Fourier transforms of the generalized tensor fields {E j0···jm k } as follows where repeated indices are summed over, k is a discrete momentum variable, and we set E j0···jm