Diffusive hydrodynamics of out-of-time-ordered correlators with charge conservation

The scrambling of quantum information in closed many-body systems, as measured by out-of-time-ordered correlation functions (OTOCs), has lately received considerable attention. Recently, a hydrodynamical description of OTOCs has emerged from considering random local circuits, aspects of which are conjectured to be universal to ergodic many-body systems, even without randomness. Here we extend this approach to systems with locally conserved quantities (e.g., energy). We do this by considering local random unitary circuits with a conserved U$(1)$ charge and argue, with numerical and analytical evidence, that the presence of a conservation law slows relaxation in both time ordered {\textit{and}} out-of-time-ordered correlation functions, both can have a diffusively relaxing component or"hydrodynamic tail"at late times. We verify the presence of such tails also in a deterministic, peridocially driven system. We show that for OTOCs, the combination of diffusive and ballistic components leads to a wave front with a specific, asymmetric shape, decaying as a power law behind the front. These results also explain existing numerical investigations in non-noisy ergodic systems with energy conservation. Moreover, we consider OTOCs in Gibbs states, parametrized by a chemical potential $\mu$, and apply perturbative arguments to show that for $\mu\gg 1$ the ballistic front of information-spreading can only develop at times exponentially large in $\mu$ -- with the information traveling diffusively at earlier times. We also develop a new formalism for describing OTOCs and operator spreading, which allows us to interpret the saturation of OTOCs as a form of thermalization on the Hilbert space of operators.


I. INTRODUCTION
The question of how quantum information spreads in a closed quantum system as it approaches equilibrium via unitary time evolution has appeared in various guises in the literature of the past decade 1-3 . While many studies focus on the buildup of entanglement between spatially separated regions, in recent years a great deal of attention has focused on different measures of the "scrambling" of quantum information, coming from the fields of high energy physics, condensed matter physics and quantum information theory. The problem of scrambling is related to the spreading of operators in the Heisenberg picture, and to the definition of "many-body quantum chaos" as put forward by Refs. 4 and 5. These effects are captured by so-called out-of-time-ordered correlation functions, or OTOCs 6 . The OTOC exhibits an initial exponential increase in time, in analogy with the exponential divergence of trajectories which defines classical chaos, in certain weakly coupled field theories [7][8][9][10][11] and in the Sachdev-Ye-Kitaev model 5,[12][13][14] . There is, however, no clear indication of such an exponent appearing in generic local lattice systems [15][16][17][18][19] . In Ref. 20 it was shown that the growth rate of the OTOC has an upper bound which is linearly increasing with temperature and which is saturated by models which are dual to black holes. Moreover, in cases where the dynamics is local, it was found that OTOCs show a ballistic spreading with the so-called "butterfly velocity".
While there is a profusion of valuable numerical work on these questions, and various, often uncontrolled forays in quantum field theory, exact results are few and far between. Recent work by some of the authors 18 (and others in Refs. 17 and 21) set to examine these questions in the context of local random unitary circuits, where a number of exact results can be derived for the average behavior of OTOCs and other relevant quantities. Most prominently, the OTOCs in these circuits were found to obey a "hydrodynamic" equation of motion, given in terms of a biased diffusion equation, which was conjectured to apply also for generic, non-integrable systems in 1D.
While the random circuits considered in these works have the advantage of giving some handle on systems with strong interactions and locality, they do not possess any conserved quantities -in particular, no conserved energy -and therefore a number of interesting questions, for example regarding the temperature-dependence of OTOCs, cannot be formulated in them. Here we address part of this shortcoming by considering circuits which conserve some local density (a U(1) charge) but which are otherwise Haar random and local, in the hope that such circuits give insight into OTOCs in more traditional energy conserving systems. Due to the presence of the conserved charge, a rich additional structure emerges, which we explore in detail below.
We show explicitly that the charge density operator evolves diffusively under the effect of the random circuit, analogously to the behavior of Hamiltonian systems in the regime of incoherent transport. As we detail below, this charge diffusion is associated with the appearance of both a power law (as opposed to exponentially quick) relaxation (Figs. 3 and 4), and a ∼ 1/ √ t − x spatial profile (Fig. 8) for the OTOCs of those operators that have a non-vanishing overlap with the local conserved charge. We will use the words "power law tails" to refer to these phenomena. Moreover, the presence of a conserved quantity allows us to formulate a finite chemical potential form of the OTOC as an analogue of finite temperature OTOCs considered in Hamiltonian systems, with the chemical potential µ playing a similar role as the inverse temperature, limiting the size of the effective Hilbert space.
Using the time-evolving block decimation (TEBD) algorithm 22 to time evolve matrix product operators 23 , we find numerical evidence that at zero chemical potential µ ("infinite temperature") OTOCs behave similarly to the behavior observed previously in random circuits 17,18 , with the addition of the aforementioned power law tails that arise as a consequence of finite overlap with the conserved charge. Further structure emerges at µ 1 ("low temperature") regarding the short-time behavior of OTOCs. The most striking feature of this is a lack of a ballistically travelling front up to times t ∼ e 2µ , which we show by developing a perturbative expansion for the OTOCs around the µ = ∞ limit. Moreover, we find that in this limit certain OTOCs can exhibit a double plateau structure, saturating to a prethermal plateau on a O(1) time scale and only reaching their thermal values at a time scale that diverges in the µ → ∞ limit.
Furthermore, we develop a novel formalism, based on superoperators, for describing OTOCs and operator spreading in a single framework. Using this formalism we show that at long times (well in excess of any diffusive time scales) OTOCs saturate to a value determined by a Gibbs ensemble on operator space, a phenomenon we term "thermalization on operator space". There are two local conserved densities appearing in this Gibbs ensemble, which we call the left/right charge density and denote byl Q ,r Q . We use the new formalism to motivate a conjectured solution for the behavior of OTOCs in the charge-conserving circuit at sufficiently large distance scales; we find that the interplay between the conserved left/right charge densities, and the ballistically propagating operator front leads naturally to the appearance of the diffusive tails in the shape of the OTOC.
The rest of the paper is organized as follows. In Sec. II A we introduce charge-conserving random circuits and then prove charge diffusion in II B. In Sec. III we turn to the discussion of out-of-time-ordered correlators. After stating our results regarding their long-time values in III A, we discuss their behavior first at infinite temperature (Sec. III B), including a discussion of hydrodynamic tails, and than in the µ 1 limit (Sec. III C). Sec. IV contains the formulation of the problem in terms of superoperators and the hydrodynamic description arising from it. After presenting earlier results in this new language in Sec. IV A we present an approximate, coarsegrained solution for certain OTOCs and comment on the resulting front shape. We conclude in Sec. V.

