Localization in fractonic random circuits

We study the spreading of initially-local operators under unitary time evolution in a 1d random quantum circuit model which is constrained to conserve a $U(1)$ charge and its dipole moment, motivated by the quantum dynamics of fracton phases. We discover that charge remains localized at its initial position, providing a crisp example of a non-ergodic dynamical phase of random circuit dynamics. This localization can be understood as a consequence of the return properties of low dimensional random walks, through a mechanism reminiscent of weak localization, but insensitive to dephasing. The charge dynamics is well-described by a system of coupled hydrodynamic equations, which makes several nontrivial predictions in good agreement with numerics. Importantly, these equations also predict localization in 2d fractonic circuits. Immobile fractonic charge emits non-conserved operators, whose spreading is governed by exponents distinct to non-fractonic circuits. Fractonic operators exhibit a short time linear growth of observable entanglement with saturation to an area law, as well as a subthermal volume law for operator entanglement. The entanglement spectrum follows semi-Poisson statistics, similar to eigenstates of MBL systems. The non-ergodic phenomenology persists to initial conditions containing non-zero density of dipolar or fractonic charge. Our work implies that low-dimensional fracton systems preserve forever a memory of their initial conditions in local observables under noisy quantum dynamics, thereby constituting ideal memories. It also implies that 1d and 2d fracton systems should realize true MBL under Hamiltonian dynamics, even in the absence of disorder, with the obstructions to MBL in translation invariant systems and in d>1 being evaded by the nature of the mechanism responsible for localization. We also suggest a possible route to new non-ergodic phases in high dimensions.

The quantum dynamics of interacting many-body systems is a great open frontier for theoretical physics. Important advances on this front over the past decade include the development of the theories of many-body localization (MBL) [1][2][3][4] and the eigenstate thermalization hypothesis (ETH) [5][6][7], and advances in our understanding of the scrambling of information in quantum chaos [8][9][10]. While these advances have occurred in the context of closed quantum systems with time translation symmetry (either continuous or discrete), more recent advances have also shed light on our understanding of quantum dynamics subject to noise i.e. without time translation symmetry, and hence without the constraints imposed by energy conservation. These advances stem from the exploration of many-body dynamics in random unitary circuits 1 , where the time evolution is generated by the application of random gates [12][13][14][15][16][17][18]. (Similar work has also been done in the context of FIG. 1: Random unitary circuit: each site is a three-state qudit. Each gate (blue box) locally conserves S total z , the total z component of the three qudits it acts upon, and P total , the total dipole moment of the three qudits. It is a block-diagonal unitary of the form shown on the right, with each block (red boxes) of each gate independently Haar-random. The smaller 1 × 1 blocks do not flip the qudits, while the larger 2 × 2 blocks produce charge-and dipole-conserving qudit flips as indicated.
Floquet circuits, featuring periodic time evolution [19][20][21][22][23][24].) These explorations have provided an unprecedentedly detailed understanding of the onset of many-body quantum chaos. However, all random circuit models that have been studied so far thermalize 2 , and it remains unclear whether noisy quantum dynamics can lead to robust non-ergodic phases analogous to MBL.
Parallel to the developments in quantum dynamics, recent years have also witnessed an explosion of interest in fracton phases of quantum matter [26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42][43][44][45] -for a review, see [46]. Fracton phases are quantum phases of matter in which the elementary excitations exhibit restricted mobility, being either unable to move under local Hamiltonian dynamics, or able to move only in certain directions. This restricted mobility can be traced to the superselection structure of the underlying theories, which in addition to familiar constraints like charge conservation, imposes also additional "higher moment" constraints such as conservation of the dipole moment of charge. Fractons have deep connections with numerous areas of physics, such as elasticity [47][48][49], gravity [50], holography [51], quantum Hall systems [52], and deconfined quantum criticality [53]. Fracton phases are also known to exhibit glassy dynamics. These have been studied [26,[54][55][56][57] for fracton systems in three spatial dimensions, using techniques akin to locator expansions. At non-zero densities and long enough times, "locator expansion" type techniques fail to converge, likely indicating thermalization at long times. Whether the study of fractons can lead to the discovery of new dynamical quantum phases remains an open question.
In this work, we combine the techniques of random unitary circuits and ideas from the study of fractons, and discover an entirely new dynamical phenomenon. Specifically, we consider a circuit that conserves not only a U (1) charge (as in [17,18]), but also the dipole moment of that charge, as an analogue of the "higher moment" constraints that underlie the fracton phenomenon. Insofar as energy conservation is essential to immobility of fractons under Hamiltonian dynamics [56], and random circuits do not conserve energy, one might have thought that random circuits with fractonic constraints would also thermalize, just like all other random circuit models that have been studied thus far. Instead, we find that in such circuits, fractonic charge fails to spread even at the longest times -even though the circuit does contain gates that allow such charge to hop. This unexpected immobility of fractonic charges is traced to the inevitability of returns in low dimensional random walks, akin to the phenomenon of weak localization [58], but immune to dephasing and not limited to the non-interacting limit. A system of coupled hydrodynamic equations is proposed to govern the charge dynamics, and makes several nontrivial predictions which are found to be in good agreement with numerics. Importantly, this system of equations predicts that the localization we observe should persist also to two dimensions, where exact numerical simulations of quantum dynamics are unavailable. At the same time, the stationary fractonic charge emits "non-conserved operators", which do spread, but with exponents distinct to those that have been previously observed in random circuits [17,18], constituting a new dynamical universality class. These non-standard exponents are also correctly predicted by analysis of our hydrodynamic equations. The non-ergodic phenomenology is numerically verified to persist in settings with a non-zero density of fractons, including in the sector with close to maximal charge. Thus we establish that combining "fractonic" constraints with random circuit dynamics can lead to qualitatively new classes of dynamical behavior. At the same time, insofar as random circuit dynamics is a model for noisy time evolution, our results suggest that low dimensional fracton systems should behave as memories forever robust against noise. 3 Equally remarkably, however, it offers a pathway to breakthrough results in many-body localization.
The remarkable developments in MBL over the past decade have largely relied for analytical understanding on techniques related to locator expansions. A series of works have established that there exist rare region "obstructions" to many-body locator expansions [59][60][61][62][63][64][65] in any setting other than strongly disordered one dimensional spin chains. While some of these obstructions can be circumvented [66][67][68] it does appear that systems in dimensions greater than one cannot admit of a fully convergent locator expansion [64], and nor can systems with a translation invariant Hamiltonian [65]. As a result, it is widely believed that higher dimensional systems and/or translation invariant systems can support at best asymptotic or "quasi" MBL, in which a memory of the initial conditions survives in local observables for times superpolynomially long in some control parameter [64], but the system nevertheless thermalizes on the longest timescales. Our results provide an explicit counterexample to this belief, since they indicate that fractonic systems in one and two dimensions preserve forever under Hamiltonian time evolution a memory of their initial conditions in local observables, even with translation invariant Hamiltonians. That is, translation invariant one and two dimensional fractonic systems undergoing Hamiltonian dynamics should exhibit true MBL (not just asymptotic), at least in the sector of Hilbert space with non-zero fracton charge, and this should survive to non-zero energy densities and should be robust to generic local perturbations obeying the fractonic constraint. 4 Our result evades the "no-go" arguments of [64,65] because it does not rely on locator expansions for localization. As such, it provides a first example of MBL beyond locator expansions, and may open a new chapter in the exploration of that field.
The paper is structured as follows. In Sec. II we introduce the basic model that will be studied. In Sec. III we present the results of numerical simulations indicating the lack of spreading of fractonic charge. In Sec.IV we present the basic analytical understanding of random circuit dynamics with fractonic constraints, and test it against numerical simulations. In Sec. V we discuss the implications of our results, and conclude.

II. THE SPIN-1 RANDOM UNITARY CIRCUIT MODEL
We consider a one-dimensional chain of L sites, each of which has a spin-1 degree of freedom. The gates in the random unitary circuit at each step in the time-evolution are designed to preserve the total z component of the spins (S total z ), which serves as a conserved U (1) charge, as in Refs. 17 and 18. Unlike previously studied systems, however, we consider a random unitary circuit which also conserves the total dipole moment ( P total ) of this effective charge. As such, an isolated charge (i.e. a state with an S z = 1 in a background of S z = 0) cannot propagate without creating additional dipolar excitations -and as such functions as a fractonic charge.
The time-evolution is governed by a random unitary circuit (as shown in Figure 1) consisting of staggered layers of three-site unitary gates. The time evolution unitary is given by As a result of the conservation laws, each three-site unitary gate U i,i+1,i+2 is a 27 × 27 block-diagonal matrix, with block structure as shown in Figure 1. Each block is a Haar-random unitary of the appropriate size, and is chosen independently of other blocks. There are four 2 × 2 blocks in the unitary matrix, corresponding to the following charge-and dipole-conserving qudit flips: The remaining charge configurations have no nontrivial qudit flips. Despite the small number of gates, this form of time-evolution is the most generic unitary evolution consistent with conservation of charge and dipole moment. As such, any other quantity which is not conserved by design can change under the action of the circuit, indicating that charge and dipole moment are the only two conserved quantities in the system. Thus the system is not a conventionally integrable model, which would have an extensive number of conserved quantities. ). The onsite basis has the following normalization: Tr We use these matrices to form a basis of 9 L generalized Pauli strings S = i Σ µi i . The initially-local operator O o consists only of those strings that are the identity at all sites except a few near the site of initialization. With time, however, there are other basis strings with nonlocal weight in the decomposition of the time-evolved operator O(t). In the string basis, the operator O(t) can be expanded as The basis strings satisfy Tr[SS ]/3 L = δ SS , making it possible to obtain the coefficients as follows: which implies that the total weight of O(t) on all strings S is 1, which is just a statement about unitarity in the dynamics: In order to determine expectation values of operator, consider an initial density matrix where I background is the background state which is initially the identity on the full system; cO o is a traceless local onsite perturbation at the origin to this initial state.

