Operator spreading and the emergence of dissipative hydrodynamics under unitary evolution with conservation laws

We study the scrambling of local quantum information in chaotic many-body systems in the presence of a locally conserved quantity like charge or energy that moves diffusively. The interplay between conservation laws and scrambling sheds light on the mechanism by which unitary quantum dynamics, which is reversible, gives rise to diffusive hydrodynamics, which is a dissipative process. We obtain our results in a random quantum circuit model that is constrained to have a conservation law. We find that a generic spreading operator consists of two parts: (i) a conserved part which comprises the weight of the spreading operator on the local conserved densities, whose dynamics is described by diffusive charge spreading. This conserved part also acts as a source that steadily emits a flux of (ii) non-conserved operators. This emission leads to dissipation in the operator hydrodynamics, with the dissipative process being the conversion of operator weight from local conserved operators to nonconserved, at a rate set by the local diffusion current. The emitted nonconserved parts then spread ballistically at a butterfly speed, thus becoming highly nonlocal and hence essentially non-observable, thereby acting as the"reservoir"that facilitates the dissipation. In addition, we find that the nonconserved component develops a power law tail behind its leading ballistic front due to the slow dynamics of the conserved components. This implies that the out-of-time-order commutator (OTOC) between two initially separated operators grows sharply upon the arrival of the ballistic front but, in contrast to systems with no conservation laws, it develops a diffusive tail and approaches its asymptotic late-time value only as a power of time instead of exponentially. We also derive these results within an effective hydrodynamic description which contains multiple coupled diffusion equations.


I. INTRODUCTION
The nature of quantum dynamics and thermalization in isolated many-body systems is a topic of fundamental interest. Over the last few years, a remarkable confluence of theoretical progress and experimental advances in engineering and controlling isolated many-body quantum systems, especially in cold-atomic gases, has led us to re-examine our understanding of the very foundations of quantum statistical mechanics 1-3 . On the theory side, research in the field of many-body localization (MBL) [4][5][6][7][8][9] has revealed the existence of classes of generically interacting systems that do not obey quantum statistical mechanics, and understanding the nature of different MBL phases 10,11 and the transition(s) between MBL and thermalizing phases 6,[12][13][14][15][16][17][18] is an active area of research. Complementarily, the development of tools like the AdS/CFT duality 19,20 has led to new perspectives on the dynamics of thermalizing strongly-interacting systems. This duality has been used to relate the physics of information scrambling in black-holes to the process of thermalization in condensed-matter systems of interacting spins and/or particles [21][22][23][24][25][26][27][28][29][30] .
One lens on the dynamics of isolated many-body quantum systems is provided by studying the spreading of initially-local operators under the system's unitary timeevolution. In the Heisenberg picture, an initially local operator O 0 evolves into O 0 (t) = U † (t)O 0 U (t) with support on a spatial region that grows with time. This spreading of O 0 (t) is reflected in the growth of the commutator between O 0 (t) and a typical local operator W x with support near position x. If x is well away from the origin, then W x initially commutes with O 0 . We define the "OTOC" as the norm of the "out-of-time-order" commutator 31 where the expectation value is taken in an appropriate equilibrium ensemble. The OTOC is related to the commutator norm that appears in Lieb-Robinson bounds 32 , and has received a great deal of attention recently as a diagnostic of information "scrambling" in quantum chaotic systems [21][22][23][24][25][26][27][28][29][30][33][34][35] . For a spin-1/2 chain of length L, a complete orthonormal basis for all operators is given by the 4 L "Pauli strings" S, which are products of Pauli matrices on distinct sites. We can then express our spreading operator in this basis of Pauli strings: The initially local operator O 0 consists only of strings S that are the local identity at all sites except one or a few sites near the origin. But, with time, the strings that dominate this sum grow in spatial extent, containing non-identity local operators at sites out to a "front" at a distance from the origin that grows with time. The OTOC remains near zero as long as the strings that dominate in O 0 (t) contain only local identities near position x, but it becomes nonzero once the operator front reaches this position.
In this paper, we are interested in understanding operator spreading and scrambling in the physically ubiquitous setting of thermalizing systems with one or a few local conservation laws that result in diffusive transport. This includes, for example, certain Hamiltonian models, since they do conserve energy, and in many cases energy transport is diffusive. It also includes certain Floquet systems which don't conserve energy but do conserve a charge or spin. We first investigate this question in a random tensor network spin chain model where the tensors are constrained to obey a single local U (1) spin conservation law. We can derive a number of exact results in this setting, which we conjecture (and numerically verify) should also universally apply to Hamiltonian and Floquet spin chains with diffusing conserved quantities.
A set of recent papers has studied entanglement and operator dynamics in random tensor network models with no conservation laws 33,[36][37][38] . A subset of these 37,38 showed that, despite the absence of traditional hydrodynamic conservation laws, unitarity alone amounts to a type of conservation law since the operator norm, , is conserved which means that the total weight of the operator on all Pauli strings, S |a S | 2 , is conserved. This leads to an emergent "hydrodynamical" picture for describing operator spreading wherein, in one dimension, the dynamics of the operator front can be described by a distribution of biased random walkers (where the bias reflects the fact that it is more likely for an operator string to grow rather than to shrink) 37,38 , while in higher dimensions the front is modeled by a random growth model 36,37 . This leads to a picture in which the operator front propagates ballistically with "butterfly" speed v B , while in low dimensions the width of the front grows as a sublinear power of time (∼ √ t in 1D) 37,38 . This means that the bulk of the total weight of O 0 (t) is contained in the spatial region lying within the "light cone" defined by the front whose spatial linear size grows linearly with time. Further, O 0 (t) is scrambled within the spatial region defined by the light cone and, with respect to some (but not all) measures, resembles an unconstrained random operator in this region 39 .
In this work, we build on the picture above and study the interplay between the "actual" hydrodynamics governing the dissipative diffusion of conserved charges and the "emergent" hydrodynamics of unitary operator spreading in systems with conservation laws. Starting from a local operator that contains the conserved charge, we find that there is again a ballistically propagating front describing the operator spreading. However, unlike the fully nonconserved case, O 0 (t) has a significant (decreasing as a power law in time) amount of weight on operator strings whose fronts lag far behind the main operator front. This can be understood as follows: Since conserved charges spread diffusively, the parts of O 0 (t) that overlap with the conserved charges are "left behind" in a region of linear size ∼ √ t near the origin, even while the main operator front has reached distance v B t. But the total weight of O 0 (t) on these conserved operators decreases as a power-law in time ∼ t −d/2 in d dimensions as the conserved density spreads diffusively. This loss of operator weight makes this diffusion effectively nonunitary and thus dissipative. But the full dynamics of the system is unitary so the full operator weight is not lost. Instead the dissipation due to the diffusive currents of the conserved density steadily converts operator weight from conserved to non-conserved operators, thus emitting a "flux" of non-conserved operators that then spread ballistically. This emitted flux is proportional to the square of the diffusive current, as expected for "Ohmic" dissipation. One point of view is that once the non-conserved operators start spreading ballistically, they rapidly become so highly nonlocal that they stop being observables. Then they can be viewed as effectively random, functioning as the bath whose entropy increases due to the dissipative diffusion.
Moreover, in the models we study here, we can also study the unitary "inner workings" of this effective bath, following the spreading of the non-conserved parts of the operator. The bulk of these non-conserved operators are emitted almost immediately, near time t = 0, and these grow to form the leading operator front at distance v B t at time t. However, as a result of the conservation law, there is still significant weight left on the conserved parts even after the initial emission, and these continue to emit non-conserved operators at a slow rate that decreases only as a power law in time. Thus, those non-conserved parts of the operator that are emitted at a later time t e have a front that lags behind the main operator front by distance v B t e . This leads to a power-law tail in the operator profile behind the front. By contrast, the unconstrained random circuit model does not show such a power-law tail, since in that case there are no slow dissipative modes so no significant part of the operator has a front that lags substantially behind the main operator front. Fig. 2 shows a sketch of the operator profile depicting all three regimes: (i) the diffusive conserved charges remaining near the origin, (ii) the ballistically moving front, and (iii) the power law tail behind the front.
The picture that emerges from our study is of multiple coupled hydrodynamic equations: The first is the dissipative diffusion of the conserved quantity. This dissipation serves as the source of the ballistically spreading nonconserved operators. The dynamics of the fronts of these nonconserved operators is biased diffusion in 1D and a random growth model in higher dimensions 36,37 . Moreover, we find that the local operator content in the interior of the spreading operator is also governed by two coupled noisy diffusion equations, with the leading front of the operator serving as a moving boundary condition on these equations. And finally, there is at least one more "layer" of this hydrodynamics that governs the entanglement dynamics of the spreading operator 40 .
The spatial profile of the spreading operator described above is also reflected in the behavior of the OTOC C(x, t) defined in Eq (1). We show that C(x, t) increases sharply when the ballistic operator front reaches x but, as a result of the power-law tails behind the main operator front, it only approaches its asymptotic late-time value as a power-law in time. This is contrast to systems without conservation laws where there are no such power-law diffusive tails 37,38 . Our results help explain the numerical observations in Ref. 41 where this late-time power-law in the OTOC was observed in systems with conservation laws.
We note that much recent work has focused on computing the OTOC at low temperatures in systems with energy conservation, and bounds on the growth of the OTOC have been derived in this setting 42 . On the other hand, random circuit models do not conserve energy and thus the "infinite temperature" Gibbs ensemble is the only meaningful one for such models. Nevertheless, for circuit models endowed with extra conservation laws like total spin/charge, the dynamics can be resolved into different spin sectors with the "chemical potential" µ now playing a role somewhat analogous to the inverse temperature β. We mostly focus on µ = 0 in this work, while briefly addressing the µ = 0 case.
We note that several papers have recently noted that quantum information can spread ballistically in systems with diffusively relaxing conserved charges [43][44][45][46][47][48] . This, by itself, is not that surprising since even MBL systems with no charge transport can show logarithmic spreading of entanglement 8,49 . Of these papers, Refs. 46-48 study a weakly interacting diffusive metal and, using a perturbative semiclassical scattering calculation, relate the butterfly speed v B to the diffusion constant of the metal via a "Lyapunov exponent" which characterizes the exponential growth of "chaos" in this semiclassical setting as measured by the growth of the OTOC before the front arrives: On the other hand, Refs. 43-45 numerically study fully quantum spin-chains which cannot be treated semiclassically, and for which no prolonged period of exponential growth in the OTOC has been observed to date (the existence of a fully quantum Lyapunov exponent in spatially extended systems with small local Hilbert spaces and only short-range interactions, while perhaps expected, is a presently unresolved question). In this fully quantum setting where semiclassical analytical methods don't apply, the aforementioned numerical papers [43][44][45] generally treat the diffusive charge relaxation and the ballistic information spreading as two independent numerical observations that do not interface with one another.
One of our main contributions in this work is to connect the diffusive charge dynamics with the ballistic front spreading in strongly quantum systems with local conservation laws into a composite picture for the operator profile, showing how the different regimes connect at different time and length scales. En route, we elucidate how reversible unitary dynamics can still display dissipative hydrodynamic modes, wherein the dissipative process is the conversion of operator weight from locally observable conserved parts to non-local and essentially unobservable non-conserved parts at a slow rate set by the local diffusion current of the conserved densities. Thus while the von Neumann entropy of the full system is conserved under the unitary dynamics, we show how the local "observ-able" entropy can still increase at a slow hydrodynamic rate, very concretely illustrating a resolution of the fundamental tension between unitarity and dissipation in closed quantum systems.
The balance of this paper is structured as follows. We begin in Section II with a description of our constrained random unitary circuit model which has a local U (1) conservation law total for total spin. This is followed by a detailed discussion of the spatial profile of spreading operators in Section III. We show that the conserved charges evolve diffusively under the action of the circuit. The coupling between the charge conservation and unitarity produces a steady flux of operator weight from the diffusively spreading conserved components to ballistically spreading nonconserved components, leading to a power law tail in the spatial profile of the operator weight. We present the coupled hydrodynamics describing this process. Next, in Section IV, we expose another layer of structure in the spreading operator by studying the local operator content of the highly nonlocal operator strings within the light cone defined by the leading ballistic front. We find that the distribution of different local operators is itself governed by a set of coupled noisy diffusion equations that "turn on" once the front passes a given position. In Section V, we turn to a discussion of OTOC's in this model, and find that the diffusive processes governing the shape and internal structure of the spreading operators lead to a late time diffusive tail in OTOC's involving the conserved charges. We numerically verify that the universal aspects of our results also apply to more "physical" spin chains with energy/charge conservation in Section VI, and conclude in Section VII. Some additional details are discussed in the Appendices.

