University of Birmingham Operator hydrodynamics, OTOCs, and entanglement growth in systems without conservation laws

Thermalization and scrambling are the subject of much recent study from the perspective of many-body quantum systems with locally bounded Hilbert spaces ( “ spin chains ” ), quantum field theory, and holography. We tackle this problem in 1D spin chains evolving under random local unitary circuits and prove a number of exact results on the behavior of out-of-time-ordered commutators (OTOCs) and entanglement growth in this setting. These results follow from the observation that the spreading of operators in random circuits is described by a “ hydrodynamical ” equation of motion, despite the fact that random unitary circuits do not have locally conserved quantities (e.g., no conserved energy). In this hydrodynamic picture, quantum information travels in a front with a “ butterfly velocity ” v B that is smaller than the light-cone velocity of the system, while the front itself broadens diffusively in time. The OTOC increases sharply after the arrival of the light cone, but we do not observe a prolonged exponential regime of the form ∼ e λ L ð t − x=v Þ for a fixed Lyapunov exponent λ L . We find that the diffusive broadening of the front has important consequences for entanglement growth, leading to an entanglement velocity that can be significantly smaller than the butterfly velocity. We conjecture that the hydrodynamical description applies to more generic Floquet ergodic systems, and we support this idea by verifying numerically that the diffusive broadening of the operator wavefront also holds in a more traditional nonrandom Floquet spin chain. We also compare our results to Clifford circuits, which have less rich hydrodynamics and consequently trivial OTOC behavior, but which can nevertheless exhibit linear entanglement and


I. INTRODUCTION
The past decade has seen a great revival of interest in the foundations of quantum statistical mechanics. This renewed interest has been driven by theoretical advances involving the long-sought demonstration that many-body localization (MBL) exists [1], as well as by the study of quantum information theory and integrable systems [2]. It has been equally driven by experimental advances in the study of cold atomic gases, which provide examples par excellence of closed macroscopic quantum systems for which the foundational questions of quantum statistical mechanics are especially acute [3]. Perhaps the broadest question has to do with identifying possible "ergodic universality classes" in quantum many-body systems and understanding their more detailed physics. Much work for such a potential classification has focused on fully MBL systems, which are believed to exhibit a breakdown of statistical mechanics [4].
In the following, we discuss OTOC behavior and entanglement growth (building on the work in Refs. [30][31][32]) for the following three spin-chain models: (I) a random circuit [see Eq. (4)], where the two-site gates are randomly chosen; (II) an ergodic Floquet system with nearest-neighbor interactions, defined in Eq. (28); and (III) a periodic Clifford circuit defined in Eq. (33). Our approach relies on quantifying operator spreading, i.e., how the support of operators changes under Heisenberg picture evolution. We derive analytical formulas for operator spreading in model (I), which we support with additional numerics, while for model (II), we rely entirely on numerical calculations. Our numerical method is based on the matrix product operator (MPO) [33] representation. Since all three types of time evolution that we consider can be represented as a network of two-site gates (see Fig. 1), the MPO can be time evolved straightforwardly by using the time-evolving block decimation (TEBD) algorithm [34]. Our results for the three models are as follows: (I) In the random circuit model (Sec. III), we find that operator spreading can be described by a remarkably simple hydrodynamical picture, which gives rise to a biased diffusion equation. Using this, we find that the typical extent of an operator grows with velocity v B , which is less than the light-cone velocity v LC , while the width of the front broadens diffusively in time (see Fig. 2). We use these results to derive exact formulas for the OTOC and entanglement growth. The OTOC travels with the same velocity v B -in the context of quantum field theory at finite temperature, the characteristic velocity scale of the OTOC has been referred to as the "butterfly velocity" (e.g., Ref. [35]), so we adopt the same terminology for v B in this work. The behavior of the OTOC near the front is also sensitive to the diffusive broadening [Eq. (20)]. At early times, before the arrival of the main front, the OTOC grows exponentially, with an exponent that increases with the initial separation of the two operators [Eq. (19)]. At long times, the OTOC saturates to 1. This result is summarized in Fig. 3. Our front propagation results also lead to an exact formula for the entanglement growth of an initial product state, from which we can extract an entanglement velocity v E [Eq. (27)]. We find that the diffusive broadening of the operator front gives rise to the inequality v E < v B . This exact result is consistent with general nonrigorous arguments [32,36], the heuristic operator-spreading model of Ref.
[ [32]], numerous results in holography [37][38][39], and the results derived for Clifford circuits in Ref. [30]. (II) In Sec. IV, we verify numerically that, for a family of ergodic Floquet circuits, there is a similar diffusively broadening front behavior as observed in the random circuit (see Fig. 7). This leads to the tentative conjecture that the diffusive front picture is valid for generic ergodic 1D Floquet spin chains, along with the resulting consequences for OTOC and entanglement dynamics. (III) Finally, in Sec. V, we compare (I) and (II) to Clifford circuits. Within such circuits, strings of Pauli operators evolve to other particular strings, rather than superpositions of such; in particular, Clifford circuits do not exhibit a diffusively broadening operator front. We connect this fine-tuned nature of local FIG. 1. Structure of the local unitary circuits studied in this paper. The on-site Hilbert space dimension is q. Each two-site gate is a q 2 × q 2 unitary matrix. For the random circuit model of Eq. (4), each gate is randomly chosen from the Haar distribution.
For the Floquet models considered in Secs. IV and V, the two-site gates are defined by the Floquet unitaries in Eqs. (28) and (33), respectively. to show the position of the front more clearly). Almost all the weight is carried by operators with endpoints at the two fronts propagating out from the initial site with speed v B ¼ ðq 2 − 1=q 2 þ 1Þ. These fronts in turn broaden diffusively in time as ffiffi ffi τ p . The two other velocity scales, the light-cone velocity v LC and the entanglement velocity v E [see Eq. (27)], are also indicated, satisfying v E < v B < v LC . The values of ρ R and ρ L after 100 layers of the circuit are shown in panel (b). Panel (a) shows the integrated operator weightsRðsÞ (LðsÞ), denoting the total weight left (right) of site s, along with the OTO commutator Cðs; τÞ. The OTOC saturates to 1 inside the front and has the value 1=2 exactly at τ ¼ s=v B .
Clifford circuits to the fact that such circuits always exhibit trivial (in a sense that we clarify below) OTOC behavior. Nevertheless, they can still have linear entanglement growth, and their local observables thermalize to infinite temperature. This demonstrates the broader point that the presence of ballistic entanglement growth and thermalization is not sufficient to predict scrambling behavior (i.e., the behavior of the OTOC). The systems above do not have local conserved quantities (in particular, they do not conserve energy), so it is surprising that hydrodynamics arises in (I) and (II). We conjecture that such hydrodynamic behavior is universal, generically appearing in 1D ergodic Floquet systems with local unitary evolution and a bounded local Hilbert space. (In this connection, our work has obvious parallels with Ref. [30], although, as we discuss in Sec. III, the hydrodynamics in our case has a different physical origin.) Moreover, our work makes a clear and precise connection between the spreading dynamics of operators, the scrambling behavior captured by the OTOC, and other metrics of ergodicity such as entanglement entropy and the late-time behavior of local correlation functions (see Appendix D).
The content of the Appendixes is as follows: Appendix A and Appendix B contain technical lemmas used in quantifying the operator weight and OTOC functions. Appendix C presents exact results on the averages of certain operator-spread coefficients. Appendix D contains exact results on time-ordered correlation functions in random circuits. Appendix E lists Haar averaging identities, while Appendix F rigorously bounds the recurrence times in a class of translation-invariant Clifford circuits.