III. CHARACTERIZATION OF OPERATOR SPREADING
In order to study the operator spreading profile of an initially-local operator O o , we follow Ref. [18] and define the "right-weight" ρ R (i, t) as the total weight in O(t) of basis strings that end at site i, meaning that they act as the identity on all sites to the right of site i, but act as combinations of the identity and non-identity operators upto i. This gives From Eq. 3, we see that there is a conservation law on ρ R (i, t), i.e. i ρ R (i, t) = 1. This gives ρ R (i, t) the interpretation of an emergent local conserved density, the spreading of which can then be examined.
A. Diffusive spreading of local conserved dipole moment As a warmup, we first consider the case where the initial operator is a two-site conserved dipole operator, The spin-1 chain (initialized with the above local dipole operator O o ) is time-evolved under the action of the random unitary circuit. We evaluate the exact Heisenberg time evolution of the operator by multiplying it by the appropriate unitary gates, a procedure that is even simpler than exact diagonalization. A further speedup results because the constraints on the unitary operators allow us to restrict to a single charge and dipole sector, such that one does not need to keep track of all 9 L (where L is the system size) coefficients at each timestep. This is repeated for various values of t to obtain the right-weight profile as in Figure 2.
Since the initial condition involves a dipole operator, which acts as a conserved U (1) charge without any conservation laws on its higher moments, this is a situation that should map onto the analysis of [18], and serves as a sanity check on our results. As discussed in [18] the spatial structure of this operator shows a ballistic front, a power-law tail behind the front, and diffusively spreading charges near the origin. This is illustrated in Figure  2, which shows the right-weight profile ρ dipole R (x, t) for a system of size L = 21 initialized with a single conserved dipole. These profiles depict the following regimes: a slowly spreading "lump" at the origin, a "front" that moves out rapidly, and a "tail" behind the front. In subsequent sections we show that the front propagation is ballistic, and that the exponent governing the tail matches the analysis in [18]. It is also apparent that in the long time limit, the central lump relaxes to a uniform charge distribution. These results, consistent with [18], serve as a sanity check on our procedure. We do not discover anything new simply by considering an initial condition with a single conserved dipole. To obtain new results, we must consider initial conditions containing fractons, a task to which we turn in the following subsection.
B. Spreading of conserved fracton charge : memory of initial conditions We now consider the case where the initial operator is a conserved fractonic charge, i.e.
I. Now the situation is unlike [18], since not only is the monopole moment of the initial operator conserved, but so is its dipole moment. The resulting behavior is strikingly altered, as is illustrated in Figure 3. Unlike the dipole right-weight profile, the spatial structure of this operator shows an extremely prominent peak where the operator was initially localized, i.e. the system has memory of its initial conditions in a local observable at all times accessible in our simulations. Additionally, there is a "propagating front" with a lagging tail and, as we shall show in a subsequent section, the power law governing this lagging tail is distinct to the power law observed in Ref. [18] and Figure 2. It is thus clear that the evolution of an initial fracton operator lies in a different universality class. An important point to note is that circuit to circuit fluctuations in the right-weight profiles of both dipole and fracton charge operators are negligible. We show  this in Figure 4 for the case of the fracton charge operator by plotting the absolute value of the difference in the right-weights for six typical runs relative to one of the runs. We also show in the Appendix that these results hold true for a spin-S (S = 1) chain.
The persistent memory of the initial condition is remarkable, since the circuit does contain gates that allow a fracton to move (by creating dipole excitations) -see Table I. One might worry that this might be an artefact of our choice to limit the circuit to three qudit gates, which block diagonalize into at most two by two blocks. To test for this possibility, we look at what happens to the fracton operator peak if we use gates of different sizes, which allow many more transitions. We consider the cases involving 4-and 5-qudit gates (see Table I), and in Figure 5 we see that the fracton peak still remains. While larger gate size results in a broader peak at large t, the integrated area under the fracton peak with respect to the background value at large t, i.e. ω peak (t) − ω peak (∞), is independent of gate size and approaches its asymptotic value as t −3/2 , as shown in Figure 7. We explain this exponent in Section IV B.
As an additional check on the robustness of this localization, we also plot ( Figure 6) the constant background of the right-weight profile (seen at large t) and the integrated weight under the fracton peak as a function of system size. As can be seen from the plot, the weight under the fracton peak persists in the thermodynamic limit, making it safe to attribute this behavior to the fractonic constraints in the system. Thus far we have been considering the "emergent" hydrodynamics of unitary operator spreading (translated into constraints on right-weight profiles). To confirm that the memory of initial conditions is manifest in the dynamics of the physical components of the system, we look at what the charges themselves are doing, and not just the right-weights. We examine what happens to S z as a function of position. The "fracton peak" still remains, and this is largely independent of time and system size, as can be seen from the plots in Figure 8. A plot of the integrated area under the S z peak as a function of system size at large t ( Figure 9) shows that there is considerable operator weight that remains at the original location, i.e. the system acts as a memory which is robust against random noise. Note the absence of any spreading front in this plot: fractonic charge does not spread under the action of the random circuit, even though there are gates that apparently allow charge to move.
Importantly, the peaked late-time S z profile seen in Figure 8 is demonstrably different from the thermal state with the same charge and dipole moment. The thermal state should simply be the state of maximal entropy. To determine what S z looks like for the maximal entropy state, we consider a classical three state model on a length L chain, where every site has charge q = +1, 0, or −1. We then consider the set of all classical states that satisfy the constraints on charge and dipole moment, i.e. circuit model has a prominent peak in the S z profile ( Figure 8). Therefore, the system at late times is indeed non-thermal.