II. RANDOM UNITARY CIRCUIT MODEL
As an explicit example where we can obtain analytical results, we consider a one-dimensional chain of length L sites where the degree of freedom on each site is the direct product of a spin-1/2 (or qubit) and a qudit with Hilbert space dimension q. For definiteness we take a finite L, but we will be interested in the behavior in the limit of infinite L. The time-evolution is constrained to conserve the total z component of the spin-1/2's, which we call S tot z , while the qudits are not subject to any conservation laws. We note that we include the qudit at each site because some results can only be obtained analytically in the large-q limit, although some of our results do apply for all q, including the case q = 1 that has no qudits.
Generalizing from Refs. 33, 36-38, the time-evolution is governed by a random quantum circuit comprising staggered layers of two-site unitary gates acting on even and odd spatial bonds at even and odd times respectively (see Fig. 1). The time evolution operator is given FIG. 1. Left: a diagram of the random unitary circuit. Each site (black dot) is the direct product of a two-state qubit and a q-state qudit. Each gate (blue box) locally conserves S tot z , the total z component of the two qubits it acts upon, and is thus a block-diagonal unitary of the form shown on the right, with each block of each gate independently Haarrandom. The smaller blocks do not flip the qubits and thus operate only on the two qudits, while the larger block also produces S tot z -conserving qubit "flip-flops". by As a result of the conservation law, each two-site unitary Labeling the spin state on each site i as (↑ a) i or (↓ b) i , where the first label is the spin state in the Pauli z basis and the second label is the qudit state, the structure of U i,i+1 looks like: (i) a (q 2 × q 2 ) block acting in the Each of these blocks is a Haar-random unitary, and each block in each two-site gate is chosen independently of all others. To characterize the time-evolution of local operators, it is useful to define a complete orthonormal basis of operators on each site. For the spin, we can use the Pauli matrices on each site to define an onsite basis as so r i and l i are suitably normalized spin raising and lowering operators, respectively. These basis operators all have a definite ∆S tot z (such as "raise/lower by one") under the U (1) symmetry that conserves S tot z and are thus more convenient than the Pauli σ x/y i matrices for characterizing the U (1)-conserving dynamics. For the qudit, one can construct higher-dimensional generalizations of the Pauli matrices {Σ µ=0,1,···q 2 −1 i } that are normalized such that Tr(Σ µ † i Σ ν i )/q = δ µν . Then, the tensor product B µν i ≡ σ µ i ⊗ Σ ν i generates a local basis for the 4q 2 operators acting on each composite site i, denoted in shorthand as (IΣ ν ) i , (rΣ ν ) i , (lΣ ν ) i and (zΣ ν ) i . Using this basis, the time evolved operator O(t) can be expanded as: where each generalized Pauli string S is one of (4q 2 ) L basis product operators, i B µiνi i . Since the basis strings satisfy Tr[S † S ]/(2q) L = δ SS , the coefficients a S can be obtained as a S (t) = Tr[S † O(t)]/(2q) L . Finally, we normalize the initial operator O 0 such that Tr[O † 0 O 0 ] = (2q) L which, by the unitarity of the dynamics, implies that the total weight of O(t) on all strings S is also normalized to 1: This sum rule is the effective conservation law due to unitarity 37,38 .
There are a few classes of operators on site i that evolve differently under the action of this conserving unitary circuit. First, (zI) i measures the local conserved charge, and (II) i is the identity operator. The conservation law implies that S tot z is conserved so that and the operators ( 2 are left invariant by the action of all local gates U i,i+1 . Further, if one starts with an operator with a definite ∆S tot z under the U(1) symmetry (for example, (rΣ ν ) i raises the spin by one), the action of the circuit preserves this ∆S tot z . Appendix A summarizes the action of U i,i+1 on all possible two-site operators.
It will be subsequently useful to separate the spreading operator into conserved and non-conserved pieces. To intuitively understand this separation, consider an initial density matrix where I all is the background equilibrium state which is the identity on the full system; AO 0 is a traceless local onsite perturbation at the origin to this equilibrium state, with O 0 normalized such that Tr[O † 0 O 0 ] = (2q) L , and A is the amplitude of this perturbation, which must be small enough so ρ remains non-negative. The system conserves S tot z (6) so that where (t) denotes expectation values in the state ρ(t).
If the perturbation injects some local charge at the origin then, on general grounds, we expect this "extra" charge to spread diffusively so that (zI) x ∼ 1 √ t e −x 2 /4Dt . We will see how this diffusion arises in operator language and explore its consequences for operator dynamics in this system. To separate the part of O 0 (t) that has overlap with the conserved charges, we define a c i (t) as the amplitude of the conserved charge (zI) i in the operator expansion of O 0 (t): The conserved charges which act as (zI) i on site i and as the identity everywhere else are a subset of the full basis of operator strings, and the "conserved part" of O 0 (t) is defined as the part of O 0 (t) with weight on these basis strings: with the "non-conserved part" being the rest, O nc The dynamics of operator spreading in this model is governed by the interplay between (i) this explicit charge conservation (11) which is a sum rule on the amplitudes of the conserved operators, and (ii) the conservation of the operator weight (which follows from unitarity) which is a sum rule on the squares of the amplitudes of all basis strings (5).

III. "SHAPE" OF SPREADING OPERATORS
A complete characterization of the spreading and scrambling of an initially local operator O 0 requires knowledge of the exponentially-many coefficients a S (t) in the expansion of O 0 (t) as in Eq. (4). Instead of doing this, we first consider a more coarse-grained approach and define the "right-weight" ρ R (i, t) as the total weight in O(t) of basis strings that end at site i, which means that they act as the identity on all sites to the right of site i, but act as a non-identity on site i. This is the weight on all strings of the form The conservation law on ρ R (i, t) which follows from unitarity (5) gives ρ R (i, t) the interpretation of an emergent local conserved "density", and Refs. 37,38 showed that the (hydro)dynamics for ρ R (i, t) is governed by a biased diffusion equation. Of course, one can analogously define the left-weight ρ L (i, t) and, together, these can be used to characterize some illuminating aspects of the spatial structure of the spreading operator (we will consider other measures probing the local operator content inside the spreading operator in the next section).
It is instructive to first understand the dynamics of ρ R (i, t) in an unconstrained random circuit model in 1D with local Hilbert space dimension q. It was shown in Refs. 37,38 that the weighted distribution of the endpoints of the basis strings in the operator time-evolve as biased random walkers, where the bias reflects the fact that it more likely for the strings to grow rather than to shrink. If we define the "right-front" of a string as the location of the rightmost unitary gate that sees a nonidentity, one can obtain the probabilities for the front to move forwards or backwards by noting that under the action of the front gate, all the (q 4 − 1) non-identity operators at the gate are produced with equal weights, on average. If the front gate is U i,i+1 , only q 2 − 1 of the operators that can be produced by this gate act as the identity on the right site of the gate (i + 1), and each of these outcomes results in the front moving a step backwards, an event with probability p = 1 (q 2 +1) . Then, from the theory of random walks, the mean location of the front after time t is which defines the butterfly speed v B , and the width of the front grows as ∼ )/2 is the diffusivity of the resulting biased random walk. Evidently, the front location can be described by the emergent random walk hydrodynamics 37,38 , and the fact that the right side of this equation is a total spatial derivative reflects the conservation of the emergent density ρ R (12) .  Fig 2(c) shows a sketch of the Haar-averaged front for this random circuit evolution which, in the scaling limit t, x → ∞ for x ≈ v B t takes the form 37,38 In the limit q → ∞, the front becomes sharp and v B ∼ 1 − 2 q 2 → 1 so that the front deterministically moves forwards for all basis strings at each time step in this limit (Fig 2(d)). Note that the geometry of the circuit limits the operator growth to at most one site on the left and right ends per time step. This imposes a strict upper bound on v B set by the "causal limit speed" v CL = 1, and v B goes to this limit as q → ∞.
Let us now turn to how this picture changes in the presence of an explicit U (1) conservation law with diffusively spreading charges.