II. QUANTIFYING OPERATOR SPREADING
Consider a one-dimensional chain of L sites, for which the Hilbert space of a single site is H site ¼ C q . There exist operators X, Z on the single-site Hilbert space obeying These generate a convenient complete basis for all operators on H site , namely, fσ μ ≡ X μ ð1Þ Z μ ð2Þ ∶μ ∈ Z ⊗2 q g. Here, μ is shorthand for the doublet μ ð1Þ , μ ð2Þ ∈ f0; 1; …; q − 1g ¼ Z q . This basis is orthonormal, such that trðσ μ † σ ν Þ=q ¼ δ μν . The operators σ μ can be regarded as generalizations of Pauli matrices, where the usual Pauli matrices correspond to the q ¼ 2 case. Generalizing this to the Hilbert space of a 1D chain, H chain ¼ ðC q Þ ⊗L , a complete orthonormal basis of operators is given by the q 2L Pauli strings, defined as σ μ ≡ ⊗ L r¼1 σ μ r r , where each string is indexed by a vector μ ∈ ðZ ⊗2 q Þ ⊗L . Our goal is to quantify how an initial Pauli string spreads over the space of all Pauli strings under local unitary time evolution. At time τ, the Pauli string σ μ becomes This defines a set of "operator-spread coefficients" c μ ν ðτÞ ≡ tr(σ ν † U † ðτÞσ μ UðτÞ)=q L . The full set of coefficients fc μ ν ðτÞg encodes all information regarding the unitary time evolution. However, as we show below, the values of particular coefficients are not so useful for accessing the most physically interesting quantities, such as entanglement entropies [30][31][32] or out-of-time-order commutators. It is more useful to consider "coarse-grained" quantities like the total weight on all operators with the right endpoint s appearing in σ μ ðτÞ, i.e., where RHSðνÞ denotes the rightmost site on which ν is nonzero [40]. Note that the "density" ρ μ R is conserved; i.e., as the operator density of the time-evolved Pauli string σ μ . In this paper, we consider systems where the time evolution can be represented as a circuit of two-site unitary gates, arranged in the geometry shown in Fig. 1. Sites of the 1D chain are indexed by s ¼ 1; …; L, while the layers of the circuit are indexed by the variable τ. It is useful to introduce coarse-grained coordinates x and t that label pairs of sites and pairs of layers, respectively, as defined in Eq. (8) and illustrated in Fig. 1. Note that, because of the geometry of the circuit, all such models have a well-defined light-cone velocity v LC ¼ Δs=Δτ ¼ Δx=Δt ¼ 1, corresponding to the fact that, with each successive time step, a local operator can spread at most one additional site in each direction.
In the following, we investigate the three models (I)-(III) discussed in the Introduction, all three of which can be represented by local circuits of the kind shown in Fig. 1. In Sec. III A, we show that for model (I) the average of ρ μ R obeys a classical biased diffusion equation, and we use this to derive exact formulas for the behavior of OTOCs and the dynamics of entanglement. In Sec. IV, we investigate model (II) and find that it shares many features with the random circuit model, such as a broadening of the propagating wavefront, which suggests that the random walk description obtained in Sec. III has applications in a wider class of ergodic systems. In Sec. V, we contrast this with model (III), where there is no diffusion and ρ μ R remains a delta function at all times, which corresponds to a nongeneric behavior of the OTOC.

III. RANDOM CIRCUIT MODEL
In the present paper, we are mostly concerned with one-dimensional local random unitary circuits, with the geometry shown in Fig. 1. Random circuits were also investigated in Ref. [30] with regards to the growth of entanglement from an initial product state, albeit with a different geometry, where the location of the unitary gates is randomly chosen, instead of the regular arrangement used here. In that work, it was argued that the evolution of entanglement obeys an equation belonging to the KPZ universality class, which determines certain universal exponents that appear in the average value and fluctuations of the entanglement entropy. Here, we shift our focus from states to operators and derive exact results for their spreading, for arbitrary on-site Hilbert space dimension. In Sec. III D, we relate our operator-spreading results to the dynamics of bipartite entanglement, as captured by the second Rényi entropy, and we find no sign of the universal KPZ fluctuations observed in Ref. [30]. It could be the case that the fluctuations mentioned only become apparent on time scales longer than those accessible in our numerics. Alternatively, they might not be present at all, and the KPZlike behavior observed by Nahum et al. could be specific to their setup [41].
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 twosite gate is chosen independently from the Haar distribution over q 2 × q 2 unitary matrices. The time evolution after an even number of 2t layers is given by where n τ ¼ ½1 þ ð−1Þ τ =2 and Wðs; τÞ is a Haar random two-site unitary acting on sites s, s þ 1. The product Q 2t;← τ¼1 is defined to be time ordered. Such a circuit is graphically illustrated in Fig. 1.
The primary goal of this work is to quantify the spread of operators under random circuits, Eq. (4), and to relate operator spread to entanglement growth. A related question is how correlation functions of local observables behave in this random circuit model in the thermodynamic limit. As we confirm in Appendix D, such correlations tend to their infinite-temperature values at long times, similarly to the case of Floquet ergodic systems [42][43][44][45][46]. This result holds for any random realization of the circuit.
Focusing on the problem of operator spreading in this random circuit model, we find that the average of the operator density ρ R , defined in Eq. (3), performs a biased random walk, independent of the internal structure of the operator considered. Solving the random walk problem allows us to derive exact formulas for the averages of outof-time-order commutators and entanglement growth, which we detail in Secs. III B and III D, respectively. In Sec. III C, we numerically quantify the fluctuations between different random realizations of the circuit.