IV. OPERATOR GROWTH IN THE FRACTONIC CIRCUIT
We now develop a basic analytic understanding of the various features uncovered in the previous sections, which we test against numerics. In the case with only local charge conservation (and no dipole conservation) [18], the charge executes a random walk under the action of a random unitary circuit. In our model likewise, the dipoles execute random walks under the action of a random unitary circuit. When we consider a system having fractons, however, the dipole conservation constraint means that these can only move by absorbing/emitting a dipole (eg., + − + ↔ 0 + 0). At first glance, one might think that since gates exist that move fractons (at the cost of creating/annihilating dipoles), the fractons FIG. 7: Log-log plot of the integrated area under the fracton peak in the right-weight plots (ω peak ) as a function of time for three different gate sizes (averaged over 10 runs). On the y axis, we plot the logarithm of the area under the central peak at time t minus the area under the central peak at t → ∞. should also execute a random walk under the action of the random circuit. However, there is a crucial distinction: the fracton remains displaced only if the dipole stays absorbed/emitted. If a fracton moves by emitting a dipole, then a subsequent reabsorption of that dipole will return the fracton to its initial position. The dipoles undergo a random walk. It is well known that random walks in spatial dimension d ≤ 2 always return to the origin. This can be understood in terms of the diffusion propagator, G(x, t) = (4πDt) −d/2 e −x 2 /Dt , where D is diffusion constant of the dipoles, and d is the spatial dimension. In order to determine the return probability of the dipole to its starting point, we consider FIG. 9: Area under the S z peak as a function of system size at large t(= 200) (data averaged over 10 runs). The peak is defined as the region bounded by central site ±2. The uncertainty in the fitting parameter is given in the legend.
the following integral: In d = 1, this integral diverges (∼ √ t), therefore the probability that the dipole returns to its initial position is 1, i.e. the dipole always diffuses back to its starting point. When the dipole returns, this causes the fracton to return to its initial position. Thus, although the fracton can temporarily hop at any given time step, the dipole always eventually diffuses back, returning the fracton to its original position. This explains the persistent peak that we see in the right-weight and S z profiles even as t → ∞.
This simple argument suggests that even though fractonic circuits contain gates that move fractons, nevertheless fractons should be immobile under random circuit dynamics not only in d = 1 (which we have studied numerically), but also in d = 2 (which is beyond reach of our numerical techniques). In d = 3, however, the behavior will qualitatively change, since three-dimensional random walks do not inevitably return to the origin, and dipoles can actually escape to infinity. So in d = 3, fractons should be able to execute a true random walk via permanent emission of dipoles. This suggests that the behavior of fractonic circuits in d = 3 is very different from that in d = 1, 2. Namely, in d = 3 a single fracton evolving under a random unitary circuit should be able to diffuse, but not in d = 1, 2.
These arguments suggest that the fractonic circuit should be well-modeled by a picture of a conserved dipole density coupled to a mobile source/sink (the fracton). The hydrodynamic equations that model the system should then take the form where η is the local dipole charge density at the location of the fracton, D is the dipole diffusion constant, A(t) is a stochastic force (taken to be delta-correlated white noise), and R is the position of the fracton (i.e. the location of the single isolated S z = 1 state). In one spatial dimension, the quantities R, η, A are all scalars.
In higher spatial dimensions, they will be vector quantities. We have written the equations in terms of vector quantities for maximum generality. These equations are expected to describe equally well the system in any spatial dimension, but (as we will show), the character of the solutions is highly sensitive to spatial dimension. We now explore the predictions of this simple picture, and test it against numerics in one spatial dimension.