A. Spreading of the local conserved charge
Let us begin with the case when the initial operator is a conserved U (1) charge located at the origin, O 0 = (zI) 0 . In this case, the spatial structure of the operator shows three regimes: (i) a ballistic front, (ii) a power-law tail behind the front, and (iii) diffusively spreading charges near the origin.

Diffusion of conserved charge
We start by discussing the diffusive dynamics of the conserved parts of the operator. For the initial operator O 0 = (zI) 0 , the initial conserved amplitudes are a c i (0) = δ i0 and the conservation laws (8), (11) give When viewed as a perturbation to an equilibrium state (7), O 0 = (zI) 0 creates an excess of charge at the origin which should spread diffusively. Thus, on general grounds, at late times we expect (zI) where D c is the charge diffusivity. We now see how this arises.
To understand the evolution of a c i (t), consider the action of a particular gate U 12 on the superposition [a c 1 (zI) 1 + a c 2 (zI) 2 ]. It can be shown (Appendix C) that the action of a gate makes the average amplitudes a c equal on the two sites acted upon by it, while preserving the sum of amplitudes. That is, after the action of the gate, the Haar-averaged amplitudes are and the gate has produced a current between these two sites of the conserved charge that is ∼ (a 1 (t) − a 2 (t)) ∼ ∂ x a(x, t). Note the "smoothing" action of the circuit which locally makes the averaged amplitudes equal, and thus reduces ∂ x a(x, t) in time. This action gives a binomial charge distribution (Appendix C): which is an exact result true for all q including q = 1. If we coarse-grain, in the scaling limit x, t → ∞ we get showing diffusion of the conserved charge with diffusion constant D c = 1/2. If we consider similar random circuits acting on a system in higher dimensions d > 1, a similar diffusive behavior will be present, just with diffusion along all directions. It is noteworthy that as a consequence of the conservation law, the Haar-averaged amplitudes a c i (t) are non-zero (16). In an unconstrained random circuit, only the squared weights |a S | 2 survive Haar averaging while all amplitudes average to zero. Likewise, in our model with S tot z conservation, certain off-diagonal products of amplitudes a S a S can also survive Haar averaging (Appendix C), while all such averages are zero in the unconstrained model.

Ballistic front and power law tails
We now turn to the effect of the diffusive dynamics on the shape of the spreading operator as measured by ρ R (i, t). It is important to note that ρ R (i) measures the weights, |a S | 2 , of all operator strings ending on i, while the preceding discussion was about the diffusive dynamics of the amplitudes, a c i , of only the conserved charges in the expansion of O 0 (t); these conserved operators are strings of length one site. Of course, the weight of O 0 (t) on a conserved charge at site i contributes to ρ R (i, t) and it is convenient to separate out this contribution: where ρ nc R denotes the right-weight from all "nonconserved" operator strings that are not one of the con- We will show that the total weight of O 0 (t) on the conserved charges, ρ c tot , decreases as a power law in time; these then act a source, continuously emitting a flux of non-conserved operators that spread ballistically and contribute to ρ nc tot . 50 Let us see how this comes about. First, we show in Appendix C that the circuit-to-circuit variance in the conserved amplitudes is suppressed both in the large q and the late time t limit: while the leading term |a c i (t)| 2 scales as ∼ 1/t (17). In this limit, ∆ a i (t) ≈ 0 and |a c i (t)| 2 ≈ |a c i (t)| 2 . Then, the power law decrease in ρ c tot is obtained as where we coarse-grain at long times to obtain the penultimate equality, using (17). Likewise, in higher dimensions, this power law decrease will scale as ρ c tot (t) ∼ t −d/2 as can be easily seen by considering the higher dimensional generalization of (17). Because of the conservation of total density (19), this decrease in ρ c tot has to be compensated by a corresponding increase in ρ nc tot . Indeed, from the previous discussion, we see that the local decrease in ρ c This is the increase in ρ nc tot that is locally generated at this gate, and this quantity is proportional to the square of the conserved quantity's local current. These showing the spreading of an initially local conserved charge (zI)0(t) in a random circuit model with S tot z conservation in a system of size L = 1000 at different times t. These profiles depict three regimes: (i) a "lump" in the region |x| √ Dct reflecting the weight of the operator on diffusively spreading conserved charges (shaded purple). This lump emits ballistically spreading nonconserved operators at a slow power-law rate. This emission creates (ii) the leading ballistic "fronts" near |x| ∼ vBt within which the majority of the operator right-and left-weight is contained (shaded red for the latest time). These leading fronts are from nonconserved operators emitted at early times and they are perfectly sharp at q = ∞ where vB = 1 (b), and have a width Dρt for finite q (a); Finally, the slow emission also leads to (iii) diffusive tails ∼ (vBt − |x|) −3/2 behind the leading fronts which reflect the operator weight in "lagging" fronts of nonconserved operator strings that were emitted at later times (shaded blue for the latest time). The curves in (a) are obtained via a simulation at q = 3 which takes into account the different processes (diffusion of charges, emission of nonconserved operators and the biased diffusion of the nonconserved right-and left-weights) to order 1/q 2 . The red dashed curve is the exact infinite q answer for the "tail" (24). (c,d): For comparision, ρ R/L (x, t) in an unconstrained random circuit model 37,38 where z0(t) isn't "special". Regimes (i) and (iii) do not exist in an unconstrained circuit, and the ballistically spreading operator fronts describe the entire right-and left-weight profiles. The fronts are again infinitely sharp at q = ∞ (d) and have a finite width ∼ Dρt for q < ∞ (c).
newly-produced non-conserved operators then evolve under the unitary dynamics and get converted to other nonconserved operators with increasing size. As in the case of the evolution of ρ R (x, t) for an unconstrained random circuit, the front of non-conserved operators (once generated) spreads ballistically, with v B = 1 for q = ∞ since the likelihood of the front moving backwards is again suppressed by ∼ 1/q 2 (Appendix B). Further, in the q = ∞ limit, there is no "backflow" of density from non-conserved operators to conserved charges; this backflow only appears at order 1/q 4 . Thus, at q = ∞ and in the scaling limit, these considerations imply that: This expression tells a very natural story. The total nonconserved right-weight at position x at time t is the integrated weight of all "fronts" emitted at locations y at times t y = t − (x − y)/v B such that the front travelling with velocity v B makes it to position x at time t. Moreover, since the conserved charges are primarily spread within a distance ∼ √ D c t near the origin, the emission of non-conserved flux is only significant at locations within this diffusively spreading "lump of charge". Then, at late time t, we see the three pieces in the shape of ρ R (x, t) mentioned earlier: 1. The diffusive "lump": In the spatial region |x| √ D c t near the origin, the right-weight comes almost entirely from the diffusively spreading conserved part of the operator, ρ R (x, t) |a c (x, t)| 2 . As a result, the spreading operator has significant weight that is "left behind" near the starting position of the operator, and this weight decreases only as a power-law in time ∼ t −d/2 in d dimensions (21). By contrast, a spreading operator in the unconstrained circuit model has negligible right-weight (exponentially decreasing with t) near its initial location at late times.
2. The ballistic front: The leading ballistic operator front (at the right end) is at x = v B t, and the weight at the leading front is from non-conserved operators that were emitted at early times. At q = ∞, the leading front is sharp, the right-weight is strictly zero for x > v B t with v B = 1, and the sharp front is due to those non-conserved operators that were emitted by gates acting at the precise edge of the causal light cone. At finite q, as in the unconstrained circuit model 37,38 , the front distributions execute biased diffusion instead of strictly moving forwards at each time step. This leads to an order 1/q 2 correction in v B : v B 1 − 8/(9q 2 ) (Appendix B) and gives the main operator front a nonzero width ∼ D ρ t ∼ t/q 2 . Thus, at finite q, the leading front is mostly due to nonconserved operators that were emitted at early times t e D ρ t/v B . For systems in d > 1, the broadening of the front at finite q is given by a random growth model and grows (if at all) with a smaller power of time 36,37 .
3. The diffusive tail: Finally, the fronts of nonconserved operator strings that were emitted at later times t e D ρ t/v B lag behind the leading front by v B t e , leading to the development of powerlaw tails in the right-weight behind the main operator front. Consider positions x well separated from both the leading ballistic front and the diffusively spreading charge "lump" such that which has a power-law tail in both space and time. Indeed, for x far separated from the lump, the diffusive charge lump emitting non-conserved flux can be approximately treated as a point source with the same integrated weight, so that ρ nc R (x, t) is given by the rate of change of the total conserved density in the lump ∂ρ c tot /∂t ∼ t −3/2 at time t e = (t − x/v B ), which explains the 3/2 power (21). The leading functional dependence of this tail is the same at both infinite and finite q. For d > 1 the exponent in this power law is 1 + (d/2). By contrast, there is no such tail in the spatial profile of operators evolving under unconstrained circuits, since such circuits have no mechanism for generating fronts that significantly lag behind the main operator front. At q = ∞ for the unconstrained circuit, ρ rand so there is strictly no tail; at finite q, the right-weight behind the leading front for the unconstrained circuit falls off exponentially in time at locations x behind the front scaling with The left-weight shows the same three regimes, by reflection symmetry. Fig. 2 shows ρ R/L (x, t) for both the un-constrained and charge-conserving random circuit models at infinite and finite q.