II. RELAXATION OF TIME ORDERED CORRELATIONS
Before attacking the problem of OTOCs in systems with a U(1) conserved charge, we consider time ordered expectation values. We first define the charge-conserving random circuits studied in this paper and then go on to show rigorously that charge diffuses in this system when averaged over many realizations of the circuit.
A. Random Local Unitary dynamics with a conserved charge Consider a spin system with a q-dimensional on-site Hilbert space H = C q . We define a global conserved chargeQ as the sum of the local charge densityQ r , given byQ = rQ r ;Q r = diag (0, 1, . . . , q − 1) . (1) However, we stress that many of the following results do not depend strongly on the definition ofQ r e.g., the charge diffusion result below goes through for any translation-invariant definition.
The random circuits we discuss are defined as follows. Consider a discrete time evolution, consisting of layers of two-site unitary gates acting on pairs of neighboring sites. Odd numbered layers act on all the odd bonds of the chain while even numbered layers act on even bonds. Each two-site gate is chosen independently from the Haar distribution over q 2 × q 2 unitary matrices which commute withQ. In practice, this means that the two site unitary U r,r+1 , acting on sites r, r + 1, is block diagonal with respect toQ r +Q r+1 , and each of the blocks is Haar random. With the definition ofQ given above, the block structure of such a two-site unitary is U r,r+1 = 2(q−1) Q=0 U Q where U Q is a Haar random unitary acting on the d Q ≡ q − |Q + 1 − q| dimensional space of states on sites r, r + 1 that have total charge Q. For example for q = 2 it has the form where the first and last blocks are 1 × 1 and the middle block is 2 × 2. In the following we will denote this charge conserving Haar measure as d H, The time evolution after an even number of 2t layers is given by where n τ = 1+(−1) τ 2 and each of the unitaries U r,r+1 (τ ), labeled by the pair of sites they act on as well as the layer/time label τ , is an independent random matrix chosen from the charge conserving (i.e., block diagonal) random ensemble defined above. The product FIG. 1. Structure of the local unitary circuits studied in this paper. The on-site Hilbert space dimension is q. Each twosite gate is an independently chosen q 2 × q 2 unitary matrix commuting with the U(1) chargeQ, defined in Eq. (1).
is defined to be time ordered. The geometry of such a circuit is graphically illustrated in Fig. 1.

B. Charge diffuses on average
In this section we show that for time-ordered correlators the random circuit reproduces the diffusive behavior observed 15,24,25 in generic interacting many-body systems with a few (in our case a single) global conserved quantity in the regime of incoherent transport. We show the diffusive spreading by directly considering the time evolution of the local charge operatorQ r in the Heisenberg picture. More generally, letÔ be an operator acting on the two-site Hilbert space of sites r and r + 1. After applying a single two-site random charge-conserving gate on these two sites the operator becomes, on averagê where P Q projects onto the sector of the two site Hilbert space withQ = Q. This result will be sufficient to derive the diffusion ofQ.
In a system with many sites we imagine acting with a series of two-site gates using the regular gate geometry shown in Fig. 1, but where each of the two site gates locally conservesQ = rQ r . We will be interested in the evolution of the local charge densityQ r . Using Eq. (2) and the fact that the measure d H,Q U is invariant under multiplication by a swap gate, we can see that Q r (∆τ ) =Q r+1 (∆τ ). This allows us to writê Thus we can think of the local charge operator as performing a random walk, such that at each application of a two-site gate it ends up on either of the two sites with equal probabilities. It is readily verified that after an even number 2t layers of the circuit it becomeŝ . At large times, the right hand side behaves like an unbiased diffusion kernel (cf random walk generating function). Note that summing the equation over all r givesQ(t) =Q(0), which is the global conservation law.
An approximate continuum formulation of the above discrete operator equations is where D is a constant independent of q 26 . Hence, on average, the local charge density obeys diffusive dynamics. In this sense our random circuit model can be thought of as a toy model for a many-body system in regime of incoherent, diffusive transport. Such behavior is expected also in clean systems at times longer than the coherence time of quasi-particle excitations (which can be very short, for example at high temperatures 15 ), or in systems that do not possess well defined quasi-particles at all 27 . Note that in our case the diffusion of charge appears directly at the level of operators, without having to refer to any particular state, indicating incoherent charge transport over all time scales. This is also shown by the single-particle Green's function, σ − 0 (t)σ + r , where σ + r is the operator creating a single charge on site r. The Green's function vanishes on average after only a single time step, independently of the state chosen, which is another way of saying that there is no coherent charge transport.
We have shown that the local charge density relaxes diffusively. As a result, at the longest times (t > L 2 /D) the charge density becomes uniform in the system. Offdiagonal operators, on the other hand, equilibrate immediately to zero on average. Both of these statements are consistent with thermalization to a Gibbs ensemble of the formρ where µ is determined by the charge density of the initial state. It is similarly possible to argue that more complicated many-body operators eventually also equilibrate to the same ensemble on average (we leave the proof of this to future work). Using this ensemble we can also make contact with more conventional definitions of the diffusion constant 15 , given by the autocorrelation function Q r (t)Q r (0) µ − Q r 2 µ in the above Gibbs state. This correlator captures the relaxation of charge to the equilibrium value. Applying the solution Eq. (4) we find that it behaves at long times as The last equation defines an effective diffusion constant D(µ) which singles out an effective time scale for charge relaxation, t D ∝ 1/D(µ) with D(µ) = 4 cosh 4 µ/2.

III. OUT-OF-TIME-ORDERED CORRELATORS
We now turn to the description of out-of-time-ordered correlators (OTOCs, for short), defined in Eq. (6), in the charge conserving random circuit. Such OTOCs have emerged in recent years as a measure of the spreading of quantum information in many-body systems. [5][6][7][8][9][10][11]15,28 For translation invariant systems they have been studied in weakly coupled 29 local quantum field theories [8][9][10]30,31 , in models for black hole scrambling 5,20,32 and more recently in local random circuits 17,18 . In all these studies it was found that the OTOC exhibits ballistic behavior with a linearly moving front, behind which it saturates to an O(1) value, even in cases where conventional (i.e., time-ordered) correlators behave diffusively. In this regard the OTOC is more similar to measures of quantum information, such as entanglement, rather than to usual correlation functions. In Refs. 17 and 18 it was shown that this can be understood in terms of a hydrodynamic description, taking the form of a biased diffusion equation in 1D, which gives rise to the aforementioned ballistic front, albeit with a front that itself broadens in time diffusively. This description was shown to hold exactly on average for random circuits without symmetries and it was conjectured to remain valid as an effective description in other chaotic systems at sufficiently large time and length scales.
One question of great interest is how the behavior of OTOCs changes with temperature. On general grounds it is expected that at higher temperatures many-body systems behave more chaotically as there is effectively more states to scramble over. In Ref. 20 an upper bound of 2πT was derived for the growth rate of OTOCs which is known to be saturated in certain holographic models. In more generic systems, however, not much is known about the detailed dependence of out-of-time-ordered correlators on temperature.
While temperature is not well-defined for the random circuits investigated in this paper, due to lack of energy conservation, we expect that the chemical potential µ can play a similar role, effectively limiting the size of the Hilbert space available for the dynamics (for example the µ → ∞ limit projects it down to a single stationary state, analogously to T → 0 in conventional systems). We therefore define the out-of-time-ordered correlator of operators V and W as whereρ µ = e −µQ /tr e −µQ is the Gibbs ensemble de-fined in Eq. (5). By expanding the commutators we get We will refer to the last term as the out-of-time-ordered part of the OTOC and to the first two terms as its timeordered part. The interesting physics of the OTOC are captured by the out-of-time-ordered part 33 , which we denote Therefore, in the following we will mostly focus in this quantity (a notable exception in Sec. III A where we discuss the long-time limit of the full OTOC, which is mostly dominated by its time-ordered part). It will be convenient to consider operators V, W with particular charge i.e., operators which are eigenstates under the adjoint action [Q, V ] = λ V V . For example, in the q = 2 case which we focus on, the one-site operators σ + , σ − , Z, 1 1 have charges +1, −1, 0, 0, respectively (in the following, Z r denotes the Pauli z operator on site r). As we show below, at non-zero chemical potential the behavior of the OTOC can depend strongly on the charges λ V and λ W . The relation holds on general grounds, decreasing the number of independent OTOCs we need to consider. Moreover, in the q = 2 case we discuss below, we can also make us of the relation Therefore we will focus solely on the OTOCs between operators ZZ, Zσ + and σ + σ + . Note that we can exchange Z r with the local charge densityQ r as which means that any correlator of the form (7) has the same behavior if we replace all occurances of Z r witĥ Q r , up to some unimportant contributions that are either time-independent or decay diffusively, as in Eq. (4).