A. Fracton localization and sensitivity to dimensionality
We now solve the system of equations introduced above, and show that they imply fracton localization in one and two spatial dimensions (but not in three dimensions). We initialize the problem with a single fracton at position zero, and zero dipole charge. The total dipole charge˜ η is related to the position of the fracton R(t) viã However, what enters into Eq. (8) is not the total dipole charge, but rather the dipole charge density at the position of the fracton. To estimate this, we need to estimate how much the dipole density spreads out. We can do this as follows. First, we estimate a typical "return time" (i.e. the typical time that a dipole wanders freely before it returns to the fracton to get reabsorbed). This can be estimated from the equation t 1 G(0, t )dt = C where C is some arbitrary large but finite constant, and G(x, t) is the diffusion propagator. This gives where D is the dipole diffusion constant. The length scale ξ over which the dipole charge density is spread out can then be estimated as ξ ∼ √ Dt return (being cutoff in d > 2 by the scale of the system size L). The density Now Eq. (11) can just be recognized as a Langevin equation, which has the well-known solution where P (R) is the probability distribution for the fracton position, and T is the strength of the stochastic kick, . Now, in d = 1, 2 the coefficient γ is finite and independent of system size, so this predicts exponential localization of the fractons in spatial dimensions d = 1, 2. Meanwhile, in d = 3, γ ∼ L d so the "localization length" is greater than system size, and thus the fracton is effectively spread over the whole system in the steady state. Thus, solving the system of hydrodynamic equations yields the prediction that fracton charge should be localized in one and two dimensions, but not in three dimensions. Note also that the localization length in two dimensions is considerably larger than in one dimension (exponentially large in the diffusion constant of the dipoles, instead of power law large). The localization length is obtained by looking at S z (r) in the long time steady state, and is defined simply as the half width at half maximum of the S z peak. We find a scaling ξ fracton ∼ (gate size) 2 , consistent with the predictions of our analysis.
We now compare the predictions of this analysis to direct numerical simulation in one dimension. Figure 11 shows a semi-log plot of the tails of S z as a function of site position at large t, and the exponential behavior seen is consistent with the above analysis, although we are not able to distinguish between exp(−x) and exp(−x 2 ) behavior numerically.
A nontrivial prediction coming from our analysis of the hydrodynamic equations is that the localization length for fractonic charge should be sensitive to the diffusion constant of the dipoles. This sensitivity is particularly strong in d = 2, where the localization length is exponential in the dipole constant, but even in d = 1 we predict that the localization length for charge should scale as ξ fracton ∼ √ DT . We now test this prediction directly against numerics.
We remind the reader that our model of quantum dynamics is minimally structured, i.e. we are considering a random circuit, made out of gates of a particular size satisfying the charge and dipole moment constraint. Once the gate size is fixed, the parameters D and T are fixed also. The only way to change these parameters is to change gate size. Gate size translates directly into "step size" for the dipoles and fractons, and we therefore expect that D ∼ (gate size) 2 and also T ∼ (gate size) 2 (remember T is just the correlator of the stochastic steps taken by the fractons). This expectation can be straightforwardly checked in numerics, by varying the gate size and looking at the spreading of dipole charge in an initial condition containing a dipole only (see Figure 12) and indeed we find the expected behavior. We therefore predict that charge localization length ξ f racton ∼ √ DT ∼ (gate size) 2 , where the charge localization length could be extracted simply as the half width at half maximum on a plot of S z as a function of position. This is a prediction that can be directly tested against numerical simulations, and we find excellent agreement with our prediction (see Figure 13).
Finally, the analysis can be straightforwardly adapted to initial conditions containing a non-zero dipole chargẽ η 0 . We can think of this initial condition as having "descended" from an initial condition with no dipoles, where the fracton started at position −˜ η 0 . The steady state will then be obtained simply by shifting (13), and will take the form i.e. in steady state the fracton shifts position by −˜ η 0 , thereby absorbing the excess dipole density, but remains localized about this shifted position. This prediction can also be tested numerically. For simplicity we work with periodic boundary conditions, such that there is only a single dipole charge characterizing the system. (Here, left and right are defined with respect to an arbitrary choice of origin, with the anticlockwise direction being "right", and dipole moment is defined mod L). Our solution of the governing hydrodynamic equations (14) predicts that for an initial condition with a fracton and dipole charge η 0 , at late times, the fracton peak will be found shifted left by η 0 i.e. the fracton will end up absorbing all the dipoles. This basic picture is confirmed by direct numerical simulations Figure 14. Note that, in configurations such as Figure 14, the identification of which charges make up dipoles versus individual fractons is inherently ambiguous. However there is no ambiguity in the location of the final peak. We thus conclude that our analysis of the coupled hydrodynamic equations (Eq.8) makes a number of non-trivial predictions which are all in agreement with numerics in d = 1, where high quality numerical tests are possible. As we will show, these equations also accurately describe the universal dynamics of the spreading non-conserved part of the operator (next subsection). This excellent agreement between analytics and numerics gives us high confidence in our analysis, even in two dimensions, where direct numerical simulation is not possible. We recall that our analysis predicts localization in two dimensions (but not in higher dimensions).