Hydrodynamic description
We now turn to a long-time hydrodynamic description of the coupled processes involving diffusion of the conserved "charge" and the propagation of the fronts of the nonconserved operators emitted from this diffusing conserved charge.
For specificity, we restrict our attention to systems that are statistically translationally invariant and inversion symmetric. In such systems the amplitudes of the conserved part of the operator obey an unbiased diffusion equation. In more than one dimension, the diffusivity may not be isotropic; if that is the case, we rescale distances along the eigendirections of the diffusivity to make it isotropic: where D c = 1/2 in our d = 1 random circuit model (17). This diffusion is also subject to noise which will lead to fluctuations in a c across different realizations of the random circuit but, as we discussed above (20), for initial operators that do contain the conserved quantity the noise is a parametrically subleading correction to the amplitudes a c in the late time limit that is relevant to the hydrodynamics. Diffusion is a dissipative process, and this is reflected here in the decrease with time of the total operator weight of the conserved part of the operator, ρ c tot (t) = dx|a c (x, t)| 2 . Since the full system is undergoing unitary dynamics, this loss of conserved operator weight means that weight is being converted to nonconserved operator weight. The density of local rate of this emission of nonconserved operators is 2D c |∇a c | 2 (22), where we show that this coefficient is 2D c below. Note that this dissipation is proportional to the conserved quantity's local current squared, as in Ohm's law. These "emitted" nonconserved operators then spread rapidly, becoming highly nonlocal. This is an example of a presumably much more general picture of how dissipative processes happen within closed systems undergoing unitary dynamics: Correlations captured by low-order observables are moved by the unitary dynamics to highly nonlocal, and thus effectively non-observable operators. Dissipation in such closed systems is thus the "hiding" of correlations in highly nonlocal operators so that the correlations that remain detectable to low-order observables are reduced and thus the "observable" entropy increases, even though the von Neumann entropy of the full system remains unchanged.
The dynamics of the spreading of the front of the emitted nonconserved operators is given by a random and nonlinear growth model for d > 1, which we will not discuss further here 36,37 . For d = 1 the weighted distributions of the fronts of the spreading strings move by biased diffusion 37,38 . We will focus on the right-moving front, but the left-moving front is doing the same thing, just spatially reflected so it moves in the opposite direction.
Thus the leading order continuum hydrodynamics in d = 1 of the right-density ρ nc R (x, t) is described by a diffusion equation with drift and a source term representing the emission of non-conserved operators from local gradients in a c (x, t): where we show in Appendix B that the drift and diffusion constants are v B 1 − 8 9q 2 and D ρ 8 9q 2 neglecting corrections of order 1/q 4 and higher. Further, it was shown in Refs. 37,38 that the circuit-to-circuit fluctuations in ρ R (x, t) (and hence ρ nc R (x, t)) scale with a parametrically smaller power of t, so we can ignore the noise in ρ nc R in the leading order hydrodynamics. Thus, the coupled diffusion equations (25) and (26) describe the leading order hydrodynamics for those aspects of the shape of the operator that are described by ρ R (x, t).
It is instructive to also directly consider the hydrodynamics for the total right-weight, ρ R (x, t) = ρ nc R (x, t) + (a c (x, t)) 2 , an exercise that explicitly reveals the form of the source term for ρ nc where we've left the functional form of the source term S nc for ρ nc R undetermined, and we've used (25) in evaluating the time-derivative for the squared amplitudes of the conserved charges. Note however that since dxρ R (x, t) is conserved, the continuity equation requires that ∂ t ρ R (x, t) is the total spatial gradient of a current J (x, t). Thus, the two terms in (27) that are not spatial derivatives must cancel each other (up to a total derivative), giving which gives the hydrodynamics for ρ nc R (26). With this relation, which nicely shows that the non-conserved part of ρ R exhibits diffusive dynamics with a drift, while the conserved part simply diffuses with no net drift.

B. Spreading of local non-conserved operators
We now briefly describe the spreading of initial operators O 0 that are orthogonal to S tot z such that Tr(O 0 S tot z ) = 0. These operators come in two categories: those with ∆S tot z = 0 (such as the raising and lowering operators r 0 and l 0 with ∆S tot z = ±1 respectively) and those with ∆S tot z = 0 (such as r 1 l 2 ). The first category of operators only have weight on basis strings with the same ∆S z = 0 for all times, and thus remain orthogonal to each conserved charge (zI) i so that a i (t) = 0 ∀ i, t. Thus, in this case, ρ R (x, t) = ρ nc R (x, t) and these operators do not have any left-or right-weight that is "left behind" near the initial location due to slow diffusive dynamics. This means that their right-weight profile does not have diffusive tails and looks essentially the same as ρ R for an operator spreading under the action of an unconstrained random circuit (13), albeit with different drift and diffusion coefficients. This can also be seen from the hydrodynamic equation (26) since the source term will be zero for this case. However, we will show in the next section that the local operator content within the spreading operator, i.e. internal to the light-cone, still shows power-law correlations due to the conservation law.
Let us now turn to the second category of initial operators with ∆S tot z = 0. These include "dipole" operators of the form O dip 01 = (zI)0−(zI)1 √ 2 which still contain the local conserved charges but with amplitudes that sum to zero. In this case, one can show that the coarse-grained Haar-averaged conserved amplitudes a c (x, t) look like the spatial derivative of the amplitudes obtained for the case starting with a monopole source (17): Indeed, this is also easily obtained from the hydrodynamic diffusion equation (25) with a dipole initial condition. Then, the total conserved weight again decreases as a power law in time, but with a faster decay as compared to the monopole case: ρ c tot (t) ∼ t −3/2 in 1d. The conserved densities again act a source of non-conserved operators leading to a power law tail in ρ R (x, t), but with scaling 1/(x − v B t) 5/2 in 1d as is seen by considering the rate of change of ρ c tot (t). For q < ∞, charge-neutral operators like r 1 l 2 create a dipole with nonzero probability at early times which then spreads as just described, again giving such a power-law tail in the asymptotic operator shape.
This brings us to a technical aside about the large q limit which is relevant for the spreading of non-conserved operators. An initial operator (rI) 1 (lI) 2 becomes a superposition of order q 4 operators, each with ∆S z = 0, under the action of the first gate. Only operators that act as the identity on the qudit spins (and as z on the spin 1/2) have a chance of making a charge dipole and thereby contributing to the conserved weight ρ c tot -but these are suppressed in probability by 1/q 4 . This illustrates one aspect of the dynamics that is suppressed by the infinite q limit, namely the likelihood for certain operators to make dipoles and thus to pick up power law tails in ρ R (x, t). As a related point, if one starts directly with the dipole operator O dip 01 as before, then the q = ∞ limit is sensitive to whether the dipole is acted upon by a single gate U 01 at the first time step, or whether it is acted upon by the two gates U −10 and U 12 . The latter case proceeds exactly as described above since each initial gate sees a conserved (zI), while in the former case the single gate immediately converts the dipole to order q 4 non-conserved operators leading to a 1/q 4 suppression in the power-law tail. This strong sensitivity to microscopic initial details in the spreading of this class of operators is a peculiarity of the large q limit.