A. The saturation values of OTOCs
Before examining the detailed time dependence of the OTOCs defined in Eq. (6), we derive some analytical results on their expected long-time behavior at different chemical potentials. It is natural to assume that over vast time scales, well in excess of the system size, our local random unitary dynamics becomes indistinguishable from non-local dynamics with the same conserved quan-tityQ. Thus we will estimate the long time value of Eq. (6) by taking where U is a unitary that conservesQ, but which is otherwise completely Haar random, without any notion of locality. Averaging the OTOC over such unitaries is expected to yield the saturation values C V W µ (t ∞ ). In the limit of large system size, provided V, W are operators with subextensive charge (which automatically holds in the case of interest where V, W are local) this approach yields: where W ⊥ ≡ W − W and W ≡ Q P Q d Q tr (P Q W ) is the part of W that is diagonal in charge. The behavior of Eq. (9) as a function µ for different OTOCs of interest is shown in Fig. 2. This result can be derived even more straightforwardly through a direct application of our formula for the Gibbs ensemble on operator space Eq. (E2) in App. E (see discussion in Sec. IV B). Note that Eq. (9) indicates that for µ = 0, if one of the operators involved in the OTOC has non-zero overlap witĥ Q, then the out-of-time-ordered part does not saturate to zero, i.e., F V W µ =0 (t ∞ ) = 0. This fact might also be of relevance for Hamiltonian systems if the operators considered overlap with the local energy density.
The above result relies on averaging over all possible charge-conserving time evolutions without restrictions of locality, which is a valid approximation at time scales long compared to the system size. One might expect that this saturation value is in fact approached on a much shorter, possibly O(1), time scale. As we discuss below, this is indeed the case at low chemical potentials, for example at µ = 0 where the OTOCs relax to the above predicted long-time values either exponentially or as a power law, depending on the overlap between the operators involved and the conserved charge. In the limit of µ 1, however, we find that the saturation of certain OTOCs can take a time which grows exponentially with µ and in the limit µ → ∞, the long-time value of the σ + σ + OTOC in an infinite system deviates from the above prediction by an O(1) value. For a finite system this means that the OTOC first saturates to a prethermal plateau and only approaches its final value on a time scale that grows as ∼ L 2 (see Sec. III C).

B. µ = 0 and power law tails
At µ = 0 we find, using TEBD numerics, that OTOCs exhibit ballistic behavior much like that which has been analytically described for random circuits without conserved quantities 17,18 . In particular, there exists a velocity scale v B such that the OTOCs are very small at v B t − |r| > 0, increase to O(1) value near the so-called "butterfly front" r = v B t, and appear to saturate for v B t − |r| < 0, as shown in Fig 3. In line with previous work, our numerics indicate the regime over which the OTOCs obtain an O(1) value (the "width of the front") broadens diffusively (∼ √ t) in time (see the last panel of Fig 3).
There is, however, an important qualitative difference compared to the case with no conserved charge. While in that case, the OTOC saturated to its final value exponentially quickly 17,18 , here we observe that OTOCs involving certain operators (i.e., ones that have a non-vanishing overlap with the total charge operator) show a slow relaxation, consistent with a power law, similarly to the hydrodynamic tails that appear for time-ordered correlators. 24,25 . In the following we give a heuristic explanation of this phenomenon in terms of operator spreading coefficients.
Consider two local Pauli operators σ α=x,y,z 0 , σ β=x,y,z 0 on the same site. At time t, σ α 0 evolves into a superposition of operators where σ ν denotes a Pauli string of the form L r=1 σ νr r , where ν r = 0, x, y, z. The out-of-time-ordered part of the OTOC at zero chemical potential than takes form where S(ν 0 , β) = ±1 depending on whether σ ν commutes or anti-commutes with σ β 0 . In the case without symmetries it was found that the measure |c µ (t)| 2 is strongly dominated by Pauli strings ν which fill most of the spatial region [−v B t, +v B t], while the total weight contained in strings of a fixed length decays exponentially, an observation that follows simply from the fact that there are exponentially more long operators than short ones. 18 . Since the operator norm is conserved, the average weight on a single string is |c ν | 2 ∼ q −4v B t . Assuming |c ν | 2 is largely independent of ν 0 for typical strings, the sum in Eq. (11) would average to 0, as the positive and negative contributions cancel. In practice not all strings have the same weight, but we expect such fluctuation to be suppressed (law of large numbers) as O( 1/q −4v B t ). Accounting for these fluctuations, and exponentially small contributions from Pauli strings well behind the front, we were able to prove 18 that F αβ µ=0 ∼ O(e −at ). The presence of conserved charges puts important constraints on some of the operator spread coefficients and thus alters the above argument significantly. In particular, expressions of the form tr(f (Q)σ α 0 (t)) are independent of time due to charge conservation for any function f . On immediate consequence for the operator Z 0 is that the operator spread coefficients of single-site Z r operators satisfy l c Z l Z0 (t) = 1 at all times, where l ranges over all l in the forward light cone of Z 0 . Indeed, as we have shown in Sec. II B, the coefficients decay on average as t −1/2 , rather than exponentially as they would without conservation laws. Since |c ν | 2 ≥ |c ν | 2 this means that these particular weights cannot decay faster than 1/t. Summing over all sites l this implies that the total weight on single-site Z operators is lower bounded by a value that decays as t −1/2 . We observe numerically that while this , compared to the exact results for the squares of the average coefficients |c The weights approach the lower bound at times t ≈ 7. Inset: the total weight contained in single-site Z operators decays in time as t −γ with an exponent γ ≈ 0.8 initially, but approaches the lower bound ∼ 1/ √ t at the longest times attainable.
weight initially decays faster (approximately as t −0.8 ), it approaches this lower bound at times ≈ 10 (see Fig. 4), and we expect that at longer times l |c Z l Z0 | 2 (t) ∼ 1/ √ t. The scrambling of quantum information in closed many-body systems has received considerable recent attention. Two useful measures of scrambling have emerged: the spreading of initially-local operators, and the related concept of out-of-time-ordered correlation functions (OTOCs). Recently, random circuits have been used to give these quantities an effective hydrodynamical description. We extend these results by considering local random unitary circuits with a conserved U(1) charge and argue, with numerical and analytical evidence, that the presence of a conservation law slows relaxation in both time ordered and time-out-of ordered correlation functions; we show that both can have a diffusively relaxing component or "power-law tail" at late times. Moreover, we consider OTOCs in Gibbs states, parametrized by a chemical potential µ, and apply perturbative arguments to show that for µ 1, the ballistic front of information-spreading can only develop at times exponentially large in µ -with the information traveling diffusively at earlier times. We also develop a new formalism for describing OTOCs and operator spreading, which allows us to describe the saturation of OTOCs as a form of "thermalization on operator space", and leads to a conjecture for their long-time behavior, wherein the conservation law results in a particular profile for OTOCs, showing a slow decay behind the main front. Our results may be relevant to systems with energy conservation, with the role of chemical potential played by the inverse temperature, even though our models explicitly break time translation symmetry.
For the OTOC (11) the above argument means an anomalously large positive contribution coming from the single site Z operators. Since we do not expect similar enhancement for any of the negative contributions, the OTOC itself should relax to its long-time value as ∼ 1/ √ t at leading order in t. The same argument also suggests a power-law relaxation for the OTOC between Z 0 and σ + r . We expect similar behavior for the relaxation of OTOCs in Hamiltonian systems for operators that have a nonvanishing overlap with the local energy density, an effect already observed numerically in Ref. 34. For OTOCs between operators that have no overlap with the conserved charge we expect that the exponentially fast saturation, seen in random circuits without conservation laws 18 , remains valid, which is in rough agreement with the TEBD data shown in Fig. 3 for the σ + σ + OTOC. We return to the issue of power law tails in the OTOC in Sec. IV, where we outline a complementary way of understanding such slow relaxation and discuss the shape of the OTOC near the butterfly front (see Fig. 8 in particular).