B. Propagating fronts and power law tails
For either an initially localized fracton or dipole, the operator spreading profile contains a ballistically propagating front. The front itself behaves in the same way in the fracton and dipole cases, i.e. it propagates ballistically with a similar velocity, featuring front broadening that shows a power-law exponent of 1/2 when plotted against time t. This is due to the fact that the front consists of nonconserved operators emitted at the earliest times and is not influenced by what happens near the origin at late times. This behavior is numerically verified in the plots in Figure 15 and Figure 16.
In addition to the ballistic front in the operator spreading profile, there is important information contained in the structure of the tail behind the front, which behaves differently in the cases of dipole and fracton charge operators (see Figure 17). We observe that the diffusive tails show power-law behavior in both cases and extract the exponents. In the case of the dipole operator, the tails show power-law behavior of the following form: (v B t − x) −3/2 , where v B is the velocity of the ballistic front. This is consistent with the discussion in Ref. [18] and can be understood analogously, i.e. the nonconserved parts of the operator are sourced by the diffusing charge, with strength ∼ dρ c /dt, where ρ c is the weight in the conserved charge sector. But this weight goes as ρ c ∼ 1/ √ t in one dimension, so this naturally gives the tails due to the nonconserved parts emitted at later times the form dρ c /dt ∼ t −3/2 . For the case of an initially localized fracton operator, the tails satisfy power law behavior too, but with exponent −5/2, i.e. the tails have the form (v B t − x) −5/2 . This exponent can also be naturally explained in terms of emission of nonconserved operators from the dipole sector, which in turn is stochastically sourced/sinked by the fracton. That is, fracton motion will generate dipole density, which will then be governed by diffusion close to a sink (the fracton itself, which can reabsorb the dipole). The problem of diffusion close to a sink has been studied previously [69], and it is known that, close to the sink, the resulting diffusion propagator is given by a spatial derivative of the free-space diffusion propagator: (15) which gives the amplitude of a dipole to be at a specific time and place. The total weight in the central peak is then given by: We see that the weight in the conserved sector decays faster, as ρ c ∼ t −3/2 (consistent with our prior observation in Figure 7), and thus by an argument analogous to [18] produces tails that go like dρ c /dt ∼ t −5/2 . 6

C. Entanglement diagnostics
A question closely related to operator hydrodynamics is the problem of entanglement growth. An initial product state develops spatial entanglement as it evolves with time. In systems without quenched disorder, the 6 It is interesting to note that in the absence of dipole conservation [18], an initial condition containing a monopole charge produces tails characterized by exponent −3/2, while an initial condition containing dipole charge produces tails characterized by exponent −5/2. In the presence of dipole conservation (our work) the exponents are inverted -an initial state with monopole charge produces tails characterized by an exponent −5/2, whereas an initial state with dipole charge produces tails characterized by exponent −3/2. entanglement entropy is expected to grow linearly in time, with a growth rate characterized by the "entanglement velocity", v E . This entanglement growth generically continues until the system has reached a maximally entangled thermal state. We here investigate entanglement growth in the fractonic circuit using three different entanglement metrics -the observable entanglement entropy, the operator entanglement, and the entanglement spectrum. We note that the conventional von Neumann entropy is not a suitable measure, because it is identically zero for a pure state, and nor is the conventional entanglement entropy (von Neumann entropy of a bipartition), because according to this measure the initial condition we consider is close to maximal entanglement.

Growth of local observable entanglement entropy
We first work with "observable" entanglement entropy [18], arising from the conserved piece of the initially localized operator. We compute the difference between the typical value of the local observable entropy between two sides of a spatial entanglement cut in our system at time t and the value of local observable entanglement entropy at t = 0, and extract the entanglement velocity from it by looking at how the growth of observable entropy of the spin-1 system scales with time. To define the observable entropy more explicitly, consider the density matrix ρ(t) = (I + O(t))/3 L . The von Neumann entropy of this state is independent of time by unitarity. Now, if we consider the conserved part of this state, i.e. ρ c (t) = (I+O c (t))/3 L , then the observable entropy of this state is S c (t) = −Tr[ρ c (t)logρ c (t)]. We note that while the growth of local observable entanglement entropy scales linearly with t at intermediate times, it saturates at large t. We determine the growth of the local observable entanglement entropy as a function of time for four cuts at location x with x = 5, 7, 10, 15 in an L = 21 spin-1 chain (Figure 18), and see that while the initial scaling of the growth of local observable entanglement entropy is the same in all three cases, the entanglement velocity (∼ 0.16) does not match the front velocity (∼ 0.4) (Figure 15). This is normal and well understood -see e.g. Ref. [17]. Our numerical results also show a distinct difference in the late-time growth of local observable entanglement entropy of the fractonic system versus a system of ordinary conserved charges. For ordinary charge conservation, the final saturation value of the growth of local observable entanglement entropy scales linearly with the size of the partition, indicating a volume law for growth of entanglement, as expected for a thermalizing state. Specifically, the saturation value is very close to maximal observable entanglement. In contrast, the growth of local observable entanglement entropy in a fracton system stops well short of its maximal value, and the saturation value is largely independent of partition size, indicating an area law for the growth of entanglement. Such an area law for the growth of local observable entanglement entropy is consistent with the fact that a low-dimensional fracton system fails to thermalize under random unitary dynamics.