C. Spreading of multi-local conserved operators
We now briefly address the spreading of initial Pauli strings that are a product of N conserved charges, with the sites ordered so that i 1 < i 2 < · · · < i N . The operator dynamics of such strings is also directly constrained by the conversation law since the conservation of S tot z implies the conservation of (S tot z ) N = j1j2···j N (zI) j1 . . . (zI) j N , and thus "multi-localized" strings of charges act as the "conserved densities" of this higher order conservation law. Then, as in (14), the sum of all amplitudes a c j1j2···j N of the strings of z that appear in the expansion of (S tot z ) N is conserved in time. Further, if O(t = 0) has N conserved charges, then only amplitudes a S on basis strings with exactly N conserved charges survive Haar averaging, even though the expansion (S tot z ) N involves basis strings with less than N charges (in fact the relevant conserved operator is the appropriate weighted sum of (S tot z ) N , (S tot z ) N −2 , etc., such that it only contains these basis strings of exactly N charges).
We can now ask about the time evolution of amplitudes of the form a j1j2···j N with j 1 < j 2 < · · · < j N . Note that if all the j s lie on different gates at a given time t, then we have independent diffusion of charges on each of the gates and action of the circuit simply averages all the inter-gate correlations and makes them equal. Indeed, if one starts in an infinitely large system with the charges in the initial string well-spaced such that the "diffusive cones" of the individual charges never intersect, , and these amplitudes decay with time parametrically faster than those for a single conserved charge. On the other hand, if we start with a finite density of charges such that the independently diffusing charges encounter each other on a gate, then we have to account for the "hard-core" interaction between these charges which is encoded in the invariance of (zI) i (zI) i+1 under the action of a gate. As a result of this, the diffusing charges can never pass each other and the coarse-grained problem is that of "single-file diffusion" of N particles with a hard-core contact interaction. This problem has been solved in the literature 51 , and predicts that where the sum is over all N ! permutations σ for the particle labels at time t, the initial locations of the charges are i 1 · · · i N , and we require j 1 < j 2 · · · < j N .
Finally, note that initial strings that are products of conserved operators on some sites and non-conserved operators on other sites do not have any overlap with any of the conserved charges or their higher moments. All Haaraveraged amplitudes for all basis strings will be zero in this case (although two-point correlations of the amplitudes can still survive averaging). Further, even if there is initially a diffusive "cone" of locally conserved charges in the vicinity of a conserved operator in the initial string, these quickly lose their coherence upon encountering the ballistic cones emanating from the initial locations of the non-conserved operators.

IV. INTERNAL STRUCTURE OF SPREADING OPERATORS
So far we have discussed the shape of the spreading operator in terms of the right-weight ρ R (x, t); we've used this to describe the ballistic "light cone" within which the operator has spread, and the power law tails in the distributions of right-and left-weights within the light cone. We now turn to another layer of structure in the "shape" of the spreading operator that is relevant, for example, for the resulting OTOC's. This is the weighted local distributions of the different qubit operators (I, r, l, z) within the strings S that make up the operator. For an unconstrained random circuit, the distribution of these local operators is uniform between the four operator types after the leading operator front has passed a given location, but the distribution is highly biased towards identities before the front reaches. On the other hand, different operators are inequivalent in the conserving random circuit, leading to imbalances in these local densities that persist after the front passes. Indeed, the evolution of these local densities and their correlations are constrained by two further conservation laws that encode the fact that the unitary circuit conserves the action of the initial operator under the U (1) symmetry. These are the conservation of the total "spin" s if one starts in a state with definite S tot z , and the conservation of the net "raising charge" R ≡ ∆S tot z of an operator. Further, we will see that these conservation laws are coupled to each other through suitably-defined two-point correlations of charge/spin within the operator.

A. Conservation of raising charge
The action of the conserving circuit preserves the net "raising charge" R ≡ ∆S tot z of an initial operator that starts with a definite ∆S tot z action. Concretely, for a given basis string S, the local R i = +1 if the qubit operator in S at site i is a raising operator r, while a lowering operator l has R i = −1, and z and I have R i = 0. Note that R i is neither sensitive to the qudit operator content, nor the operator content on all other sites. The unitary circuit conserves R tot = i R i for operators that start with a definite R tot , and the operator's local "raising charge" R i moves diffusively. Nevertheless, the circuit still converts a given basis string into a superposition of many strings each with different local patterns of R i (but with the same R tot ) and this choice of strings means that the diffusion of "raising charge" is subject to noise. For example, if a unitary gate U 12 acts on a two-site operator with R tot = +1 such as (rI) 1 (zI) 2 , the action of the gate produces a superposition of (4q 4 ) two-site strings of the form {(ra) 1 (Ib) 2 , (ra) 1 (zb) 2 , (Ia) 1 (rb) 2 , (za) 1 (rb) 2 }, where a and and b are arbitrary operators on the qudit and, within this ensemble of strings, the raising charge r is equally likely to be on either of the two sites acted upon by the gate (Appendix A). Of course, there is also noise from circuit-to-circuit fluctuations. If we average over the weighted ensemble of all strings within an operator that is evolved by a particular circuit and across circuits, then where denotes the joint average, and the notation S i = (ra) means that the basis string acts as r on the spin-1/2 on site i and acts arbitrarily on the qudit. The last equality defines ρ r/l (i, t) as the weight on all basis strings in the operator expansion of O(t) that locally act as r/l on the spin-1/2 on site i (one can analogously define ρ I/z (i, t)). As discussed above, the action of a gate U 12 makes R i equal on the two sites and produces a "raising charge" current ∼ ∂ x R(x, t): which is identical to the structure we had previously obtained for the dynamics of a c i (t) (15). Thus, any initially local spatial distribution of R i , for example if O 0 = (rI) 0 , spreads diffusively with diffusion constant D r = D c .

B. Conservation of spin
Likewise, there is a noisy diffusion process governing the conservation of total spin s which measures the projections onto definite S tot z states. If we rotate the local spin-1/2 operator basis to the projection "up", u = (I + z)/ √ 2 and "down" d = (I − z)/ √ 2 operators, then these are charged under s as s = ±1 respectively, while r, l have s = 0. In this basis, the operators (uI) i (uI) i+1 , (dI) i (dI) i+1 and (uI)i(dI)i+1+(dI)i(uI)i+1 √ 2 are special and left invariant by the action of all two-site gates U i,i+1 , while the others mix between themselves in a manner that locally conserves s tot = i s i on each gate. Thus, as before, the circuit conserves total s tot while generating noise, so that after averaging over circuits and strings The action of a gate makes the average spin equal on the two sites: This again means that an initial spin polarization spreads diffusively with D s = D c . Note that this equality between diffusivities and the one for raising charge above (33) are not completely trivial statements, since D c is for the amplitudes of the conserved operators that are a nonidentity only at one site, while D s /D r are for the local operator weights (amplitudes squared) of u/d's and r/l's within all strings, both conserved and non-conserved.