C. µ 1 and OTOC diffusion
We next turn our attention to the behavior of OTOCs at low fillings, or large chemical potentials, and argue perturbatively that there is an additional structure arising in this limit, wherein the ballistic OTOC front does not appear at times that are short compared to the Boltzmann factor e −µ . As discussed above, in the "infinite temperature" ensemble, i.e., at µ = 0, OTOCs are closely related to the problem of operator spreading and essentially sample over all coefficients appearing in Eq. (10) with equal measures. This explains the ballistic spreading of OTOCs, which in this language is a simple consequence of the fact that there are exponentially more long Pauli strings than short ones. However, when µ is increased the OTOC will be more and more dominated by states with a few charges. Here, we set out to explain how this affects their space-time structure and saturation behavior by considering the limit µ 1. One advantage of the approach in this section is that it allows us to access large systems, and longer time scales than TEBD. While the method is only well controlled in the µ 1 limit, there is excellent qualitative agreement between the results of this section and those from TEBD even when µ ≈ 3. Thus, we believe the results in this section could be very useful in developing a qualitative description of the early time behavior of OTOCs in low temperature strongly coupled systems (in particular systems not permitting a quasiparticle description).
In this section we develop a perturbative approach to compute OTOCs as a power series in the small parameter e −µ . We find that the terms in this perturbative expansion describe a diffusively, rather than ballistically, spreading OTOC. This diffusive behavior is exhibited by the three lowest orders of the expansion, and we conjecture that it survives up to a time scale t ∼ O(e 2µ ), at which point the perturbation theory breaks down.
As we detail in Appendix B, the average value of the OTOC, introduced in Eq. (7), can be expanded as a power series in e −µ where the O(e −N µ ) term corresponds to a process where 2N + λ V + λ W particles diffuse in the system, interacting with each other through specific interactions that can be derived by averaging over unitary gates (see Fig. 11), and the operators V and W impose certain boundary conditions on the possible initial and final configurations of the particles. Using the precise form of the interactions, derived in App. B, each term can be formulated as a transfer matrix problem of the form is a sum over all states with N + λ V + λ W particles that are compatible with the boundary conditions at time 0 (t). This transfer matrix problem than can be evaluated numerically for the first few orders of the expansion, giving us insight into the behavior of the OTOC for times significantly longer than those obtainable by TEBD numerics.
Here we detail the behavior of the σ + σ + OTOC which has a non-trivial behavior even at zeroth order, leaving the discussion of other OTOCs to App. D. This zeroth order contribution can be computed from a random walk problem involving a pair of particles that annihilate upon meeting each other, which has an exact solution as we show in App. C. Let us highlight two properties of the solution. First, as suggested by the formulation of the problem as a two-particle random walk, the spreading of the OTOC front is diffusive, rather than ballistic. This is in contrast to the behavior seen at small µ (and general expectations of ballistic propagation), and as we argue below is a property of the perturbative expansion that is in general valid up to a µ-dependent time scale. Second, if we consider the OTOC between two σ + operators at a fixed distance r as a function of time, we find a double plateau structure: it first saturates to the value 1 2 − 1 π on an O(r 2 ) timescale and only goes to zero, as predicted by Eq. (9), when the particles reach around the whole system, at times O(L 2 ). Both of these properties are illustrated in Fig. 5. As we show in App. C, this latter result, the non-commutativity of the L → ∞ and the t → ∞ limits can be understood from the fact that in an infinite system two random walkers have a finite probability of avoiding each other forever, while they have to meet eventually if the system is finite. Moreover, the mapping from the OTOC to the above random walk problem also holds if we consider similar random circuits in higher dimensions, in which case the probability for crossing paths is smaller, and the deviation from the thermal expectation value in the thermodynamic limit is even larger.
Computing the next term, which is of order O(e −µ ), we find that it increases as √ t up to an O(L) value, as shown in the inset of Fig. 6. Similarly, we find that the ratio F σ + σ + (2) also increases as √ t. This suggests that the perturbative expansion is valid up to a time scale of order t ∼ e 2µ , at which point all terms become of comparable size. Moreover, while the second order contribution  7), between σ + 0 and σ + r at infinite chemical potential. Left: The OTOC first saturates to a prethermal plateau at 1 2 − 1 π (shown by the dashed horizontal line) as 1/ √ t. Then at a later timescale t ∼ L 2 it decays to zero. At late times its value decreases as exp(−π 2 t/L 2 ). The red dashed line shows the next order prediction at µ = 5, which indicates that the plateau survives up to a time scale that diverges with µ. Right: The OTOC as a function of initial distance r for times t = 50, 100, . . . , 500 in an infinite system. The OTOC spreads diffusively and saturates to the prethermal plateau behind the front. The inset shows that the data collapse into a single curve when the position is rescaled as r → r/ √ t does lead to a speed-up of the spreading of the OTOC compared to the µ = ∞ result shown in Fig. 5, it is still diffusive as also illustrated by the same inset. This suggests that the diffusive behavior persists up to the aforementioned O(e 2µ ) time scale. Therefore, the σ + σ + OTOC will saturate to the prethermal plateau seen in Fig. 5, if µ is sufficiently large, meaning that the scrambling time 20,35 associated to saturation of the OTOC can be exponentially large in µ. As shown by Fig. 6, the expansion up to O(e µ ) agrees well with numerical TEBD results even for µ = 3 at short times t ≤ 7.  . We observe an approximate collapse of the data when F t is plotted as a function of r/ √ t, indicating that the OTOC is still diffusive in nature.
We find similar behavior for the out-of-time-ordered part of the OTOC between operatorsQ 0 andQ r (note that the OTOC F Q0Qr µ is related to F Z0Zr µ via Eq. (8)). While both the first and second order contributions decay in time as a power law and than saturate to an Ldependent value, their ratio, F Q0Qr , increases in time as t 1/2 until it saturates to a value which is linear in system size. The OTOC betweenQ 0 and σ + r on the other hand shows a more complicated, non-monotonic behavior. The data for these two cases is presented in App. D.
In summary, we find that a variety of intriguing phenomena can occur in OTOCs at early times when the available space for the dynamic is restricted by a finite, large chemical potential. The most robust of these seems to be the initial diffusive spreading of the OTOC at early times. Whether this initial behavior has some bearing on the shape of the OTOC front at later times is an interesting question for further study.