Athermal operator entanglement
We now investigate a different metric for the production of operator entanglement within the framework of spreading operators in the fractonic circuit, which does not separate out the conserved piece of operators. It had been suggested that, according to such a metric, a generic operator rapidly becomes maximally entangled within the region where it is present [70], but we find that this is not the case for fracton operators: both the dipole and fracton charge spreading operators are volume-law entangled, but the late-time entanglement entropy of the fracton charge operator is well below that of a maximally entangled operator.
Our system has 3 states per site: |S z = 0, ±1 . Following the discussion in [71], we may view an operator as a state in a doubled system with 9 states per site. In order to see this mapping between operators and states, we consider a single site operator O ab |a b|, with a, b labeling the three states at each site. The corresponding state is ||O = O ab |a ⊗|b . The notation ||... indicates that this state lives in the doubled system. We can now construct two complete basis sets {Â i } and {B i } which are orthonormal, and consist only of operators with a support on the subsystems A and B respectively (taken to be the two halves of the spin chain in this case). For any linear operator O in the full tensor product Hilbert space of the two operator spaces spanned by {Â i } and {B i }, we have a unique decomposition for O as follows: where O ij = (Â i ⊗B j , O) is given by the inner product on the operator space. For a unitary time evolution operator U (t) (like in each time step of our random unitary circuit), we have U ij = (Â i ⊗B j , U (t)). From this, we construct the operator reduced density matrix We use this to define the operator entanglement entropy [72] S op = −Tr(ρ A,op lnρ A,op ). We note that while the operator entanglement for the fracton charge operator shows volume-law growth as a function of partition size, its late-time value is demonstrably lesser (see Figure 19) than the operator entanglement of the thermal state (a flat distribution over all S z product states with given charge and dipole moment). This entanglement metric thereby serves as yet another independent confirmation that a system initialized with a single fracton fails to reach thermal equilibrium.
We note that the volume-law operator entanglement of the late-time configuration is in sharp contrast with the case of Anderson localization (or many-body localization), for which operator entanglement will remain in an area law at arbitrarily late times. Furthermore, conventional localized phases would not have a ballistically   FIG. 19: Operator entanglement (averaged over 10 runs) for the fracton charge operator in a spin chain with L = 9. The late-time operator entanglement of the corresponding thermal state and for the state evolved without the dipole constraint are also plotted, and these are well above the corresponding value for the fracton charge operator. The "thermal" state here is the equal weight superposition over all S z product states with a given charge and dipole moment. propagating front, since there is no analogue of the emission of nonconserved operators. Rather, conventional localized phases have a logarithmically growing light cone [3], resulting in a much slower spread of quantum information.

Entanglement spectrum
Another important question within a localized phase is the level statistics of the entanglement spectrum, i.e. of the entanglement HamiltonianH defined by ρ = e −H , where ρ is the reduced density matrix for part of the system. We explore the entanglement spectrum by physically bipartitioning our system into a "left" and "right" half, constructing the reduced density matrix for the left half, and then obtaining the spectrum of the entanglement Hamiltonian. In Figure 20, we plot the level spacings of the late-time state for a configuration initialized with a single fracton. In a thermal spectrum, these level spacings would be expected to follow random matrix theory, whereas in a conventional integrable system they would be expected to be Poisson. We find a behavior that is neither random matrix nor Poisson, but which rather follows the semi-Poisson distribution, just like in many body localized systems [73]. The entanglement spectrum level statistics provide further evidence that this system is not thermalized, and nor is it conventionally integrable.