C. Coupling between spin and raising charge
If the initially local operator contains a nonzero net value of either of these two "charges", then that charge will only spread diffusively. Thus, if we look near the ballistically-moving operator front at long times, none of this initial charge can be there yet. Moreover, the action of the circuit is symmetric with respect to interchanging r ↔ l and u ↔ d (Appendix A), so that no new imbalances of the average charges is created due to the action of the circuit. Thus, both these average charges are zero far from the initial location at late times so that ρ r (x, t) = ρ l (x, t) and ρ u (x, t) = ρ d (x, t) for |x| √ D c t. However, for all ballistically spreading operators, there is a next layer of structure in the two-point correlations of these charge densities that evolves diffusively near the front and influences some of the OTOC's. This structure stems from the initial condition in which the local operators are the identity on all sites except those near the origin, and these identities only contribute to ρ u/d and not to ρ l/r . That is, before the front comes through, the operator is locally the identity which is an equalamplitude linear combination of all strings that contain only u's and d's. On the other hand, in the final "local equilibrium" state of the operator, all local strings are equally likely a long time after the front has passed through (up to corrections that decay as a power law in L). This equilibrium state is reached via a noisy diffusion process that turns on at each position when the front passes through that position. The two point correlations before the front passes through are R i R j = 0 and s i s j = δ ij , and these serve as the initial conditions on the diffusion process. This diffusion relaxes the initial imbalances between the relative density of local operators charged under spin and raising charge, in particular the imbalance in (ρ u (x, t) + ρ d (x, t)) − (ρ r (x, t) + ρ l (x, t)). We now see how this comes about.
To start, we look at the dynamics of inter -gate correlations. Consider two different gates acting at the same time t on sites i, i + 1 and j, j + 1. These two gates simply set all inter-gate correlations of each type equal, just as they do for the average densities: and 7 other similar equations for the other inter-gate correlations of s and R of the form R i R j+1 (t + 1) etc. For concreteness, s i s j = +1 for strings that locally look like u i u j /d i d j on sites i, j and s i s j = −1 for strings that look like u i d j /d i u j , and is zero otherwise (analogous expressions for R i R j ), and the equation above can be derived by looking at the action of the circuit gates on the different types of operators (Appendix A). The zero inter-gate correlations before the front comes through are indeed invariant under this dynamics. The coarse-grained diffusion equation for this two-point correlation for x = y is thus with R(x)R(y) (t) obeying the same equation. This is for d = 1, but the generalization of this and the results below to higher d seems straightforward. Thus far we have two separate diffusion processes governing the diffusion of "raising charge" R and "spin" s correlations. In fact, these two processes are coupled once we consider the action within one gate on "dipoles". The "dipoles" of R are net raising-neutral operators like The coupling between spin and raising charge stems from the fact that while the total density of all such dipoles within one gate is locally conserved, the different dipole species mix among each other under the action of the gate. And it is this process that allows an initial operator that contains no r's or l's to relax to the final equilibrium where all local strings are equally likely. The precise "boundary conditions" that couple the above diffusion equations within one gate i.e. at x = y are q-dependent and involve other correlations that can initially be well out of equilibrium just after the front passes.
The coupling between the charges can be simplified by looking at particular linear combinations of the R i R j  3. Densities of local operators charged under "spin" and "raising action" behind the ballistic front of a spreading operator at q = ∞ and late times (41). ρ u/d (x, t) measures the local operator weight on site x on u/d operators that are charged under spin, while ρ r/l (x, t) measures the local weight on r/l operators that have raising charge. Outside the ballistic operator front (at |x| > vBt), the spreading operator locally is purely identities which contribute equally to ρ u/d = 1/2, but do not contribute to ρ r/l = 0. The arrival of the front at a given site turns on a noisy coupled diffusion process between the spin and raising charges which relaxes the initially imbalanced densities of these charges to the final "equilibrium" value where all local operators are equally likely. Note that ρu = ρ d and ρr = ρ l in this regime since any initial imbalances in raising/spin charges spread only diffusively and are not present near the ballistic front at late times. Dashed lines plot the "coarse-grained" densities in the scaling limit (42). and s i s j correlations. If we consider R i R j + s i s j , this has the initial condition R i R j + s i s j = δ ij . This initial condition is identical to the equilibrium final state where all strings are equally likely. And if we consider the action both within one gate and between two different gates, it does not change this quantity, which is therefore time-independent. Thus, although we do have two coupled diffusion equations, it reduces to just one diffusion equation for for this initial state, which also has G ij (0) = δ ij as its initial condition. The local effective diffusivity due to the processes within one gate is, in general, q-dependent and different from that when i and j are in different gates, but at long time this affects a "region" of G ij that is a fraction ∼ 1/t of the distance (i − j) over which it has spread. Thus the leading coarse-grained long-time behavior of G ij is where the time t here is measured from the time when the front passes, say, the point (x + y)/2. Note that the effective diffusivity here is 2D c = 1, since we have independent diffusion at the locations x and y. In fact, at q = ∞, the evolution of G ij (t) within one gate is identical in structure to the intra-gate evolution and we can solve for G ij exactly in this limit (Appendix D), and the coarse-grained answer agrees with (39). Finally, note that this form for G(x, y, t) when evaluated at x = y shows how the initial imbalances in the populations of u, d vs. r, l decay diffusively away from the location of the front: For q = ∞, where v B = 1, these imbalances can be evaluated exactly and they take the form where we've taken the late-time coarse-grained limit in the last line to evaluate the sum via a saddle-point approximation (Appendix D), and restored units. At finite q, we expect the leading functional dependence to still be the same, just with a reduced v B . This equation, coupled with the fact that ρ l ≈ ρ r and ρ u ≈ ρ d away from the diffusive lump allows us to solve for the local densities of the different operators, and this is sketched in Fig. 3, showing how all local densities diffusively approach the equilibrium value of 1/4 after the front passes through. This diffusive relaxation of the imbalance leads to a power-law diffusive tail in certain OTOC's, as we will see in the next section.

V. OUT-OF-TIME-ORDER COMMUTATORS
We now turn to measuring the out-of-time-order commutator between different classes of operators (1). We find that, in several cases, this quantity is sensitive to both the "shape" and the internal structure of spreading operators. If we consider the OTOC between two generic (initially) local operators, then this quantity is close to zero until the arrival of the ballistically spreading operator front and shows a sharp increase upon the arrival of the front. However, the OTOC for systems with local conservation laws also generically develops a "diffusive tail" and approaches its late time asymptotic value only as a power-law in time. Let us now see how this arises. Define where σ {α=1,2,3} x = {r x , l x , z x } as before, we suppress the identity operators acting on the qudit spins on sites 0, x for notational simplicity, and is the equilibrium "grand-canonical" ensemble at a particular chemical potential µ. Note that while the infinitetemperature average is the only equilibrium for a random circuit model with no energy conservation, for a circuit model with an extra conservation law (like ours), we may also consider a non-zero chemical potential which weights the different spin sectors as above. We will focus on the equal-weight ensemble with µ = 0 for the majority of this section, briefly commenting on µ = 0 towards the end. We now study the OTOC for several different choices of α, β.
A. OTOC between z0(t) and rx Let us start with µ = 0, α = 3 = z and β = 1 = r. Then and it is clear that only basis strings in the operator expansion of z 0 (t) that act as (la) or (za) on site x contribute to C 0 zr (x, t). Then, using the commutation relations [z x , r x ] = 2r x and [l x , r x ] = −2z x , and the orthonormality of basis strings, we get where the notation a z S (t) denotes the amplitude of basis string S in the expansion of z 0 (t). A bit of algebra shows which is a useful way of writing this for x > 0; an analogous expression involving the left-weight ρ L can be written for x < 0. Note that the first two terms in (46) are sensitive to the overall operator shape, namely the right-weight profile of the spreading operator, while the latter two terms are sensitive to the local "internal" operator content after the front has passed site x. In random circuits with no conservation laws, all operator types are equally likely after the front has passed through and thus only the first two terms contribute substantially to the OTOC 37,38 . Notice that for this particular pair of operators, the second line measures the local raising charge R x (32) and the difference between the local densities of z's and I's at site x after the front has passed x. As discussed in the previous section, there is no average raising charge away from the diffusing "lump" near the origin, and the difference in local densities of z's vs. I's is also zero except near the origin and right at the front. So, to leading order the second line in (46) vanishes in the region between the diffusive lump and the front. Further, since we start with zero raising charge, the only contribution to the second line comes from the imbalance between ρ z (x) and ρ I (x). Indeed, all conserved operators (zI) i with i ≥ x contribute to ρ I (x), while only (zI) x contributes to ρ z (x). These conserved operators are the leading source of the imbalance between ρ I (x) and ρ z (x) near the origin. Nonconserved strings do not substantially contribute to the imbalance between ρ z (x) and ρ I (x) away from the front. Then, putting it all together, where we've separated the contributions from conserved and non-conserved operators (18) and used also the leftweight ρ nc L to produce an expression that is valid at late time for all x. Since ρ c tot (t) ∼ t −1/2 ρ c (x, t) ∼ t −1 at late times, the contribution of ρ c (x, t) is also subdominant at late times and will be neglected below, although it is visible as a weak "dimple" in the diffusive regime near the origin in Fig. 4. Thus, away from the fronts and the diffusive regime near the origin, the OTOC for this pair of operators receives its primary contribution from the total right-or left-weights alone, just as in the unconstrained circuit 37,38 . We saw above that the right-and left-weight profiles for z(t) show diffusive power-law tails, and these translate into diffusive tails for this OTOC as well.
Let us now examine this OTOC in different regimes at late times such that the leading ballistic front is well separated from the diffusive "lump" near the origin: 1. Outside the "light cone" |x| > t: Due to the locality of the circuit, a spreading operator O 0 (t) acts as the identity outside the light cone defined by |x| ≤ t. Thus, the commutator C(x, t) = 0 in this regime.