IV. SUPEROPERATOR FORMALISM FOR GENERALIZED HYDRODYNAMICS
In this section we outline a general formalism that allows us to extend the hydrodynamic description of operator spreading and OTOCs, introduced in Refs. 17 and 18, to the charge conserving case considered in this work. In order to achieve this, we develop a description of OTOCs in terms of superoperators that act on the space of operators in the original problem. We first use this formalism to re-derive the results of Ref. 18, before generalizing them to the case with U(1) symmetry. The presence of the symmetry leads to two conserved superoperators,l Q andr Q , which feature prominently in the new hydrodynamic description. We arrive at this description after coarse-graining the dynamics by considering longerrange circuits, and applying certain approximations, to be detailed below, which we expect to hold at long times and at sufficiently large distances. The resulting OTOC has both a ballistically spreading term, similar to the one observed without symmetries, and a diffusive part which leads to both the hydrodynamic tails, discussed in Sec. III B and the initial diffusive spreading we found in Sec. III C. The same formalism can be used to argue that the saturation of OTOCs, discussed in Sec. III A, can be thought of in terms of ''thermalization on operator space", an idea we sharpen in App. E.

A. Operator spreading without symmetries
We begin by showing that the central quantities of interest considered in Ref. 18, the operator weight and the OTOC, can both be thought of as expectation values of superoperators, i.e., operators acting on the Hilbert space of operators. Throughout this section we will work in the general case where the local Hilbert space dimension is q. Let {o ν } be a normalized basis of operators acting on a single site. We can think of these as states in the q 2 dimensional Hilbert space of operators L(C q ), written as o ν →|| o ν , equipped with inner prod- The space of operators on this operator Hilbert space is denoted L(L(C q )) -one should think of this q 4 dimensional space as the space of superoperators. Note that || o ν o ν || is a basis for such superoperators. In order to declutter the notation, we denote superoperators of form || A B || as A B † ; in this notation o ν o ν † is a basis for superoperators. Where it might not be immediately clear, superoperators in L(L(C q )) are distinguished from elements of L(C q ) through the use of a small circle superscriptO. Useful examples of such superoperators are where a is an arbitrary operator, andr a ,l a can be thought of as right/left multiplication by operator a so that, e.g., l arb || o =|| aob . The superoperatorp is a projection onto the identity, not to be confused with1, which is the identity operator on operator Hilbert space. The superoperators defined so far are associated with a single q dimensional site in the original Hilbert space, but they extend naturally to situations with more sites. For instance, the identity operator on a spin system with L sites decomposes as1 ≡ ⊗ L r=11r , and similarly the projectionp ≡ ⊗ L r=1pr . Using these, we can re-express the operator weight R(s), defined in Ref. 18 as the total weight of an operator left of site s, as the expectation value of the superoperator The operator weight of an operator O = ν c ν σ ν , where σ ν denote Pauli strings, is then given by Under time evolution with the unitary U , superoperators evolve asŮ † adR (r)Ů ad , whereŮ ad || o ≡|| U oU † . Using the Haar averaging formula (A3) for each two-site gate gives an exact equation of motion forR, which is indeed the correct biased diffusion found in earlier works 17,18 , with the continuum analogue The same formalism clarifies the relation between the OTOC andR. Let us concentrate on the infinite temperature case of the out of time ordered part of an OTOC, F vw µ=0 , as defined in Eq. (7), where WLOG we take w to be a traceles, unitary local operator on site r. We can recast this as an expectation value in operator space Applying Eq. (A3) on sites r, r + 1 toO w (and making use of the identities (A2)) gives The first part of this superoperator has no subsequent dynamics and can henceforth be ignored. The second part looks very similar to the form ofR, as one can see by comparingp r,r+1 =p r,r+1 ⊗ s =r,r+11s to Eq. (13). This motivates the construction of a generalization of both the operator weight and OTOC superoperators, where we assumed x < y. This superoperator measures the weight of an operator contained outside the interval [x, y]. Note thatR(y) =W (y + 1, ∞), and that after a single Haar random time stepO w (∆τ ) =W (r, r + 1), up to an unimportant constant term, and proportionality factor, so that understanding the behavior ofW (x, y) gives us access to both the OTOC and the operator weight. Using the two site dynamics in Eq. (A4), we find thatW (x, y, t) performs a left/right biased diffusion in its co-ordinates x, y respectively, as shown in Fig. 7 In a continuum limit this dynamics is described by an equation of form with the additional requirement that x < y which we keep implicit. This biased diffusion behavior leads at long times to solutions of the formW (x, which explains common origin of the diffusively broadening front behavior observed for both R(y) and the OTOCO w .

B. Hydrodynamics in the presence of a conserved charge
We now consider random circuits with the conserved chargeQ = xQ x , as defined previously in Sec. II A. The presence of the symmetry results in two new conserved superoperators on the operator Hilbert space corresponding to left and right multiplication by Q namelẙ l Q andr Q . Both superoperators are local densities, in that they can be written as sums of local superoperators, e.g.,r Q = xr Qx . To see that these superoperators are conserved, note that where the equality follows from U † QU = Q. Note that the relationl Qx (t),r Qx (t) =l Qx(t) ,r Qx(t) , together with the diffusion of the local charge Q x] , derived in Eq. (4) imply thatl Qx ,r Qx diffuse on average as well. Thus, the presence of a diffusing conserved quantity in the original many body problem leads to the presence of two new diffusing conserved quantities in operator space. In the continuum approximation we have These conserved superoperators can be used to give a new interpretation to the saturation values of OTOCs, previously discussed in Sec. III A. As we show in App. E, these saturation values, and indeed the long-time value of any superoperator, can be understood as a form of thermalization, wherein the initial state (in operator space), ρ v (0) =|| v v ||, becomes during time evolution locally indistinguishable from the "thermal" state The latter part of this expression in Eq. (18) is none other than the Gibbs ensemble with respect to conserved quantitiesl Q ,r Q . This result suggests (in a manner we detail in App. E) that when considering objects like OTOCs, or operator weights, the usual notion of thermalization should be supplemented by considering the equilibration in operator space, as defined above.
We can also make us of the conservation ofl Q ,r Q to shed new light on the diffusive saturation of certain OTOCs, discussed in Sec. III B. There, we attributed the power law relaxation of the ZZ andZσ + OTOCs to the conservation of expressions like tr (QZ r (t)). An alternative explanation of the diffusive relaxation of the Z, σ + , using the superoperator formalism, goes as follows. For concreteness consider evolving the σ + 0 opera-tor in time. Note that the corresponding state on operator space || σ + 0 has a non-uniform distribution of l Qr ,r Qr -in particularl Qr −r Qr = δ r0 initially. But we know that on averagel Qr −r Qr relaxes diffusively from Eq. Eq. (17), so we expect the initially inhomogeneous charge profile to relax as 1/ √ t in time. On the other hand, the Zσ + OTOC, up to unimportant factors and constants, is given by the closely related expression σ + 0 (t) || (l Qr −r Qr ) 2 || σ + 0 (t) ; therefore we expect this OTOC to show diffusive relaxation becausel Qr −r Qr does.
After this initial discussion we now turn to the main point of the present section, which is understanding in more detail how the OTOC front evolves in the case with charge conservation. In order to compute OTOCs at finite chemical potential, it is useful to consider a generalization of the OTOC operator introduced in Eq. (14), given byO where we introduced the notation ω µ ≡ e − µ 4 Q .Note that the expectation value of this OTOC operator on some local operator v takes the form of the out-of-time ordered part of an OTOC, i.e., v ||O µ,Z1 || v = e µ 2 (λo+λv) F ov µ (this is similar to the regularized version of the OTOC introduced in Ref. 20). Our goal is to time evolve such OTOC operators, using the charge-conserving form of the Haar averaging formula, given in Eq. (A5). However, an exact solution is considerably more difficult in the presence of symmetry, and in what follows we instead adopt aggressive approximations in order to motivate an approximate description which is qualitatively in line with our numerical results and captures the fact that the OTOC has both ballistic and diffusive components.
To arrive at this approximate description, we consider a "coarse-grained" version of the circuit, by increasing the range of the random unitary gates such that they acts on 2M consecutive sites. We label physical sites by r and super-sites (consisting of M copies of the q = 2 Hilbert space) by x. Time evolution is then described by a network of these longer range random symmetric unitaries, with a geometry similar to the original M = 1 case illustrated in Fig. 1. We find that in the limit M 1 the dynamics simplifies considerably, allowing for a closed approximate expression for the OTOC, which we detail below. We consider the case o = Z 0 here and leave the o = σ + 0 case for later work. Let us start by noting that the Z OTOC operator can be decomposed as This already suggests that is should have a diffusive component, since we already showed thatl Qx andr Qx diffuse on average. The main technical aim of this section is to understand the average behavior of the non-linear term l Q0rQ0 . Let us first evolveO µ,Z0 with one unitary gate on r = 1, 2, . . . , 2M , corresponding to super-sites x, x+1, using Eq. A5. A straightforward calculation yields a sum of two terms where we introduced the quantitiesp Q1Q2 ≡ P Q1 P Q2 , b Q ≡ 1 − (1 − Q M ) 2 and ζ x ≡ M r=1 Z r and we neglected terms that are suppressed by at least 1/M 2 . As we argue in App. F, the first term in Eq. (20) grows ballistically upon applying further layers of unitaries. For the second term, on the other hand, we find that superoperators of the forml Qx ,r Qy diffuse, unless super-sites s and y are acted upon by the same gate in the circuit, in which case extra contact terms arise. Summing up these different contributions, and applying some further approximation which we detail in App. F, we arrive at the following form of the OTOC operator at time t: where Here K x,y denotes the diffusion kernel defined by Eq. (4), and we have defined α µ ≡ 1−2M 2M cosh −2 (µ/2). Note that the time evolution of the OTOC involves summing up contributions from diffusion processes starting at different times t . This is a consequence of the aforementioned contact terms, wherein the diffusively spreading densititesr Qr ,l Qr can be converted into ballistically propagatingp µ superoperators. While Eq. (21) is only heuristically motivated, rather than exactly derived, it passes a number of consistency checks. By contruction, the OTOC saturates to the correct value predicted by Eq. (9). Moreover, it satisfies the condition that when evaluated on the identity operator 1 1, the resulting OTOC F Z1 1 is one at all times.
Applying the approximate solution (21) for the ZZ OTOC, we get where Y(t, x) stands for the double sum appearing in Eq. (21) evaluated at the operator Z x , which reads The resulting OTOC for µ = 0 is shown in Fig. 8. We find that it exhibits a tail behind the front where the OTOC is proportional to (t − x) −1/2 . These hydrodynamic tails were also studied in more detail by Khemani et. al. 36 finding a similar scaling near the front. Near the origin, the function Y relaxes as at −1 + bt −1/2 , consistent with our earlier discussion of power law tails in Sec. III B. Evaluating Eq. (21) on the local operator σ + x on the other hand gives, at µ = 0, the result The main contribution, determining the shape of the OTOC front, is given by the same function, Y(t, x), as for the ZZ OTOC, meaning that the shape of the tail behind the front is the same as the one seen in Fig. 8. Notice, however, that there is a factor of 1/2 in front of Y, such that the value of the tail at a given distance is half the value for the ZZ OTOC: this is also in agreement with the results in Ref. 36. The equation and the formalism in this section marries two notions of hydrodynamics: The first was explored in our previous work 18 , and arises due to the local conservation of quantum information and is associated withW in Sec. 15. The second arises due to the presence of local conserved densityQ, which leads to the two locally conserved, diffusing, superoperator densities. It is satisfying that our approximate solution for the hydrodynamics couples these two types of quantities: Eq. (21) includes terms which locally resembleW (namely thep terms), as well as the conserved densitiesl Qr ,r Qr , and couplings between these terms.