D. Non-zero fracton density
We now turn to configurations having multiple fracton charge operators localized at different sites in the system. Now dipoles emitted during fracton motion do not have to return: they can be absorbed by the other fracton instead. In consequence, the fracton operators can now permanently change their position. For a system with multiple fractons, the diffusion constant for the fractons could be estimated as follows: for fractons a distance l apart to move, they need to exchange dipoles. The diffusion time for dipoles is Dl 2 . Thus the fractons move an amount of the order 1 in a time of order Dl 2 . Thus, the diffusion constant of the fractons would be lower than the diffusion constant of the dipoles by a factor of l 2 , where l is the initial separation between fractons. This identification is supported by the numerical results in Figure 21a.
This basic picture has an interesting corollary: the closer together the fractons are, the more mobile they will be. The random circuit should thus generate an effective "attraction" between fractons (analogous to the discussion in [50]). In the long time limit, therefore, the fractons should agglomerate at their center of mass. This is confirmed by Figure 21b. We further note an interesting "coarsening" dynamics to the system, in which nearby fracton peaks first coalesce together, before all peaks eventually join together at the center of mass of the system. We expect such a coarsening process to lead to an infinite hierarchy of timescales for relaxation, for a finite density system in the thermodynamic limit. We also note that the state being relaxed to is one where the fractons all agglomerate near their center of charge, and this is very different from the classical thermal state (i.e. an equal weight superposition of all S z product states with a given charge and dipole moment, see Figure 10). Thus, the system does not relax to an ergodic state, at least according to the intuitive definition of ergodicity.
The agglomeration of fractons at their center of charge also resolves a conceptual puzzle that may have worried the skeptical reader. Namely, random circuit dynamics allows for the creation from vacuum of dipole-anti-dipole pairs, and their subsequent dissociation into fractonic charges. Why then cannot fractons move by exchanging dipoles with such "spontaneously generated" fractons? The resolution is simple: even while fractonic charges can be created from vacuum, these creation processes cannot affect the total fracton charge or dipole moment, and hence cannot affect the center of charge. Since fracton charge ends up agglomerating at the center of charge, and the center of charge is a conserved quantity, such spontaneous pair production cannot affect the long time steady states of the system.
We now propose a possible mechanism for the interfracton attraction. While the analogy to the "gravitational" attraction between fractons discussed in [50] is appealing, a precise connection is difficult, since the discussion in [50] relies on the existence of a conserved energy. Instead, we propose that this attraction may be of entropic origin, and may be driven by the spreading of the non-conserved part of the operator. The appeal to an entropic mechanism may be surprising, since of course the charge sector does not thermalize, but the "non-conserved" parts of the operator do spread, and presumably tend towards their maximum entropy configuration for a given "charge" profile, and we propose that having the charge agglomerate at its center of mass increases the entropy of the "non-conserved" part of the  operator, even though it does not correspond to the maximum entropy in the charge sector 7 . If we take the total weight (of late-time operator right weights) held in two well-separated fracton peaks (≈ 1, see Figure 3) and consider the situation where the two fractons have come together (Figure 21b), and look at the total weight in the operator right weight peak (≈ 0.4), we see that the total weight in the resulting peak is reduced when the two fractons come together. This suggests that fractons coming together allows the system to "emit" more right weight and thus have higher entropy in the "non-conserved" part of the operator. To further support our claim that clustering at finite fracton density occurs for entropic reasons, we compare late-time operator entanglement values for the different states in question in Table II, and we see that the state where the fractons agglomerate at their center of mass has higher operator entanglement than twice the operator entanglement of a single fracton, but less than that of the corresponding thermal state (consistent with lack of thermalization in the charge sector). A complete understanding The late-time profile shows a peak indicating localization even at finite fracton density.
of the inter-fracton attraction and its consequences, however, is beyond this scope of this work. Finally, we have numerically verified that the fracton localization persists at finite fracton density. In Figure  22 we plot the evolution of S z (r) starting from an initial condition that is close to maximal charge (i.e. packed full of fractons). We still observe localization of fracton charge at late times, despite the fact that this initial condition is most certainly at non-zero fracton density. Similarly, in Figure 23 we consider an initial condition that is packed full of dipoles, and again, we observe localization of fractonic charge at long times. We thus conclude that the localization phenomenon we have discovered here persists to systems at non-zero density, with the long time limit of such systems being a state where charge is localized near the center of charge.