A. Random walk dynamics of operator density
In the following, we quantify how operators spread under the time evolution generated by the random circuit defined above in Eq. (4). We focus on the average of the operator density, defined in Eq. (3), for which we derive an exact equation of motion. Upon solving this equation, we find that the operator density moves in a front whose velocity v B is an increasing function of the on-site Hilbert space dimension and with a front width increasing diffusively in time.
We start by noting that under Haar averaging, the operator-spread coefficient c μ ν ðτÞ vanishes for any time τ ≥ 1, provided that μ is nontrivial, since c μ ν and −c μ ν have equal probability. However, the average of its modulus squared, jc μ ν ðτÞj 2 , can be nonzero. (An explicit expression for this quantity is written down in Appendix C, using a mapping to a classical Ising model, but we will not require it for the subsequent discussion). Following Eq. (3), we define the average operator density as As we show below, this quantity satisfies an equation of motion, Eq. (9), which does not depend explicitly on μ.
Preempting this, we drop the explicit μ dependence ρ μ R → ρ R to declutter notation. In fact, μ will enter considerations only as an initial condition on the operator density which is the same for all initial operators sharing the same right endpoint.
To understand how ρ R evolves in time, consider the effect of applying a single two-site gate on sites s and s þ 1.
There are q 4 − 1 nontrivial operators acting on this twosite Hilbert space. Of these, q 2 − 1 contribute to ρ R ðs; τÞ (the ones that are trivial on site s þ 1), while the other q 2 ðq 2 − 1Þ contribute to ρ R ðs þ 1; τÞ. Under a two-site Haar random unitary transformation, all the possible transitions between any of these q 4 − 1 operators have, on average, the same probability [6]. The upshot is that after the application of the unitary gate, the density ρ R evolves as with probabilities p¼ðq 2 =q 2 þ1Þ and 1 − p ¼ ð1=q 2 þ 1Þ.
To apply a similar argument for two subsequent layers of the circuit, it is useful to redefine the density by grouping together the pairs of sites on which the first layer of the circuit acts. We abuse notation and denote this quantity as where we now only consider the value of the operator density at even time steps τ ¼ 2t. Applying Eq. (7) for two layers, we arrive at the equation Thus, the right endpoints of Pauli strings perform a biased random walk on the lattice, where in each step, they move to the right with probability p 2 , to the left with probability ð1 − pÞ 2 , and stay on the same site otherwise. A feature of the above equation is that the time evolution of ρ R is independent of the internal structure of the operator, and thus the solution ρ R ðx; tÞ will be the same for all initial Pauli strings, modulo a shift x → x − x 0 , where x 0 is defined by the right endpoint of the initial string. The result of the random walk process outlined above is a front that propagates to the right from its initial position Thus, the speed at which the operator weight travels is smaller than the light-cone velocity: v B < v LC ¼ 1. This resonates somewhat with the result of Ref. [35]. The width of the front increases in time as Note that in the limit q → ∞, the "particle" described by ρ R ðx; tÞ hops to the right with probability 1 in each step, and consequently, the front becomes infinitely sharp with velocity v B → v LC ¼ 1.
The total weight of left endpoints, ρ L ðx; tÞ, obeys a similar equation except that it propagates to the left with velocity −v B , while diffusing at the same rate, as shown in Fig. 2. This means that at time t, the vast majority of quantum information initially stored in σ μ with left (right) endpoint x l (x r ) is carried by operators with support ½x l − v B t; x r þ v B τ, but where the precise position of either endpoint can be uncertain within a region of width Δx ∼ ffiffiffiffiffi ffi Dt p . We can find the full distribution of ρ R ðx; tÞ using a standard generating functional method. In the rest of this section, we use coordinates relative to the initial position of the front, i.e., x − x 0 → x. The solution to Eq. (9) then reads In the scaling limit, t, x → ∞, but keeping x=t ≈ v B fixed, this becomes (using Stirling's approximation) Thus, the traveling front has the shape of a Gaussian, as one would expect from the solution of the continuum limit of the lattice diffusion equation (9).
As we see in the next section, it is also useful to compute the total weight of all Pauli strings contained entirely to the left of position x. This quantity, which we denote by RðxÞ, is given by Around the position of the front, where x ≈ v B t, we can integrate Eq. (11) to obtain where erfðxÞ is the error function. Later in Sec. III B, we will also need an approximation for R well away from the front. Using the fact that in the large t limit ρ R ðx; tÞ increases sharply with x for x=t < v B , the sum in Eq. (12) is dominated by its largest term (i.e., y ¼ x). Using the fact thatRðx ≥ v LC tÞ ¼ 1, we can similarly approximateR for x=t > v B to obtain where Θ is the Heaviside step function. This result is accurate up to multiplicative Oð1Þ constants when jx − v B tj=t ¼ Oð1Þ in the large t limit. See the discussion in Appendix A for a precise statement and derivation of Eq. (14). Using our results for the coarse-grained density ρ R ðx; tÞ, we can also write a formula for the density in terms of the site coordinate s. Note that because of Eq. (7), the ratio Using this, the density on site s after applying an even number of layers becomes where n ¼ 0, 1. We can use Eq. (15) to derive the total operator weight to the left of site s, i.e., RðsÞ ¼ P r≤s ρðrÞ, which, as we see in the next section, is closely related to the OTOC between sites 0 and s.

B. Behavior of out-of-time-order commutators
We relate our results for the time evolution of operator weights to another oft-used measure of information spreading in many-body systems, the so-called OTOC [8,9,18,27,[47][48][49]. For concreteness, consider the following OTOC between two Pauli operators separated by distance s (in this section, we work in a shifted coordinate system where one of the Pauli operators resides at site 0), where s and τ are the original time or lattice coordinates [as opposed to the coarse-grained coordinates t, x below Eq. (8)]. Here, Z r denotes the generalized Pauli operator introduced in Eq. (1), situated on site r. We show how the OTOC Cðs; τÞ behaves in the scaling limit τ → ∞, with κ ≡ s=τ held fixed. As we detail below, for s outside of the light cone (1 < κ), it is zero. As s enters the light cone (κ < 1 and close to 1), it increases exponentially. When s is near the operator front (κ ¼ v B < 1), the OTOC becomes Oð1Þ. After the front has passed (κ < v B ), the OTOC exponentially saturates to the value 1 with an exponent that is independent of s. See Fig. 3 for a summary. LetC denote the average of the OTOC over all unitary circuits, with the geometry shown in Fig. 1. Note that because of the averaging, this quantity is independent of the choice of Pauli operator; i.e., it would be the same if we replaced either or both operators in Eq. (16) with another local Pauli that is different from Z. We are concerned with the second term in Eq. (16), which equals where e θ μ;Z s is a qth root of unity arising from commuting σ μ past Z s , and c μ 0 ðτÞ are the operator-spreading coefficients of Z 0 ðτÞ. Notice that the Haar average forces μ ¼ ν, which removes all dependence on the particular initial state ψ [50]. In particular, the average OTOC value in any state, pure or mixed, will be identical to the average OTOC value at infinite temperature, i.e., trð 1 2 j½Z 0 ðτÞ; Z s j 2 Þ=2 L . At this point, we can use Eqs. (17) and (C2) to write an exact closed-form expression for the OTOC. However, instead of doing that, we write a more manageable asymptotic expression for Eq. (17) using simpler results from Sec. III A.
To perform the sum over Pauli strings in Eq. (17), we first need to prove the following statement: jc μ 0 ðτÞj 2 depends only on the position of the two endpoints of the string μ. The proof goes as follows. First, it is easy to verify that under Haar averaging, jc μ ν ðτÞj 2 ¼ jc ν μ ðτÞj 2 for any ν, and indeed for "ν ¼ 0" in particular. This implies that the average probability of the one-site operator Z 0 evolving into a specific string μ is the same as the probability of string μ evolving into Z 0 . In the random walk picture, this latter process corresponds to both left and right endpoints, which end up on site 0 at time τ during their respective random walks. As we noted previously, these random walks are independent of the internal structure of the initial string. Thus, jc μ 0 ðτÞj 2 depends only on where the two endpoints of μ are located. We confirm this argument more concretely with an explicit expression for such operatorspread coefficients in Appendix C.
The above statement has important consequences for the OTOC. If site s lives in the support of μ, then the contribution to Eq. (17) coming from the strings with the same support as μ has an equal weight for each possible value θ μ;Z s ∈ ð2π=qÞf1; …; qg, so the cosine term averages to zero. The remaining part is the total weight due to Pauli strings that are supported on intervals that do not contain site s (along with some corrections for Pauli strings that border on site s). Deferring the full justification to Appendix B, the upshot is that, provided κ > 0 in the τ → ∞ limit, the OTOC behaves as up to exponentially small corrections in τ. Hence, the OTOC physics is directly related to the operator density, and it changes appreciably at the operator front s ¼ v B t, as we show in Fig. 3.
Let us now summarize the behavior of the OTOC as a function of space and time, as parametrized by the ratio κ ≡ s=τ and taken in the limit τ → ∞. We distinguish four regimes of OTOC behavior, which we illustrate in Fig. 3.
(1) Trivial OTOC at early times (1 < κ): In this regime the events ðτ; 0Þ, ð0; sÞ are causally disconnected, so the commutator in Eq. (16) (and hence the OTOC) is exactly zero. (2) Early OTOC growth (v B < κ < 1): This regime describes the behavior after site s has entered the light cone but before it encounters the main operator front. Here, we approximate the OTOC using Eq. (14), soCðs;τÞ ≈ c 1 ρ R ðs − 1;τÞ þ ρ R ðs;τÞ, where c 1 > 1 is bounded in the s; τ → ∞ limit. Fortunately, a simple closed-form expression already exists for ρ R , namely, Eq. (15). We obtain a more convenient expression for the initial OTOC growth by expanding Eq. (15) near the light cone in the δ 2 =s → 0 limit, where δ ≡ τ − s is the distance between s and the light cone, up to multiplicative Oð1Þ constants, where . This formula demonstrates that the OTOC will increase with an exponent λ ∼ log s for 0 < δ ≪ s. Because of its dependence on s, and its limited range of validity, it is unclear whether this should be viewed as a Lyapunov exponent as in Ref. [49]. Note that the exponential increase occurs in a regime where the overall scale of the OTOC is still exponentially small in s. In the regime where the OTOC increases to an Oð1Þ value (i.e., when the operator front hits, see next point), its behavior is not exponential. Furthermore, we note that γ ∼ 1=q 2 for large q, such that the regime in which the exponential behavior can be observed becomes smaller in the large q limit.
(3) Near the front (jκ − v B j ¼ Oð1= ffiffi ffi τ p Þ): As mentioned, the above approximation breaks down when the main front, which we recall travels at speed v B and has width ∼ ffiffi ffi τ p , arrives at site s. In this intermediate regime, we estimate the OTOC by combining Eqs. (18) and (13),C This formula describes the behavior of the OTOC in the regime when it increases from a value exponentially small in s to an Oð1Þ number. (4) Late times (0 < κ < v B ): After the main front has passed, the OTOC relaxes exponentially to 1. Expanding Eq. (14) for fixed s − v LC τ and large τ, we find that the OTOC in this late-time regime is Thus, the OTOC decays to its equilibrium value with an exponent logð1 þ q 2 =2qÞ.

C. Fluctuations from circuit to circuit
The results discussed above concern quantities averaged over many different random circuits with the same geometry but different choices of two-site gates. The question remains as to whether these average quantities are also "typical"; i.e., how large are the fluctuations between different realizations of the random circuit? In this section, we investigate this problem numerically. Our numerical method relies on representing the operator Z 0 ðtÞ as a MPO, which allows us to apply the two-site unitary gates efficiently. Two layers of the random circuit can be applied by just a single step of the TEBD algorithm, which allows us to go up to bond dimension χ ¼ 20000. Both the infinite-temperature OTOC and the total operator weight contained in an arbitrary subregion can be extracted straightforwardly from the MPO representation [both calculations are similar to computing the overlap of two matrix product states, but in the computation of RðsÞ, only the legs corresponding to sites ≤ s are contracted].
To quantify the fluctuations, we look at an ensemble of 100 random circuit realizations (which is enough to reliably reproduce the exact average quantities; see Fig. 4) with onsite Hilbert space dimension q ¼ 2 and compute (a) the OTOC Cðs; τÞ defined in Eq. (16) and (b) the total operator weight Rðs; τÞ of Z 0 ðτÞ contained within the region to the left of site s. Both Rðs; τÞ and Cðs; τÞ are functions of the distance s and the number of layers τ. We find that for both quantities, the circuit-to-circuit fluctuations are largest at the traveling wavefront and become smaller deep behind it, as shown in Fig. 4. We also see that there is a well-defined front for the information propagation in each individual circuit.
We also find that the fluctuations decrease in time. Figure 4(c) shows the standard deviation of the weight RðsÞ for different times. We find that the maximum of this standard deviation over all values of s decreases in time, approximately as ∝ τ −β with an exponent 0.4 < β < 0.5.

D. Relationship to entanglement spreading
Another question closely related to operator hydrodynamics is the problem of entanglement growth. An initial product state develops spatial entanglement during time evolution. In systems without quenched disorder, the entanglement is expected to grow linearly in time, with a growth rate characterized by the "entanglement velocity" v E . In this section, we use our results for operator spreading to compute the time evolution of the second Rényi entropy S ð2Þ between two sides of a spatial entanglement cut and extract the entanglement velocity from it. (More precisely, we calculate a closely related quantity, namely, the logarithm of the average of the exponentiated second Renyi entropy S ð2Þ exp ≡ − log e −S ð2Þ , which would coincide with the second Renyi entropy in the absence of statistical fluctuations). We find that the resulting entanglement velocity is smaller than the butterfly velocity for any finite q, and it approaches the light-cone velocity logarithmically slowly, so that v LC − v E ∝ 1= log q for large q. At long times, S ð2Þ exp saturates to its maximal value with the saturation becoming increasingly sharp as q is increased.
Consider an initial ferromagnetic product state of the 1D chain, where the state on site s is an eigenstate of the local Pauli operator Z s with eigenvalue þ1. (Note that for the average behavior of the random circuit, the choice of initial product state is unimportant.) The density matrixω corresponding to this state is then a sum over all possible Z strings, i.e., Pauli strings that only contain powers of the operator Z on each site: The density matrix at time t is obtained by replacing each Pauli string σ ν in Eq. (22) with its time-evolved counterpart σ ν ðtÞ.
Let us now divide the system into two regions, A and B, the first of which corresponds to sites 1; …; L A . Generalizing the formula of Refs. [31,32], the second Rényi entropy S ð2Þ ≡ − log trðω 2 A Þ of the reduced density matrixω A ¼ tr B ðωÞ is related to the operator-spreading coefficients by where the strings ν and ν 0 are both Z strings and μ has support entirely in subsystem A. In the last equality of Eq. (23), we assumed that the off-diagonal contributions are negligible, which becomes exactly true in the random circuit model once we average over different realizations. Let us assume that L A is even. Reverting back to the coarse-grained position x [see Eq. (8)], we recognize Eq. (23) as the total operator weight in region A, Rðx ¼ L A =2; tÞ, as defined in Eq. (12), summed over all initial Z strings. As we noted previously, this quantity is, on average, the same for all initial strings with the same endpoints x 0 . The number of different Z strings with right endpoint x 0 is q 2ðx 0 −1Þ ðq 2 − 1Þ. After averaging over random circuits, and assuming an even number of layers, Eq. (23) thus becomes where we have used that ρ R ðxÞ [and consequentlyRðxÞ] only depends on the position x relative to the initial endpoint x 0 . The first term in Eq. (24) is the contribution coming from the identity operator, which is responsible for the saturation of the entanglement at long times.
Using the exact solution Eq. (10), one can perform the sum over initial positions to find The sum over binomial coefficients can be expressed in terms of a hypergeometric function.
Equation (25) describes an entanglement that initially increases linearly with time and saturates to the maximum value L A log q at long times. For τ ≪ L A , we find e −S ð2Þ ðτÞ ≈ 2q from which we can identify the entanglement velocity [51] v Recall that S ð2Þ exp ≡ − log e −S ð2Þ . Note that the entanglement velocity approaches 1 logarithmically slowly for large q, i.e., v E ∼ 1 − logð2Þ= logðqÞ. This is a separate velocity scale, distinct from and smaller than the front speed v E < v B . This difference comes from the diffusive broadening of the operator wavefront. First, it is straightforward to verify that if the wavefront is sharp, i.e.,Rðx; tÞ ¼ Θðx − v B tÞ, then Eq. (24) gives v E ¼ v B . Second, we have checked that Eq. (24) gives v E ¼ v B even if the wavefront has a width that is finite but independent of time [52]. Hence, we attribute the difference between v B and v E to the fact that the operator front broadens in time.
In the right panel of Fig. 5, we compare the exact formula Eq. (25) to the second Rényi entropy as computed numerically (using a matrix product state representation), averaging over 100 realizations of the circuit, and we find extremely good agreement. Moreover, the numerical calculation allows us to compare S ð2Þ exp and the mean Renyi entropy S ð2Þ avg ≡ S ð2Þ , respectively. We find no significant difference between the two values, showing that there are no strong circuit-to-circuit fluctuations in the entropy and both are captured well by our exact formula, at least for the time scales accessible in the numerics. We also found numerically that replacing the Rényi entropy with the von Neumann entropy leads to a slightly larger entanglement velocity.
The entanglement saturates when the contribution of the identity becomes significant (i.e., when all other operators have essentially left the subsystem). Note that the saturation softens, compared to the prediction of the simple operatorspreading model of Ref. [31], which is another consequence of the diffusive broadening of the front. This intermediate saturation regime becomes smaller with increasing q, as shown in the left panel of Fig. 5.

IV. COMPARISON WITH THE KICKED ISING MODEL
A natural question that emerges in relation to the results stated above is to what extent they are representative of other, more generic, thermalizing, quantum many-body systems. To address this question, we investigate a system with a periodically driven, nearest-neighbor Hamiltonian. Our model has the same geometry as the random circuit, shown in Fig. 1, and it similarly does not conserve energy. However, unlike random circuits, it is periodic in time, and its two-site (and one-site) gates take a specific form, rather than being randomly chosen. Despite these two significant differences, we find that several details of the operator spreading described in Sec. III, such as the diffusive broadening of the wavefront, remain approximately valid. These features have also appeared in more recent numerical studies of ergodic Hamiltonian spin chains [26,29].
For concreteness, we consider a model with on-site Hilbert space dimension q ¼ 2 that consists of switching back and forth between two Hamiltonians, such that one period of the time evolution (with period time T) is given byÛ This system-which we refer to as the "kicked Ising model"-is known to be ergodic, provided that both g and h are sufficiently large. Since at any given time the terms in the Hamiltonian all commute with each other, the time evolution can be represented as a circuit of two-site unitaries (with the one-site rotations included in the twosite gates) with the same geometry as in Fig. 1. As such, it is in fact contained among the ensemble of random circuits considered before. The question is to what extent the properties of this specific circuit coincide with the average quantities discussed above. At first, operator spreading in the Floquet system seems very different from the case of the random circuit. There is no inherent randomness, and the evolution is completely determined by the internal structure of the initial operator σ μ ; however, for the random circuit, the average behavior is independent of the internal structure. However, it is well known that ergodic systems can behave as baths for themselves [53,54] and, in effect, generate their own noise. If the same reasoning can be applied to the question of operator spreading in ergodic systems, then it is plausible that the noise-averaged equation (9) holds in deterministic ergodic systems, and as a result, the diffusively broadening ballistic front picture continues to apply.
Another more hydronamical motivation for the result is to note that the random unitary and deterministic Floquet systems both involve evolution under a local unitary circuit. Locality puts strong constraints on the evolution of ρ R : Not only does it obey the global conservation law R dxρ R ðx; tÞ ¼ const (we revert to a continuum notation for ease of presentation), it should also obey a local conservation law for some local current density Jðx; tÞ. This conservation law puts severe restrictions on the equation of motion of ρ R . For example, one can imagine that in a coarse-grained picture, on long enough time scales, the constitutive equation J ∼ vρ R þ D∂ x ρ R þ Á Á Á becomes valid; note that the discretized version of this constitutive relation is exactly the random walk equation we derived for the random circuit averaged ρ R [see Eq. (9)]. Therefore, it seems plausible that in a sufficiently coarse-grained picture, the dynamics might be well approximated by a biased diffusion similar to the one described in Sec. III for the Floquet circuit, with hopping probabilities depending on the microscopic couplings. Here, we present numerical evidence in support of this conjecture. Our results can be summarized in three points: (i) The butterfly velocity v B depends strongly on the coupling g and can be tuned to be much smaller than the light-cone velocity v LC . (ii) When tuning the couplings to decrease v B from its maximal value v B ≈ v LC , the front also becomes wider, as expected for a random walk when increasing the probability of hopping to the left at the expense of the probability of hopping right. (iii) The operator wavefront gets wider during time evolution, with the width increasing in time as ∼t α , with an exponent 0.5 ≲ α ≲ 0.6. We find numerically a linearly propagating wavefront for the time-evolved operator Z 0 ðtÞ, which shows up in both the OTOC and the weight RðsÞ, with the OTOC Cðs; tÞ saturating to 1 behind the front. While for the random circuit the speed of the front was set by the on-site Hilbert space dimension q, for the kicked Ising model we find that this speed can be tuned continuously by changing the value of the transverse field g [55], as shown in Fig. 6. Note that changing g does not affect the light-cone velocity, which is Δs=ΔT ¼ 1 because of the geometry of the circuit that represents the Floquet time evolution. For g ¼ 0, an initial operator Z 0 remains localized on the same site. As we make g larger, the butterfly velocity gradually increases, and it reaches v B ≈ v LC for g ≈ 0.9 with period time T ¼ 1.6 [56]. Looking at Fig. 6, we notice that decreasing v B from its maximum corresponds to an increased front width at any given time. This is consistent with a coarse-grained random walk description, wherein increasing the probability of hopping to the right results in both a larger butterfly velocity and a suppression of the diffusion constant. We note that the dependence of v B on microscopic parameters is expected on general grounds (see, e.g., Refs. [35,58]), although, to our knowledge, this is the first study of this dependence in infinite-temperature ergodic systems.
The most important evidence in support of a hydrodynamic description comes from examining the front width as a function of time. Similarly to the random circuit model, we find that the wavefront of the operator spreading broadens as we go to longer times. To quantify the width, we look at the standard deviation of ρ R ðsÞ ¼ RðsÞ − Rðs − 1Þ, i.e., As shown in Fig. 7, at long times, the width grows algebraically in time as σðtÞ ∝ t α , with an exponent 0.5 ≲ α ≲ 0.6. This is roughly consistent with the random walk description of operator spreading put forward in Sec. III. This diffusive broadening is expected to result in the strict inequality v E < v B for the entanglement velocity, according to the arguments put forward in Sec. III D. We confirmed numerically that this indeed holds in this model for various values of g.
Finally, one might wonder whether the above story continues to hold when we consider a system with energy conservation, i.e., a time-independent Hamiltonian. Recent unpublished work [26,29] shows evidence of a diffusively growing front in energy-conserving ergodic spin chains. The longitudinal field is fixed at h ¼ 0.809, while the period time is T ¼ 1.6. The butterfly velocity shows a strong dependence on the coupling g, with the front width increasing as one moves away from the limit of maximal velocity.

V. FRACTAL CLIFFORD CIRCUITS
In this section, we compare the results obtained for the random circuit model of Sec. III to another set of circuits that do not, in general, exhibit energy conservation, namely, Clifford circuits. We show that despite the fact that Clifford circuits can be "ergodic" in certain senses-they can exhibit linear entanglement growth and correlations heating up to infinite temperature-both their spectrum and their OTOC behavior are anomalous and nonchaotic.
General Clifford circuits have a particularly simple structure to their operator-spreading coefficients: Under time evolution by t steps, a simple Pauli string becomes another Pauli string, where M t is a linear endomorphism acting on the set of strings ðZ ⊗2 q Þ ⊗L . Thus, a Pauli operator evolves to a single Pauli operator, rather than the superposition of Pauli operators allowed by Eq. (2). Here, M t has to obey a number of constraints. In particular, time evolution should preserve the commutation relations amongst the Pauli strings [59][60][61]. (Incidentally, using these constraints, it is possible to classify all translation-invariant Clifford circuits into three types called fractal, glider, and periodic [60,61].) In line with the stringency of the constraint Eq. (31), it is unsurprising that Clifford circuits have pathological properties distinguishing them from more generic ergodic systems. Calculating the OTOC Eq. (16) for initial Pauli strings σ μ , σ ν for a Clifford circuit gives where θ is the phase obtained by commuting σ M t ðμÞ through σ ν , and the final result is independent of the state ψ. For generic Clifford circuits, this result shows persistent Oð1Þ oscillations at late times and does not settle to a specific limit even in the thermodynamic limit. We refer to this as "trivial" OTOC behavior because we do not see the regime of persistent decay and eventual saturation characteristic of "chaotic" quantum systems (e.g., the SYK model or the models of Secs. III B and IV). Additionally, Clifford circuits tend to have pathological spectral properties not associated with ergodic systems. For instance, one can prove that translation-invariant Clifford circuits have exact recurrences U t rec f ∝ 1 on time scales linear in system size t rec ¼ OðLÞ (see Appendix F). This directly implies that the eigenvalues of U t rec f are t rec th roots of unity, which, in particular, implies that the average level degeneracy is Oð2 L =LÞ. This spectral structure does not exhibit the level repulsion we expect in systems with ETH.
Although operators obey the stringent condition Eq. (31), Clifford circuits can still exhibit many ergodic properties usually associated with "ergodicity." Indeed, it can be rigorously proven that the above-mentioned "fractal" Clifford circuits exhibit the following: (a) Linear entanglement growth starting from certain so-called stabilizer initial states [62] and relatedly (b) starting from an initial product state, at long times, all local observables tend towards their infinite-temperature expectation values [63]. These fractal Clifford circuits are periodic in time, in addition to being spatially translation invariant. An explicit example of such a circuit has q ¼ 2 and takes the form Note that the resulting circuit has the geometry illustrated in Fig. 1 (with the one-site rotations merged into the two-site gates), and the circuit elements repeat every two layers.
Operators evolve under this circuit in a particularly simple way: Note that for certain strings, there is a possibility of cancelation, e.g., Y s Z sþ1 → Y s−1 X s 1 sþ1 . Figure 8 shows that for this particular circuit, even more generic initial states show near linear entanglement growth, with a rate that is mostly independent of the initial state at long times. We can explain this using the operator-spreading picture of entanglement growth, discussed in Sec. III D. Looking at Eq. (34), we notice that a string with a Pauli operator X or Y at its right endpoint will keep moving to the right at a fixed speed of 1 site per 1 period forever, leading to the fixed-rate linear entanglement growth seen for stabilizer states. The deviations from this behavior for a random product initial state come from strings that, up to time t, failed to start growing to the right. However, such operators need to have a very specific structure, and their number is exponentially small in t (going roughly as 4 −t ). Consequently, after the first few periods, most product states settle to the same entanglement velocity exhibited by stabilizer states. This is in contrast to the behavior of the In summary, Clifford circuits can be ergodic in the sense that they can exhibit linear entanglement growth and thermalization of local observables but do not have the spectral or OTOC behavior expected of generic ergodic systems. This suggests that linear entanglement growth is a rather coarse measure of quantum information spreading, sensitive only to the fact that operators tend to grow in extent over time. The OTOC, on the other hand, is sensitive to both the fact that operators grow in extent and the fact that they become complicated superpositions of many Pauli strings (i.e., operator entanglement in the sense of Ref. [64] increases). Note also that the entanglement entropy tests the behavior of a large ensemble of different Pauli strings, while the OTOC characterizes the evolution of a single initial Pauli operator, thus giving more detailed information on the dynamics.

VI. CONCLUSIONS
We considered the spreading of quantum information in one-dimensional systems with local unitary time evolution but no conservation laws. Our key results are as follows. For random circuit systems, we derived an exact hydrodynamical description of operator spreading. According to this description, operators grow into superpositions of Pauli strings, which tend to be supported at a front, propagating with velocity v B , where v B is a velocity scale distinct from the light-cone velocity (see Fig. 2). An important consequence of the hydrodynamic equation is that the front itself undergoes a diffusive broadening in time. We proved that the velocity v B also determines the characteristic scale of change of the OTOC between two local operators, while at very early times, long before the arrival of the front, we find a regime of exponential growth for the OTOC, with an exponent that depends on the initial separation of the two operators-it remains to be determined whether this early-time exponential growth represents a quantum analogue of the Lyapunov behavior present in classical chaotic systems. The exact description of operator spreading in this model also allows us to give an exact formula for evolution of entanglement across a cut as measured by S ð2Þ exp , starting from an initial product state. We find that the entanglement grows with a third distinct velocity scale, the entanglement velocity v E , which is strongly affected by the diffusive behavior, making v E < v B .
Comparing our exact results for the random circuit to an ergodic Floquet spin chain, we verified the presence of a diffusively broadening operator front in the kicked Ising model. This leads us to a tentative conjecture that this behavior is present for generic ergodic interacting quantum systems in 1D, evolving under local unitary Floquet evolution, allowing for a universal coarse-grained hydrodynamic description in these systems. We contrasted this with the fine-tuned behavior observed in Clifford circuits, for which operators can spread ballistically but do not become entangled superpositions of many operators. The ballistic spread of operators leads to the linear entanglement growth and the thermalization of local observables seen in certain Clifford circuits, but the fine-tuned nature of the operator growth shows up in the pathological behavior of OTOCs. This demonstrates that linear entanglement growth and thermalization are not good predictors of the scrambling behavior captured by the OTOC.
We note that within existing field-theoretic calculations [19][20][21], the broadening of the operator wavefront is not present, in contrast with our results. This discrepancy could be connected to the weak-coupling limits implicit in these works and the unbounded local Hilbert spaces associated with continuum field theories; in any case, the discrepancy should be resolved. Another possible direction for future work is finding the exact range of validity of the hydrodynamic operator-spreading picture proposed here. It is known, for example, that the propagation of the OTOC becomes slower than ballistic in strongly disordered (but still thermalizing) systems [65], and it is an interesting question whether this effect can be incorporated into some form of hydrodynamic description.  (33). The plot shows the growth of entanglement for seven different initial random product states. After the first few steps, the growth rate of entanglement becomes roughly the same for all initial states, set by the fact that the majority of operators travel with a fixed velocity of 1 site per 1 period.
Note added.-Recently, we became aware of a related work [66]. Where they overlap, our results appear to be consistent with this work. Additional numerical work also appeared recently, apparently confirming the presence of diffusively broadening fronts in Hamiltonian systems [29].
APPENDIX A: VALIDITY OF EQ. (14) It suffices to work in coarse-grained coordinates t, x. Recall thatR Our task is to justify Eq. (14), which can be more carefully phrased as follows: In the limit t → ∞, with jx − v B tj=t ¼ ε a fixed nonzero number, the integrated operator density obeysR where c 0;1 are positive numbers, bounded in the t → ∞ limit, and we defined κ ≡ x=t and work in units The quotient Q x is always positive. It is easy to verify that for −1 < κ < v B , the quotient is greater than 1 and is an increasing function of x. On the other hand, for v B < κ < 1, this quotient is less than 1 and a decreasing function of x. When −1 < κ < v B , we use these facts to bound On the other hand, we immediately see thatRðxÞ ≥ ρ R ðx; tÞ. Noting that in the large t limit Hence, c 0 ≤1þ½ð1þκÞð1−v B Þ=2ε is an Oð1Þ constant. Similarly, we consider v B < κ < 1. Using 1 ¼ P y ρ R ðy;tÞ, it follows thatRðxÞ ¼ 1 − P y>x ρ R ðy; tÞ. In the present case, Q x < 1, and it is straightforward to derive a similar bound, Hence, is an Oð1Þ number in the large t limit. Note that near the edge of the light cone, κ ¼ ∓1 respectively, Eq. (A3) becomes increasingly exact as each of the bounds Eqs. (A8) and (A9) becomes tighter.
On the other hand, as we approach the front ε → 0, the bounds become looser and Eq. (A3) is less reliable-in this regime, the near front expansion Eq. (12) becomes more useful.

APPENDIX B: APPLICATION TO OTOC
Following Sec. III B, the quantity fðs; τÞ ≡ 1 −Cðs; τÞ can be written as where e iθ μ;Z s is a qth root of unity acquired by commuting σ μ past Z s . We can reparametrize the sum by summing over the left and right endpoints of the Pauli string μ, There are also contributions to Eq. (B2) that arise when s is on the left and/or right edge of an interval, i.e., s ¼ l or s ¼ r.

APPENDIX C: EXACT OPERATOR-SPREAD COEFFICIENTS
In what follows, we give an exact expression for the averaged operator-spread coefficients and a sketch of the derivation. We leave a more detailed derivation to future work. The averaged operator-spread coefficients are defined as The Haar averaging can be performed explicitly using the identity (E1) in Appendix E. After averaging, each two-site gate can be represented by a classical, Ising-like variable taking only two possible values. Because of the geometry of the circuit, these Ising variables form a triangular lattice. Equation (C1) becomes a classical partition function, i.e., a sum over all possible spin configurations. The partition function consists of two edge parts, which depend on the configurations on the first (last) layer and the Pauli strings ν (μ), and a bulk part, which is independent of the Pauli strings in question. Because of the Haar averaging, the only information that remains in the partition function about the strings μ and ν is which sites they act on nontrivially. The bulk transfer matrix enforces a light-cone structure on the spin variables. A light cone with velocity v LC ¼ 1 emanates from the two endpoints of the string ν such that all spins outside of the light cone have to point up [otherwise, the configuration has zero weight in the partition function for jc ν μ ðtÞj 2 ]. In the case of an initial Pauli operator acting on a single site, the partition function for the operator spreading can be evaluated exactly. In this case, the fact that μ acts on one site only yields a boundary condition for the partition function, wherein, in the first row, there is a single spin pointing down while all others point up. The partition function then becomes a sum over all possible ways this initially one-site domain can spread within the light cone, as shown in Fig. 9. Furthermore, the bulk interaction terms are only nontrivial at the boundary between the two domains; consequently, they give the same contribution for all domain configurations with the same depth. Thus, the calculation simplifies to counting the possible domain configurations, which can be done by considering it as a two-particle random walk for the two endpoints of the domain.
FIG. 9. Example of a classical spin configuration contributing to the operator-spreading coefficients of an initial one-site Pauli operator. All such configurations have a single domain of down spins spreading inside the light cone, and each of these configurations contributes equally to jc ν ðtÞj 2 . The calculation outlined above yields all average squared coefficients jc ν μ ðtÞj 2 , where μ j ¼ 0 if j ≠ 0. In the following, we simplify the notation by dropping the first index and denoting these as jc ν ðtÞj 2 . The exact formula for these squared coefficients is given in Eq. (C2). A surprising property of this formula only depends on the positions of the right and left endpoints l, r of the Pauli string ν, and not on more detailed information concerning the internal structure of the string, i.e., jc ν ðtÞj 2 ¼ hðl; rÞ. This is a consequence of the Haar averaging, which, in each step, washes out all memory of the internal structure.
We form a circuit with an even number of D layers. Then, we number the spaces between the two-site unitaries in the last layer, −D=2; 1 − D=2; …D=2. The support of an operator string ν can be represented by x; y ∈ f−D=2; 1 − D=2; … þ D=2g, with x < y. The average square coefficient is obtained by plugging x, y into the formula.
Note that this expression depends only on x, y and q, so we also denote it as jc μ 1 ν j 2 ¼ hðx; y; DÞ, where we drop the q dependence for simplicity. These expressions are of course complicated. We note as an aside that this formula has a slightly neater expression in terms of hypergeometric functions.

Useful limits
Let us calculate the weight on an operator with endpoints x, y in the large D limit (in circuit coordinates) starting from Eq. (C2). Reexpress the weight as

APPENDIX D: LONG-TIME CORRELATIONS
In Sec. I, we anticipated that random unitary circuits should "heat to infinite temperature," much like ergodic Floquet circuits. Here, we back up this claim by examining the long-time behavior of various correlation functions, demonstrating that they relax to their expected infinitetemperature values. For simplicity, we consider one-and two-point functions for specific one-site Pauli operators of the form σ α 0 -although most of the results below follow for more general operators as well. Fix any initial state ω and times t 2 , t 1 . We find that hσ α 0 ðt 2 Þσ β z ðt 1 Þi ω ¼ 0: This follows from two observations. First, the operator in the expectation value can be written U 01 U 12 σ α 0 U 21 σ β x U 10 , where U ij ¼ U −1 ji is shorthand for the unitary evolutions