V. CONCLUSIONS
We investigated the dynamics of a U(1) symmetric local unitary circuit, which we propose as a toy model for chaotic many-body systems with incoherent transport properties, and proved that the conserved charge in this system obeys an exact diffusion equation on average. We provided both analytical arguments and numerical evidence that this leads to the appearance of hydrodynamic tails in out-of-time-ordered correlators of operators that overlap with the total charge.
We developed a perturbatve expansion which captures the short-time behavior of OTOCs at low filling, and found that it describes diffusive spreading, indicating that the ballistic behavior usually associated with such correlators can only develop at time scales that are sufficiently large compared to the chemical potential µ. Moreover, we found that in this regime a peculiar double plateau structure appears for the σ + σ + OTOC, wherein for large enough µ it initially saturates to a prethermal plateau before decaying to zero.
In the last part of the paper we developed a general formalism, involving superoperators, to describe the spatial spreading of operators and the evolution of OTOCs in a unified framework. We believe that this formalism will prove useful in other settings as well. A corollary of this new formalism is an interpretation of long-time saturation of OTOCs in terms of a generalized notion of thermalization for operators rather than states; in par-ticular we noted the appearance of conserved superoperatorsl Q ,r Q , whose diffusive behavior is connected to the power law relaxation of the OTOC. Using this formalism we provided motivate a conjecture for the large distance behavior of the OTOC, showing a slow, power-law decay behind the butterfly front.
An important direction for future research, touched upon at sevaral times in this paper, is finding out which, if any, of the results presented here are relevant for Hamiltonian systems. We believe that in certain regards energy conservation can play a similar role as the charge conservation studied here, and can lead, for example, to hydrodynamic tails in OTOCs, such as the ones found numerically in Ref. 34. Whether other aspects of the present work, such as slow scrambling at low temperatures, or a hydrodynamic approach similar to the one developed here, carry over to the Hamiltonian case is an interesting open question for further study.
Related Work: Shortly before the completion of this manuscript we became aware of closely related work by Khemani et. al. 36 , which should appear in the same ArXiv posting. While they take a slightly different approach, our results appear to agree where they overlap.