V. DISCUSSION AND CONCLUSIONS
We investigated the operator spreading dynamics of a one-dimensional random unitary circuit with fracton conservation laws, wherein we had not only a conserved U (1) charge, but also the conservation of local dipole moment. We find that even though the circuit contains gates that can move fractonic charge, fractonic charge nevertheless does not move under random circuit dynamics. This remarkable "localization" of fractonic charge FIG. 23: S z (t) for an initial condition in the single charge sector but packed with dipoles. The late-time profile shows a peak indicating localization even in this limit. under random circuit evolution can be traced to the inevitable returns of low dimensional random walks. In this sense, the mechanism is akin to weak localization, except that it is not based on quantum interference and is therefore insensitive to dephasing. There is also some family resemblance to the diffusive returns governing the Altshuler-Aronov correction to conductivity in weaklydisordered metals [74]. We note, however, that our system has neither a conserved energy nor any meaningful notion of density of states, so there is no direct parallel to e.g. [74,75]. We have proposed a set of coupled hydrodynamic equations to describe the dynamics of our problem. Solution of these hydrodynamic equations predicts localization of charge in one and two spatial dimensions, but not in higher dimensions. It also makes a number of nontrivial predictions regarding the steady state distribution of charge for various initial conditions, and its dependence on the range of the gates. These predictions are in excellent agreement with numerical simulations in d = 1, where high quality numerics is available, and we therefore have high confidence in our predictions for higher dimensions, where direct simulation of the quantum dynamics is infeasible. Even while fractonic charge fails to spread under random circuit dynamics, the non-spreading fractonic charge emits a ballistically propagating front of "nonconserved" operators, similar to the case of non-fractonic circuits [18]. The ballistically propagating front broadens as √ t, again much as in non-fractonic circuits. However, the tail behind the front is characterized by a different power law to the case of non-fractonic circuits, and hence falls in a different universality class. These non-standard exponents are also correctly predicted by our hydrodynamic equations, and can be understood as a consequence of the unusual hydrodynamics governing the fractonic circuit, in which conserved "dipolar" charges are coupled to "fractonic" charges that acts as sources and sinks for dipole charge.
We have considered the entanglement spreading in fractonic circuits. Since the initial conditions we consider are close to maximal entanglement according to the conventional entanglement entropy measure, we have to use different entanglement measures. We have considered three different entanglement measures: observable entanglement entropy, operator entanglement, and entanglement spectrum. By looking at observable entanglement entropy, we have found that while observable entanglement grows ballistically, the entanglement velocity is much smaller than the front propagation velocity. Additionally, a system initialized with fracton charge appears to saturate to an "area law" observable entanglement entropy independent of system size, consistent with the lack of thermalization. Meanwhile, the operator entanglement is volume law but subthermal, whereas the entanglement spectrum exhibits level statistics that follow a semi-Poisson distribution, similar to many body localized eigenstates [73] but distinct to either thermalizing or conventionally integrable systems. Finally, we have considered initial conditions that contain multiple fractons. If the initial condition contains multiple fractons, then the system evolves to a localized configuration where all fractons agglomerate at their "center of charge", and the rate of evolution to this configuration is determined by the initial separation of the fractonic charges. We have conjectured that this fracton agglomeration may be a consequence of feedback from the nonconserved operator to the conserved operator sector. Regardless of the origin of this attraction, however, the long time steady state attained is non-ergodic, even at nonzero fracton density i.e. it is not equivalent to a thermal state (defined as an equal weight superposition of all S z product states with a given charge and dipole moment). Such a "thermal" state does not exhibit charge localization, and has markedly different entanglement properties. This constitutes the first example that we know of in which unitary time evolution with no time translation symmetry at all leads to localization.
We comment now on the implications of our results. Firstly, our work conclusively establishes that random circuit dynamics can lead to non-ergodic dynamical phases. This has obvious implications for quantum information processing. Particularly important is our prediction that localization in quantum circuits can persist to two dimensions, since actual information processing architectures cannot be purely one dimensional, such that our work opens up the possibility of harnessing localization for future scalable quantum circuits. Insofar as random circuit dynamics is a minimally structured model for time evolution, we expect the results to be robust, and applicable to any generic model with U (1) charge conservation and dipole conservation, and no other constraints.
Our work also has intriguing implications for the Hamiltonian dynamics of fracton systems. In particular, since Hamiltonian dynamics is more constrained than random circuit dynamics (because Hamiltonian dynamics must conserve energy), the immobility of fractons under random circuit dynamics in d = 1, 2 strongly suggests that fractonic systems in d = 1, 2 should also fail to thermalize under Hamiltonian dynamics. That is, fractonic systems in d = 1, 2, prepared in a sector with non-zero fracton charge but zero fracton charge density should have the fractons be strictly immobile, preserving forever a memory of their initial conditions in local observables -the very definition of MBL [3]. Moreover, this localization should survive to non-zero energy densi- ties (unlike [54]), since while initializing the system with a non-zero density of dipoles will "break" the locator expansion style arguments of [54], it will not affect the mechanism underlying the immobility of fractons in our analysis. The localization of fractonic charges will also be robust to arbitrary local perturbations -after all, it is robust to application of a random circuit. In other words, low dimensional fracton systems can realize MBL even in the absence of disorder in the Hamiltonian, and even in d = 2. The obstructions to MBL in higher dimensions and in translation invariant Hamiltonians [64,65] do not apply, since these are essentially obstructions to construction of a convergent locator expansion, and our argument for fracton localization does not rely on locator expansions. By thus liberating the study of MBL from reliance on the crutch of locator expansions, our work may also open a new chapter in this field. Moreover, the localization in low dimensional fractonic systems should be even more robust than traditional MBL, since the localization mechanism we have discovered herein is insensitive to noise, whereas conventional MBL does not survive noise [76]. 8 Our work also opens up new possibilities for noisy fracton dynamics in three spatial dimensions (where fracton physics is best understood). Previous explorations of fracton dynamics in three dimensions (e.g. [55,56]) have essentially demonstrated the lack of convergence of "locator expansion" type approaches to fracton dynamics at non-zero energy densities. This lack of convergence has been interpreted in terms of thermalization, and the diffusion constant inferred from the breakdown of the locator expansion. However, in the present work we have identified a new mechanism for ergodicity breaking, which does not rely on locator expansions. This re-opens the possibility of a non-ergodic long time limit for three dimensional fracton systems. Prima facie, our analysis does not apply to three dimensional systems, since our analysis relies on the return probabilities of low dimensional random walks, and in three dimensions random walks need not return to the origin. However, the "dipolar" excitations in fractonic models can frequently be "subdimensional" [31], being restricted to move in a space of lower dimensionality. The localization of fractonic charge is tied to the return probability of dipoles, and the subdimensional character of dipoles in certain three dimensional fracton models may open the door to non-ergodic long time states for three dimensional fractons. This too is an interesting topic for future exploration.
Finally, it is worth critically examining the role of quantum mechanics in the localization of fractons. Throughout, we have made heavy use of the machinery of quantum mechanics, such as quantum states and operators. Indeed, our results have several intrinsically quantum-mechanical features, such as operator spreading dynamics and unusual entanglement behavior. We have also invoked feedback between the conserved and non-conserved operator sectors, which does not have an obvious classical analogue. However, our primary argument for localization, in terms of the guaranteed return of low-dimensional random walks, may well hold at the classical level as well. It is therefore plausible that there may be an analogous classical phase featuring localized charge dynamics. This would open the door to a new type of purely classical localization phenomenon, in contrast with the quantum mechanical origin of standard MBL phases. We leave this question as a topic for future investigation.

ACKNOWLEDGMENTS
The authors acknowledge useful discussions with Yang-Zhi Chou, David Huse, Vedika Khemani, Abhinav Prem, Jonathan Ruhman, and Sagar Vijay. We thank Paul Romatschke and Ted Nzuonkwelle for providing access to the Eridanus cluster at CU Boulder. This material is based upon work supported by the Air Force Office of Scientific Research under award number FA9550-17-1-0183 (SP, MP and RMN); by the U.S. Department of Energy, Office of Science, Basic Energy Sciences (BES) under Award number de-sc0014415 (SP); a Simons Investigator Award to Leo Radzihovsky (MP); by NSF Grant 1734006 (MP), and by the Foundational Questions Institute (fqxi.org; grant no. FQXi-RFP-1617) through their fund at the Silicon Valley Community Foundation (RN and MP). We also acknowledge the support of the Alfred P. Sloan foundation through a Sloan Research Fellowship (RN). We include numerical results for S = 10 and S = 100 spin chains subject to random unitary evolution via 3qudit gates. Figure 24 shows the fracton peak in the right-weight profiles for S = 10 and S = 100. In Figure  25, we show the power-law behavior of the tails lagging behind the front and provide numerical evidence for the −5/2 exponent appearing even at large S. And finally, we plot the integrated weight under the fracton peak as a function of system size, and see that it saturates to a finite value in the thermodynamic limit in both cases (see Figure 26). Altogether we conclude that our basic results survive as we take the limit of large S.