2.
Beyond the leading operator fronts so |x| v B t + D ρ t, but within the causal light cone so |x| < t: Before the main operator front gets to x, any operator weight that is not locally the identity is exponentially small in t and thus the OTOC is also exponentially small in t (for fixed x/t in this regime). This regime does not exist for q = ∞. Ref. 38 showed that the OTOC shows a nearexponential increase with time in this regime for the unconstrained circuit model, but with a positiondependent analog of the "Lyapunov exponent". One minus the out-of-time-order commutator (OTOC) between z0(t) and rx at zero chemical potential, C 0 zr , plotted against x for a system of length L = 1000 at different times t showing the different regimes discussed in the text. For |x| > t (outside the dashed vertical lines), the OTOC is strictly zero due to the locality of the circuit. In the region vBt < |x| < t, which is inside the causal light cone but before the leading front arrives, the OTOC is exponentially small (green shaded area for the latest time). The arrival of the ballistic operator front (|x| ∼ vBt) leads to a strong increase in the OTOC from a value exponentially small to an O(1) value (shaded red area for the latest time). However, diffusive tails in the operator shape or internal structure lead to diffusive power-law tails in space and time ∼ (x − vBt) −1/2 in the late-time approach of the OTOC to its final value of 1 (shaded blue area for the latest time). By contrast, for an unconstrained random circuit (not shown), the OTOC at a given site approaches one exponentially quickly after the leading front passes 37,38 . The diffusive region near the origin |x| √ Dct (shaded purple) receives a subleading 1/t contribution from the conserved charges which shows up as a "dimple" in the curves at early times which becomes weaker at late times. All curves are obtained via a simulation using q = 3 and taking into account all processes to order 1/q 2 . The dashed red curve is the q = ∞ prediction for the functional form of the tail (48).
We expect the same qualitative behavior in this regime for our model, since this operator "edge" just comes from non-conserved operators emitted at early times whose right-and left-weights then show biased diffusion, just as in the unconstrained circuit.
3. Within the leading operator fronts, |x| − v B t ∼ D ρ t: This regime describes the growth of the OTOC from an exponentially small value in t to an order one number due to the arrival of the ballistic front. Here the operator right-or left-weight is again from non-conserved operators emitted at early times and, to leading order at long time for finite q, is a Gaussian of width ∼ D ρ t, so the leading behavior of the OTOC in this regime is given by the corresponding error function that is the integral of this Gaussian profile, just as in the unconstrained circuit 37,38 . 4. In the "tails", |x| v B t − D ρ t: This regime describes the late-time approach of the OTOC to its final asymptotic value long after the main front has passed site x. In this regime, the deviation of the OTOC from its final value of one is given by the total weight (conserved and non-conserved) of operator strings that have not yet reached site x at time t. This is obtained by considering the total conserved weight at time t = (t − |x|/v B ), since any non-conserved fronts emitted after after time t do not reach site x by time t. Then, from (47) and (21), Thus, the tails in the right-and left-weights of z 0 (t) lead to tails in the OTOC and, long after the front has passed site x, the OTOC still has a power law (in x and t) deficit from its asymptotic value. By contrast, in an unconstrained random circuit, the OTOC approaches 1 exponentially in time after the front passes. This slow power-law approach characterizes the "diffusive tail" in the late time behavior of the OTOC. In the diffusive regime near the origin, |x| ∼ √ D c t, there are further corrections to this of order ∼ 1/t from ρ c (x) that smoothly connect the OTOC's x-dependence between x < 0 and x > 0. Fig. 4 shows a sketch of C 0 zr (x, t) at different times, depicting the different regimes above.

B. OTOC between r0(t) and zx
Let us now consider the OTOC between r 0 (t) and z x so that µ = 0, α = 1 = r and β = 3 = z. Then In this case, the first line is again sensitive to the overall operator shape, while the second line cares about the relative density difference within strings between local operators that are charged vs. uncharged under R, i.e. the difference in the local densities of r/l's compared to I/z's -note that this is a different combination of local densities as compared to the one in (46). Next, note that the averge OTOC between r 0 (t) and z x should equal the average OTOC between z 0 (t) and r x since the two are related by time-reversal: C 0 rz (x, t) = C 0 zr (x, −t), and the average properties of our circuit are invariant under time-reversal. In the previous subsection, we saw that the diffusive tail in the OTOC C 0 zr (48) came from the tail in the right-weight profile of z 0 (t). However, recall from Section III B that the spreading r 0 (t) shows no tails in its ρ R , since this operator starts out orthogonal to the conserved charges so that ρ R (x, t) ∼ e −(x−v B t) 2 /4Dρt / 4πD ρ t, and ρ R (x, t) is exponentially small for x v B t. Thus, the tail in C 0 rz behind the ballistic front must come from the imbalance between the density of raising charges r/l and spin charges u/d (which are superpositions of I/z) . Indeed, we saw in Section IV that the coupling of spin and raising charge combined with the initial condition at the left and right fronts leads to a diffusive decay of the imbalance ∆ Rs between the local ρ r /ρ l and ρ u /ρ d densities away from the locations of the fronts (42) (Fig. 3). Using (42), we find that the OTOC C 0 rz in the regime , which agrees with (48), as it must.
This case provides a nice example of how the OTOC can probe the internal operator content inside the light cone and show diffusive tails, even when there are no tails in the overall operator shape of the spreading operator as measured by ρ R/L .

C. OTOC between z0(t) and zx
Turning next to the "diagonal" case µ = 0, α = β = 3 = z, we find that C 0 zz has the same structure as in (49), but ρ R now refers to the right-weight profile of z 0 (t). Thus, this commutator receives contributions both from the tails in the right-weight and from the tails in the imbalance between the spin and raising charges, so for |x| v B t. Thus, a commutator between two conserved charges couples more strongly to the diffusive processes in the system, doubling the amplitude of the diffusive power-law tail in the OTOC, relative to the two other cases discussed above.

D. OTOC between r0(t) and rx
Finally, consider the commutator between r 0 (t) and r x so that α = β = 1. Note that C 0 rr has the same structure as in (46) (with ρ R referring to the right-weight of r 0 (t)). This commutator again shows a sharp increase when the ballistic front reaches site x, but now both the contributions from the right-weight and the "internal" operator content are exponentially suppressed away from the front so this OTOC does not show power law tails in the regime √ D c t x v B t and approaches its asymptotic value exponentially after the front passes through. Instead, if we probe the commutator at short distances x ∼ √ D c t, then the OTOC displays a weak effect within the "diffusive cone" from the diffusion of the initial raising charge i.e. the diffusion of the imbalance between ρ r and ρ l . The same is true for OTOC's between r(t) and l x .
This result emphasizes that the diffusive tails in the OTOC at long distances and late times are tied to the operator dynamics of the conserved charges in the system. OTOC's between operators that are both orthogonal to the conserved charges do not display diffusive tails behind their fronts. Of course, if one simply measures the OTOC between two generic local operators, then these will have some overlap with the conserved charges and show diffusive tails.
We now consider OTOC's evaluated in an equilibrium ensemble with a chemical potential µ > 0 that weights different S tot z sectors , and we will be interested in understanding the degree to which this compares with the finite temperature OTOC's that have been evaluated for Hamiltonian models. A difference between the energy conserving and S tot z conserving cases is that, on the Hamiltonian side, both the Gibbs factor e −βH and the time-evolution are governed by the same Hamiltonian. On the other hand, for the U (1) problem, the "Gibbs factor" e µN commutes with the time evolution operator U (t) but is otherwise unrelated to it. This distinction between the two becomes important when considering non-zero µ.
To start, consider again the equilibrium Gibbs ensemble for a given chemical potential: where we've switched to the (normalized) u/d basis which projects onto up and down conserved spins in the last line. In the limit µ → ∞, In this limit, ρ eq ∞ projects the conserved spins to the "all up" state so that the conserved part of the dynamics is completely "frozen out", reducing the system to just the unconstrained qudit spins. The unconstrained qudit system is chaotic and displays ballistic growth of operators with a non-zero q-dependent butterfly velocity, just as in the unconstrained random circuit model 37,38 . This illustrates how the independence of U (t) and ρ eq can lead to chaotic dynamics even in the µ → ∞ limit, while the analogous T → 0 limit for Hamiltonian systems arrests chaos. On the other hand, if one considers µ → ∞ with q = 1, then the butterfly speed does go to zero since since spin flips in this fully polarized background can move at most diffusively. Thus, the q = 1 problem as µ → ∞ is closer to the zero temperature dynamics of Hamiltonian systems. In the rest of this section, we will work at q = 1 in the µ → ∞ limit.
We start with µ = ∞ where ρ eq ∞ reduces to the projector on the |0 ≡ | ↑↑ · · · ↑ state. Then, First consider OTOC's involving z and r. Since |0 is an eigenstate of both r and z (with eigenvalues 0 and 1 respectively), any commutator involving only r's or z's must be strictly zero. Thus, the only non-trivial OTOCs are those involving [l 0 (t), z x ] and [l 0 (t), r x ] (the commutators between [z 0 (t), l x ] and [r 0 (t), l x ] will be related to these by time-reversal). Let us start with C ∞ lz . Expanding the commutator, we find: where we've made use of the fact that |0 is an eigenstate of U (t), z and r 0 l 0 in the last line, with eigenvalues 1,1, and 2 respectively. Now, every basis string that appears in the operator expansion of l 0 (t) has net raising charge R tot = −1. Since |0 is annihilated by any r i , all strings in the expansion of l 0 (t) that contain r's cannot contribute to the OTOC, and thus only strings containing exactly one l i on some site (with arbitrary I's and z's on the others) can contribute. These create single spin flips in the polarized background. Of these, only the strings which contain l x fail to commute with z x . However, due to the diffusion of "raising charge" R i discussed in Section IV, the weight of the operator on such "single l" strings decays only diffusively away from the origin, showing that this OTOC displays only diffusive (as opposed to ballistic) behavior corresponding to a zero butterfly speed. Likewise, consider the infinite µ OTOC involving [l 0 (t), r x ]: where we've made use of the fact that r x |0 = 0|l x = 0, and (lr) = (1 − z). Thus, this OTOC has an identical structure to (53) and also displays purely diffusive behavior. Perturbing ρ ∞ µ away from µ = ∞ to leading order in large but not infinite µ gives (50) Understanding the modifications to the diffusive µ = ∞ behavior of the OTOCs at large but finite µ is an interesting question that has been explored in detail in Ref. 52.

VI. PHYSICAL SYSTEMS
Let us now turn to more generic examples of thermalizing spin chains with conservation laws, and examine to what extent the universal aspects of our results continue to hold in this setting. We will first look at a periodically driven "Floquet" spin chain where energy is not conserved, but S tot z is, just like in our random circuit model. We then turn to a interacting quantum chaotic spin chain with energy conservation.