ACKNOWLEDGMENTS
We are grateful to Michael Knap, and especially Shivaji Sondhi for many useful discussions. TR and FP acknowledge the support of Research Unit FOR 1807 through grants no. PO 1370/2-1, and the Nanosystems Initiative Munich (NIM) by the German Excellence Initiative. CvK is supported by a Birmingham Fellowship.

Haar averages
The average of four copies of a unitary gate acting on operators a and b, which we use to derive equations of motion for OTOCs and operator weights can be written, in the superoperator language of Sec. IV, aŝ Combining this with the identities (A2) yields, for the effect of a single two-site gate in the case without symmetries, the relations A modified version of the Haar averaging formula, valid in the presence of charge conservation, readŝ We use this equation in deriving both the perturbative expansion in the following appendix, and the OTOC equations of motions in Sec. IV.
Appendix B: Perturbative expansion for µ 1 In this appendix we detail the derivation of Eq. (12), showing explicitly how to construct each term in the perturbative scheme. The starting point of the expansion is the identity where we definedṼ = e − µ 8Q V e − µ 8Q . It is useful to think of the correlator on the right hand side of this expression as a analogous to a Keldysh path integral with four time contours, corresponding to four copes of the random circuit U , as shown in Fig. 9. Each operator insertion than connects two such layers.
The way the perturbative expansion proceeds in than to write which we can think of as either having no particle or one particle on each site, with the latter option having a relative weight e − µ 8 . Expanding out these products in powers of e − µ 8 thus corresponds to boundary conditions at times 0 and t with different numbers of particles. Between these two boundaries, the particles evolve according to a set of 'Feynman rules' that result from performing the Haar averaging over each charge sector of each random gate, using the formula in Eq. (A5). Note that every matrix element of the two-site unitary has to occur twice, one with positive and once with negative time direction (i.e. complex conjugated) in order to get a non-vanishing result. As a consequence, particles are pairwise glued together: whenever there is a particle on a given site on one of the red contours (see Fig. 9) these has to be another particle present in one of the blue contours as well. This means that there are six distinct possibilities for each site at any given time: it is either unoccupied, or occupied in two of the layers, which FIG. 9. Representation of the OTO four-point function, defined by on the right hand side of Eq. (B1), as a 'path integral' involving four layers. Each layer corresponds to one of the unitary time evolution operators (red:Û ; blue:Û † ) appearing in the correlator. Each of these is given by a random circuit, and averaging over random circuits effectively 'glues' red and blue layers together in this representation.
FIG. 10. Notation of the five different 'particle types' that can occur: the first four correspond to exactly two of the four layers (shown in Fig. 9) being occupied by a charge, while the last one is a bound state, either formed by the first two or the second two particle types.
can be 11, 12, 21 or 22, or it can be occupied in all four layers (1122). In the following, it will be useful to think of the doubly-occupied cases as having a single particle which can, however, belong to four different particle types, while we will think of the case when a particle is present in all four layers as a bound state formed by two particles. We will repsecent these five cases with five different lines, one of which is a double-line signifying the presence of a bound state, as shown in Fig. 10.
Each layer of the random circuit acts a transfer matrix, evolving a configuration of particles to a linear combination of many such configurations. As long as a charge is far away from all the others -i.e., a circuit element acts on a single particle, belonging to either of the first four types shown in Fig. 10 -it simply diffuses according to Eq. (3). More complicated behavior occurs when two particles meet, or when a bound state (double line in Fig. 10) is involved in the process. A number of possible processes and the corresponding matrix elements, some of which are negative, are shown in Fig. 11. All other allowed processes can be generated from these by applying one or more of the following operations on them: a) swap the two lines on either side of the gate; b) exchange incoming to outgoing lines and vice versa; c) permute the particle types according to (1 ↔ 2) or d) permute the particle types according to (1 ↔ 2) (1 ↔ 2). All other matrix elements not in this set are zero.
Using these rules one can write down the transfer matrix of the problem in a sector with a given particle number (note that there are further restrictions since the particle number has to be conserved in all four layers individually). The insertion of the operators V and W enforces the presence of particles on certain sites at times 0 and t, which can be read off from Fig. 9. For example the choice V =Q r implies that site r has to be occupied in all four layers at time FIG. 11. Some possible one-and two-particle processes. 0, while W = σ + r forces site r to be occupied in layers 1 and 1 at time t. At order e −N µ in the expansion we have further 2N particles in the system, coming from expanding out Eq. (B2). From Fig. 9 we can deduce that at time 0 these 'vacuum-particles' have to be inserted on layers 11 or 22, while at time t they are on layers 12 or 21 (including the possibility of having them on the same site and thus forming a bound state). We can then define |ψ V W i (N ) as the sum of all initial charge configurations compatible with the boundary conditions at time 0 and |ψ V W f (N ) as the sum over all possible final states at time t. This gives us the O(e −N µ ) term in the expansion as where T is the transfer matrix whose matrix elements are given by the interaction terms shown in Fig. 11.
Appendix C: Solution of σ + σ + OTOC in the µ = ∞ limit In this appendix we show how the OTOC F σ + 0 σ + r µ=∞ (the only non-trivial OTOC in this limit) can be understood in terms of a two-particle random walk of absorbing particles, and how this description gives rise to the two important qualitative features (diffusivity and double plateau structure) shown in Fig. 5. As described in Appendix B, this OTOC, which is the zeroth order term in the perturbative expansion, is given by a process wherein the transfer matrix is sandwiched between two 2-particle states. The boundary conditions are the following: • At time 0 there is a particle on site 0 on layers 21 and a second particle on site s = 0 on layers 22 • At time t there is ap rticle on site 0 on layers 22 and a second particle on site s = r on layers 21.
As long as the two particles in the initial state do not meet they each perform a random walk process of the type described in Eq. (3). Upon meeting each other the two particles annihilate, since there is no matrix element with this specific set of incoming particles, as can be seen from Fig. 11. This means that the computation of the OTOC reduces to the following problem: Given two random walkers, one that has to start at site 0, and another which has to end up at site r at time t, what is the probability that their paths avoid each other?
The solution to this problem can be easily formulated in terms of single-particle transition probabilities, by noting that there is a one-to-one mapping between crossing paths of the two particles with the initial endpoints and arbitrary paths where the two endpoints at time t are interchanged. This mapping is simply given by reinterpreting the paths of the two particle by changing the last crossing into a reflection or vice versa (this is a simple case of the Lindström-Gessel-Viennot Lemma, see Ref. 37 and the references therein). 38 Using this trick, the solution is given by where p(r 1 → r 2 , t) is the probability of a single random walker travelling from site r 1 to r 2 in time t The problem of calculating the OTOC thus reduces to solving a single-particle diffusion problem. This is easily done in an infinite system, with the result already stated in Eq. (4). Plugging this formula into Eq. (C1) yields a solution shown in the right panel of Fig. 5, which spreads diffusively and saturates to the value 1 2 − 1 π as t −1/2 . This saturation value is non-zero because in an infinite systems there is a finite probability that the two particles avoid each other for arbitrarily long times, i.e. by travelling in opposite directions. Note that the mapping of the µ = ∞ OTOC to the random walk problem defined above is not restricted to 1D and we would end up with a similar counting of non-crossing paths for random circuits in higher dimensions. This means that the saturation value (which equals the probability of non-crossing paths) is even larger in those cases, as random walkers in higher dimensions have a larger probability of avoiding each other.
To get the full form of the OTOC, with eventual saturation to the second plateau at zero, one needs to solve the diffusion problem in a finite system with either periodic or reflecting boundaries. For a finite system of size L, and for times t L 2 /D (where D is the diffusion constant which is of O(1) in our case) the paths of the two particles have to cross eventually, and as a result the OTOC decays to zero. Here we focus on the case of reflecting boundaries, where the above way of counting crossing paths remains valid, although we checked numerically that the the results are similar for closed boundaries (the time signalling the end of the prethermal plateau is numerically larger in the case with open boundaries, reflecting that fact that particles can evade each other for longer). Instead of doing an exact solution (which is nevertheless possible), we solve the same problem in the continuum, substituting the single particle transition probabilities with the solution of the continuum diffusion equation with reflecting boundaries, where we defined p(x, t) ≡ p(0 → x, t). This equation can be solved by doing an eigendecomposition of the operator −D∂ 2 x , using eigenstates with the appropriate boundary conditions, resulting in the single-particle propagator We can than approximate the OTOC by plugging this formula into Eq. (C1). At short times, when Dt L 2 the resulting curve follows the result in an infinite system (which can be seen explicitly ba applying the Poisson summation formula the the above expression and then looking at the lowest order term in L 2 Dt ) while at times Dt L 2 it goes to zero as ∝ exp(− π 2 Dt 2L 2 ).
Appendix D: F Q 0 Qr and F Q 0 σ + r in the µ 1 limit Here we complement the results presented in Sec. III C, for µ 1 behavior of σ + σ + OTOC with analogous results on the same "low temperature" expansion of the QQ and Qσ + OTOCs. Looking at the first two terms in the expansion for F QQ µ we find that both terms decay algebraically: the O(e −mu ) term as t −1/2 and the O(e −2mu ) term as t −1 , such that the relative size of the second increases as √ t, similarly to the σ + σ + case presented in the main text. Considering the shape of the front, we find that its These results are shown in Fig. 12. shown in the left panel of Fig. 13, we observe that while the first order term has a simple algebraic decay, the second term has a somewhat more complicated structure than the ones presented in the main text. Rather than the ratio of the two terms simply increasing monotonically in time as a power law, it has an initial increase, a maximum and then a decreasing part. At an even later time scale, t ∼ L 2 , finite size effects become prominent which leads to an eventual increase to an O(L) value. Looking at the spatial structure (right panel of Fig. 13) we observe diffusive spreading of the OTOC, similarly to σ + σ + and QQ discussed in the main text.

Appendix E: Equilibration in operator space
In this appendix we use the superoperator formalism, developed in Sec. IV, to show that the expectation values of local superoperators, e.g. OTOCs, are at long times determined by a Gibbs ensemble on operator space, which reproduces the results established in Sec. III A. Consider a spin system with L sites, each with on-site Hilbert space dimension q = 2 (for concreteness). Consider evolving some pure state operator initial density matrixρ v (0) =|| v v ||. It is convenient to consider initial operators which include a Gibbs factors, e.g., take an operator of form x is a local Pauli matrix. This is useful for our purposes because the out-of-time-order part of the OTOC (the focus of our study) can be expressed as an expectation value of a local superoperator with respect to such If we apply local two site U(1) random unitaries to such a spin system for a very long time, we expect the system to scramble so completely that the time evolution is essentially a non-local random unitary with conserved U(1) charge. Hence, at long times, we expect the average density matrix to be that obtained by pluggingρ v into Eq. (A5) for a unitary that acts on the whole chain. The result is a Haar averaged density matrix of the form In what follows, we consider the expectation value of a local superoperator -for concreteness we will take a superoperator l wy r w † y where w y is a local operator. In this particular case the OTOC operator will have two separate contributions form the , ⊥ components respectively. Let us deal first with the component, tr ρ v (t ∞ ) l wy r w † y = tr ω 2γ o † w y ω 2γ o w † y = Q1Q2 e − µ 2 (Q1+Q2) tr o † P Q1 w y P Q2 o w † y .
It is readily verified by example that for o, w y local observables, this sum is for large L sharply peaked for Q 1,2 = Q + O(1) where Q = L/(1 + e µ ). (The key observation here is that local operators have O(1) charge under the adjoint action ofQ). This justifies replacingρ v (t ∞ ) with essentially any other distribution peaked in the same position. A particularly simple choice is.ρ where Z µ = tr e −µQ . The ⊥ part of the density matrix takes form We again consider the expectation values of local superoperators (e.g., l wy r w † y ). Once again, the sum is sharply peaked around Q 1,2 = Q + O(1) in the large L limit. occupies two sites, x, x + 1, where b Q ≡ 1 − (1 − Q M ) 2 . Considered as a super-operator in smaller 2 2M dimensional Hilbert space on sites x, x + 1, the individual terms give typical expectation values on local operators of size d Q e −µQ . Such summands are, for large M , dominated by Q in a small window around 2M/(1 + e µ ). The most significant term, is similarly dominated by those terms with e L,R = f L,R ≈ Q/2. Taking the approximations together gives We will probe the dynamics of Eq. F2 under unitary dynamics on x + 1, x + 2. At this point remember that the OTOC operator was originally defined on the whole Hilbert space as being (pre-) and post-multiplied by thermal factor l, r operators as in Eq. (19), and we should account for the additional factors coming from site x + 2 Another straightforward but tedious application of Eq. A5 gives two contributions to the expression for (∆τ ) = (∆τ ) 1 + (∆τ ) 2 which we detail below. The first comes from the first line of Eq. A5 making similar large M approximations as above we find (∆τ ) 1 ≈ e −3µQ α µ Z 1 µ 3 P x Q/2 P x+1 Q/2 P x+2 Q/2 P x Q/2 P x+1 Q/2 P x+2 Q/2 ≈p x,x+1,x+2 (µ), where we have Z 1 µ is the 1 super-site partition function ensemble, which is peaked at charge Q/2. The second term (∆τ ) 2 is obtained by applying the second line of Eq. A5 to ; it is sub-leading, by a factor at least O(1/d Q ), which is typically exponentially small in Q.
On net, considering the full Hilbert space, we can iterate the above procedure to argue that (t) = α µ x∈A(t)p x (µ) x / ∈A(t)l e −µQx /2r e −µQx/2 , where A(t) is a region that ballistically spreads out from initial site 1 at a velocity of 2M . We anticipate that there are O(1/M ) errors involved in this approximation associated with neglecting fluctuations in the charge arguments of the projectorsp. We leave a more thorough accounting of these errors to other works.
a. Computing (∆τ )1 Here is the straightforward but tedious application of the first line of Eq. A5 to . We act on super-sites x + 1, x + 2.
p x+1 e L f L ,p x+1 Q/2,Q/2 p x+2 e R f R ,l e −µQ x+2 /2r e −µQ x+2 /2 = e L +e R =Q1 f L +f R =Q2 χ Q/2 χ Q/2 δ e L Q/2 δ f L Q/2 e −µe R δ e R f R χ e R = δ Q1Q2 χ Q/2 Note that when we take expectation values of this quantity, we should find a value of size χ 3 where χ is the typical value of χ Q in the thermal ensemble, which is exponentially large in M for large system size. As a function of Q 1 the norm of the terms is peaked around Q 1 = Q giving where Z 1 µ is the 1 site partition function.