Conserving Floquet chain
We consider a one-dimensional chain of spin-1/2 qubits that is periodically driven in time with period τ between two separate Hamiltonians H z and H xy , each of which act for half the period. The time evolution operator for this system is with The parameters we choose for this model are J z = (10 + √ 5)/16, J x = J y = (5 + 5 √ 5)/16, and τ = 1, where we have found that this system is well thermalizing, even though the individual Hamiltonians in the drive are integrable 53 . The time-evolution conserves total S tot z : [U (τ ), S tot z ] = 0 and, just as in the random circuit problem, operators like r/l have a definite ∆S tot z action which is respected by the time evolution. We will measure time discretely in multiples of τ , and the time evolution of operators in the Heisenberg picture is given by O(nτ ) = U † F (nτ )O(0)U F (nτ ). Finally, there is no randomness in this problem, and the time evolution respects locality.
Of course, the fact that charges in this thermalizing model will display diffusive dynamics is well understood. We wish to probe the interaction between the diffusive charge dynamics and the ballistic operator growth, say as probed through the right-weight profiles ρ R (x, t) of r 0 (t) and z 0 (t). While ρ R for z(t) does show a tail for this model, it is difficult to tease apart the different regimes for small finite-sized systems (it only takes time t ∼ L 2 for the diffusive charges to reach the end of the chain). To aid with this, we consider a related quantity, theweight W (x, t) which measures the weight on all Pauli strings with maximum end-to-end separation between non-identity elements equal to x: (57), plotted for a system of length L = 12. The expected diffusive "lump" in the right-weight profile of the spreading conserved charge is manifested in the enhanced weight of W (x, t) at x = 0 for z0(t), even at late times (b). No such enhancement is observed for r0(t) which evolves orthogonal to the conserved charges (a). Notice also that the late-time decay of W (x, t) away from x = L is much slower for z0(t) as compared to r0(t), consistent with the presence of power law diffusive tails in the right-weight profile of the former but not the latter.
where, as before, the LHS and RHS of S define the locations of the rightmost and leftmost non-identity operators in S. This has the advantage that all conserved charges have the form z i and these are all mapped to x = 0. On the other hand, strings that are part of the ballistic front and stretched between −v B t, v B t contribute to x = 2v B t. Thus, the non-conserved contributions to this quantity still show ballistic dynamics. In a system with no conservation laws, W (x, t) for x < v B t decays exponentially with t and quickly approaches its asymptotic value which is itself exponentially small in L. On the other hand, in a system with conservation laws, this quantity decays as a power law and approaches its asymptotic value which is instead only power law small in L if one starts with an initial operator which has overlap with the conserved densities (at late times, each conserved amplitude ∼ 1/L so that the total weight on all conserved charges is ∼ L/L 2 = 1/L). This helps us tease apart the dynamics of the conserved charge at these small sizes. Fig. 5 shows the -weight profiles W (x, t) for both r 0 (t) and z 0 (t) for a 12 site chain, with both operators starting -weight profile W (x, t) (58) of spreading operators σ x 0 (t) and σ y 0 (t) in a thermalizing Hamiltonian model that conserves energy (59), plotted for a system of length L = 14. σ x 0 has non-trivial overlap with the local conserved energy density operators, which is visible in the enhanced weight of W (x, t) for x = 0, 1 even at late times (a). By contrast, σ y 0 starts out orthogonal to the conserved charges and its latetime -weight profile does not show a "lump" at x = 0, 1. While the late time decay of W (x, t) away from x = L is again faster for σ y 0 (t) as compared to σ x 0 (t), the difference is not as pronounced as the difference between z0(t) and r0(t) in the U (1) Floquet case. This is because σ y 0 (t) can develop overlap with the slow conserved charges, unlike r0(t) which is constrained to remain strictly orthogonal to the charges during its entire evoluton.
at the end of the chain to maximize the available distance for spreading. We clearly see the tail and "lump" with enhanced -weight at x = 0 for z 0 (t). On the other hand, the late time profile for r 0 (t) shows a simple decay of the -weight, with the most weight on strings stretched across the entire system i.e. those with x = L. This data shows that the central aspects of our results on operator shape seem to be born out for this more "physical" spin chain, although for such a short chain this can only be qualitatively tested.

B. Hamiltonian spin chain
We now turn to a thermalizing Hamiltonian spinchain with energy conservation, but no other continuous symmetries 43 : where J = 1, h x = ( √ 5 + 5)/8, h z = ( √ 5 + 1)/4. The local energy density operators are one and two-site Pauli strings of the form σ x i , σ z i and σ z i σ z i+1 . The conservation of total energy implies that Tr (O 0 (t)H) = const, and thus if we start with an operator with a non-zero overlap with a local energy density, then the sum of the amplitudes of the spreading operator on all the local energy density strings is conserved. A difference between this Hamiltonian model and the U (1) conserving model is that energy is not quantized and there are no special operators (like r) with definite algebras under H. Even if we start with a local operator that is orthogonal to all the conserved energy densities (like σ y i ), its operator expansion will, in time, develop some overlaps with the conserved charges. Fig. 6 shows the -weight profile for the spreading of σ x 0 (t) and σ y 0 (t) in a system of length L = 14 and, once again, we see the "lump" at x = 0, 1 only in the former. Notice, however, that the difference between σ x 0 and σ y 0 in the decay of W (x, t) away from x = L at late times is not as pronounced as it is in the Floquet problem (Fig. 5). We attribute this to the fact that the time-evolving σ y 0 (t) develops non-trivial overlap with the diffusing conserved charges (even though the sum of the amplitudes on all conserved charges is constrained to be zero), slowing down its dynamics relative to r 0 (t) in the Floquet model which always evolves strictly orthogonal to the conserved charges. This is one example of a difference between U (1) conservation with quantized charges and energy conservation where the charges are continuous.
As discussed in the previous section, these power law tails in operator shape translate into power law tails in the long distance and late time OTOC, and these have been observed numerically in Hamiltonian systems with conservation laws 41 .

VII. CONCLUSIONS
In conclusion, we presented an extensive study of the "scrambling" dynamics of local operators in chaotic quantum systems with a conserved, diffusing charge (or energy) density. A generic local operator in this setting has some weight on the conserved charges, and the late time spreading dynamics of such an operator is described by multiple coupled hydrodynamic processes. The first is the "physical" hydrodynamics associated with the diffusive dynamics of the conserved charges. We show that the total operator weight on these conserved charges decreases as power law in time which, by unitarity, necessitates a steady "emission" process that transfers operator weight from local conserved densities to nonconserved operators. This emission happens at a slow hydrodynamic rate set by the local diffusive currents of the conserved density but, once emitted, the fronts of the nonconserved operators spread ballistically and rapidly become nonlocal. The propagation of the nonconserved fronts is described via an "emergent hydrodynamics" that is biased diffusion for the one-dimensional case 37,38 , and this coupled diffusion-emission-propagation process reveals a composite picture for the operator profile, showing how the ballistic and diffusive processes in the system connect at different time and length scales. In particular, the presence of slow diffusive modes leads to the development of a power-law tail in the operator profile that reflects the weight of nonconserved operator strings emitted at later times that "lag" behind the leading ballistic front.
Our picture illustrates how reversible unitary dynamics in closed quantum systems can display dissipative diffusive hydrodynamic modes. The dissipation arises from the conversion of operator weight from locally observable conserved parts to nonlocal nonconserved parts at a slow hydrodynamic rate, a process which increases the "observable" entropy of the system. By contrast, in systems without conservation laws, any local operator is rapidly converted to non-local, so the dissipation does not appear in the late-time hydrodynamics of operator or entanglement spreading.
In addition to the diffusive tails in the distribution of operator weight, we found an additional layer of structure describing the local operator content behind the ballistic front, within the spreading operator. Outside the ballistic operator front, the spreading operator locally consists only of local identities. The arrival of the front at a given site turns on a noisy coupled diffusion process between different species of local operators, which relaxes the initially imbalanced local operator content to the final equilibrium value where all local operators are equally likely.
Once again, this relaxation happens at a power-law slow rate and contributes to the operator hydrodynamics, in contrast to the unconstrained random circuit case where the action of a single gate at the front erases the initial bias towards the identities. These power law tails in the distributions of operator weight and local operator content also lead to diffusive tails in the late-time approach of certain out-of-time-order commutators (OTOC's) to their asymptotic values.
In all, our work reveals several rich layers of physics in the scrambling dynamics of systems with conservation laws. We expect this approach (which builds on work by  of probing the dynamics of such systems via an analytically tractable constrained random unitary circuits will have broader applicability in understanding the many open questions about the fundamentals of thermalization and quantum statistical mechanics in chaotic quantum systems with conservation laws. It would also be interesting to ask which of our results can also be calculated (or extended) in a holographic setting, where they might connect with the dynamics of black holes with charges.
Related Work: Shortly before completing this manuscript, we became aware of related work by Rakovszky et. al. which should appear in the same arXiv posting 52 . While these authors take a somewhat different approach, our results appear to agree where they overlap.
transition is only one of of order q 4 possibilities. Likewise, we can see how this process introduces correlations in the amplitudes a c i and a c i+1 and |a c 1 (t)a c 2 (t)| − |a c 1 (t)| |a c 2 (t)| ∼ 1 q 4 1 t 2 . (C5)