Coarse-grained Entanglement and Operator Growth in Anomalous Dynamics

In two-dimensional Floquet systems, many-body localized dynamics in the bulk may give rise to a chaotic evolution at the one-dimensional edges that is characterized by a nonzero chiral topological index. Such anomalous dynamics is qualitatively different from local-Hamiltonian evolution. Here we show how the presence of a nonzero index affects entanglement generation and the spreading of local operators, focusing on the coarse-grained description of generic systems. We tackle this problem by analyzing exactly solvable models of random quantum cellular automata (QCA) which generalize random circuits. We find that a nonzero index leads to asymmetric butterfly velocities with different diffusive broadening of the light cones, and to a modification of the order relations between the butterfly and entanglement velocities. We propose that these results can be understood via a generalization of the recently-introduced entanglement membrane theory, by allowing for a spacetime entropy current, which in the case of a generic QCA is fixed by the index. We work out the implications of this current on the entanglement"membrane tension"and show that the results for random QCA are recovered by identifying the topological index with a background velocity for the coarse-grained entanglement dynamics.

In one-dimension (1D), the mathematical theory of LP dynamics is well developed [9][10][11][12][13][14][15][16].A crucial result, proven in Ref. [10], is that QCA are fully classified by a genuinely dynamical topological index.This result was recently generalized to include the more realistic case where the causal cone is only approximate [17].Importantly, this topological index is zero if and only if the evolution is generated by a (quasi-)local Hamiltonian.Otherwise, the dynamics is said to be anomalous.As a natural application, this theory led to the discovery of new dynamical topological phases in 2D Floquet MBL systems [4][5][6][7][8], including the case of protecting symmetries [18][19][20], which go beyond the cohomology paradigm [21][22][23].
Although the index was initially defined in terms of abstract operator algebras [10], an equivalent definition, which reflects an intuitive picture of quantum-information flow [7,16,17], was recently put forward.In turn, this made it possible to establish a lower bound on quantum scrambling in terms of the index, building a bridge between genuinely dynamical topological invariants and quantum chaos [16].
In this Letter, we develop a connection of a different kind: while Ref. [16] derived universal relations involving the index at the microscopic level, here we reveal its implications for the dynamics of generic systems at macroscopic (hydrodynamic) In the presence of an MBL phase in the bulk [27], it is always possible to decompose the one-period Floquet operator as [4] UF = U edge e −iH bulk T , where H bulk is an MBL Hamiltonian, and U edge is an effective 1D evolution acting on qudits within a few localization lengths from the boundary.(b): Any LP evolution U(t) may be approximated by a QCA in Margolous form, i.e. the single time-step unitary operator U admits a bilayer representation, where the local unitaries map a product of nearestneighbor Hilbert spaces into another, with possibly different individual input or output dimensions.
Following the logic of these works, where random unitary circuit (RUC) models played a key role, our approach is based on the analysis of random QCA, which we propose as minimal models for anomalous chaotic systems.Conventional EMT.-Let us begin by recalling the basic aspects of the EMT [24,25].Throughout this work, we consider a 1D lattice of qudits associated with a Hilbert space (C d ) ⊗2L (2L: system size) and a unitary dynamics dictated by the operator U(t) : (C d ) ⊗2L → (C d ) ⊗2L , where time t might be either continuous or discrete.
The main object of the EMT is the so-called membrane tension (or line tension, in 1D), which associates an entanglement cost with a given spacetime cut through the unitary arXiv:2109.07408v2[cond-mat.stat-mech]24 Feb 2022 u v operator U(t) [cf.Fig. 2].This quantity allows for an intuitive geometric picture for the coarse-grained entanglement dynamics.The local tension E(v) is a function of the curve velocity v = dx/dt, and the cost of a given curve is obtained by integrating E(v) along its length.Then, the entanglement of a given interval A in space at a given time is obtained by minimizing the integral of E(v) over all curves that separate a spacetime region that terminates on A on the temporal boundary.As an example, we may consider the growth of the entanglement after a quench, for an infinite bipartite system with open boundary conditions: assuming homogeneous spacetime dynamics, we obtain S(x, t) = min y ts eq E x − y t + S(y, 0) , where S(y, 0) is the entanglement of the initial state, while s eq is the entanglement density reached at equilibrium [25].S(x, t) here may indicate the von Neumann entanglement entropy, or (assuming the absence of conservation laws [28][29][30][31]) an arbitrary Rényi entropy [32].Holographic field theories give elegant examples of off-lattice systems where E(v) is explicitly computable [33][34][35].
The EMT may be equivalently formulated in terms of a local entanglement production rate.In the bipartite setting above, the membrane picture is equivalent to a dynamical equation ∂S ∂t = s eq Γ ∂S ∂x , where Γ(s) is a local production rate dependent on the entanglement gradient [24,25].Comparison with (1) reveals that Γ(s) and E(v) are simply related by the Legendre transformation The line tension also encodes information about the growth of local operators and, in general, must satisfy some basic constraints [25].First, internal consistency of the coarse-grained picture requires that E(v) ≥ 0 and E (v) ≥ 0. Second, one may argue that the minimization in Eq. ( 2) only involves membrane velocities within a range [−v − , v + ], where v ± coincide with the left/right butterfly speeds v L,R that govern the growth of local operators.Although v L = v R if spatial inversion symmetry is present, this is not generally true otherwise, even for local-Hamiltonian dynamics [36][37][38].The minimum of E(v) is the entanglement velocity v E , quantifying the entanglement growth rate (rescaled by s eq ) starting from a product state.Finally, one can show and This picture is believed to hold for generic local-Hamiltonian and quantum-circuit evolution; our goal is to find whether and how it can be extended to anomalous dynamics.The Margolus form for QCA.-As mentioned, the edge dynamics of the 2D systems in Fig. 1(a) is LP.This means that the single time-step unitary operator U of the discrete evolution has the following property: for any local observable O j acting on site j, the operator U † O j U is supported on a finiteneighborhood of j, up to exponentially decaying tails [39].
Any LP dynamics may be approximated arbitrarily well by QCA, by "chopping off" the exponential tails [17].In turn, it is known that any QCA may be expressed in the so-called Margolus form [40] [cf.Fig. 1(b)], where U is written as a two-layer product of two-site unitaries.This does not always define a quantum circuit because the dimensions of the local spaces associated with the "virtual" layer may differ from the physical ones: denoting by p, q two integers such that d 2 = pq, we have where u : Given this representation, the topological index reads [10] ind = 1 2 ln q p . ( The unitary operation of translation by one site is a simple example with ind = ln d.Note that for finite systems ind = 0 is only possible for periodic boundary conditions [41].For simplicity, we will always take L → ∞, so that the boundary conditions become irrelevant.The Margolus representation allows us to pinpoint the essential feature of anomalous dynamics which we have to take into account in order to generalize the conventional EMT.First, we note that the QCA in Fig. 1(b) can be viewed as a unitary tensor network (TN) and that, although the dimensions associated with given bonds may vary in space and time, unitarity requires that the input and output dimensions of each tensor must match.This gives rise to a non-trivial local conservation law, not accounted for in conventional EMT.
Physically, we can understand such local conservation law as a continuity equation of the form ∂ µ J µ = 0, in terms of a coarse-grained spacetime entropy current J µ .For a unitary TN dynamics locally equilibrating to infinite temperature, J µ has an explicit microscopic definition: Regarding the TN as a graph whose nodes are the unitaries and edges the bonds, we orient the latter in the direction of increasing time; then, along each bond we define the entropy current as a vector in 3. In the stationary regime, the divergence of the entropy current is vanishing, and the flux through any closed surface is zero.In order to compute the integral (6), we choose a surface containing the optimal membranes and exploit the prescriptions of the EMT.
the direction of its orientation, whose magnitude is equal to ln d i , with d i the associated local Hilbert-space dimension.
The coarse-grained spacetime entropy current is more general than the model above, where it can be introduced by "counting" microscopic bonds.For instance, it can be defined even when the equilibrium state is non-trivial and determined by the slow modes [42].In any case, it has important consequences on the properties of the membrane tension.In the following, we show how the EMT has to be modified in the presence of a non-trivial spacetime entropy current.The resulting generalized EMT turns out to correctly capture the coarsegrained features of anomalous dynamics, finally revealing the hydrodynamic implications of the index.
Generalized EMT.-As a starting point to generalize the EMT, we assume there is a well-defined line-tension E(v) satisfying E(v), E (v) ≥ 0 (as required for consistency of the hydrodynamic picture).We also postulate that there exists a spacetime entropy current J µ (x, t) = (J x , J t ) which governs the growth and transport of thermodynamic entropy.In particular, the density s of thermodynamic entropy (we assume local equilibrium) is equal to J t .
Let us consider a stationary regime, focusing, for instance, on a finite interval A at large times after a quench.Stationarity requires ∂ µ J µ = 0, and also implies that the thermodynamic entropy equals the von Neumann entropy, S A (t), i.e.
Using the divergence theorem, the integral over A may be obtained by integrating the current over any closed perimeter containing A [cf. Fig. 3].In order to make contact with the EMT, we choose the perimeter to be a triangle whose bottom sides have slopes given by the butterfly velocities −v L , v R .A simple computation then yields S A (t) = (J x ∆t + J t a) + (−J x ∆t + J t b), where ∆t, a and b are as in Fig. 3. On the other hand, the sides of triangle minimize the line tension for the region A, as one can show by generalizing the arguments of [25].As a consequence Identifying the individual terms coming from the two bottom sides of the triangle, and using a/∆t = v L and b/∆t = v R , we find . Defining now the background entropy velocity v * = J x /J t , and using S A (t) = s eq (a + b) (which follows from stationarity) we finally obtain This equation deviates from the conventional EMT, cf.(3), and has important ramifications.In particular, combined with convexity, it implies Using (7), one can also argue that the relation between Γ(s) and E(v) must be modified.Indeed, plugging (7) directly into (2), we see that Γ(s) < 0 for some values of s.This is clearly an inconsistency, since Γ(s) is the rate of entanglement growth.In order to guarantee positivity, one is led to the natural generalization From elementary properties of the Legendre transformation, we have the basic constraints −s eq Γ (s eq ) = v R − v * , s eq Γ (−s eq ) = v L + v * .Note that by construction Γ(0) = v E .Finally, differentiating (1) with respect to t, we obtain We see that the fundamental equation from the EMT, governing the coarse-grained entanglement dynamics, is modified by a constant velocity term.Importantly, v * is now left as a free parameter, and the conventional EMT is recovered for v * = 0. Note that without entropy production, i.e.Γ ∂S ∂x = 0, Eq. ( 9) still predicts a nonzero entropy change, which is qualitatively different from normal dynamics.
Models of random QCA.-In order to test the generalized EMT and identify the entropy-current velocity v * , we study concrete models of chaotic anomalous dynamics.We consider QCA of the form (4) where u, v at different spacetime positions are drawn independently from the Haar random ensemble.Generalizing from the special case of RUCs, we expect the model also to capture universal aspects of random Floquet evolutions [43][44][45][46][47][48][49] and translationally invariant homogeneous systems [24][25][26]50] if we restrict to the leading dynamics at large scales.We note that it is easy to construct explicit 2D models with trivial bulk dynamics that display the random edge evolution considered here.This construction is detailed in the Supplemental Material (SM), where we also define even simpler random 1D QCA which appear naturally in this context [51].
As a first step, we analyze how the support of a localized traceless operator O 0 grows under the dynamics, which allows us to extract the butterfly velocities.We focus on the out-of-time-order correlator (OTOC) and O 2x is a traceless operator supported at site 2x.Given the brickwork structure of the random QCA, A r e r C f r x X q 3 P q a t B W s 2 s w / + l P X 5 A 3 O D m z c = < / l a t e x i t > ind s eq < l a t e x i t s h a 1 _ b a s e 6 4 = " e b O N D B q H 0 M m 0 J y 9 t c N 3 a 5 E 0 q E 0 4 7 k k I p e l i p / y u Q h y + y R 7 5 T c i j s t u A h e D i W Q V 8 1 3 v j q 9 C C e c C I 0 Z U q r t u b H u p k h q i h n J 7 E 6 i S I z w E P < l a t e x i t s h a 1 _ b a s e 6 4 = " x y T P S 9 s y 5 d / 9 q G y 4 P w m D p w N w 4 c A 0 c e W A k / Y A 7 V 7 J 1 s / P + 7 y y / Y Q J Z 9 p 4 3 n d n b / / R 4 y d P a 8 / c 5 y 8 O X h 7 W j 1 4 N d Z w q Q g c k 5 r E a h 6 A p Z 5 I O D D O c j h N F Q Y S c j s L l h 7 I / y q j S L J Z f z C q h M w E L y S J G w F g p q F 8 2 s + C z a 7 e P b n M a K S D 5 c 5 k 0 l q q C T r i 6 K U Y x P j M k 8 8 Z 4 o S w 1 c W g C h m 3 4 r J F d g g j U 3 d t S H 4 2 1 / e h e F x x z / t e J 9 O G t 2 z T R w 1 9 A a 9 R S 3 k o 3 e o i y 5 Q H w 0 Q Q d f o B / q F 7 p w b 5 6 d z 7 / x e W / e c z c x r 9 F 8 5 f / 4 C j H a t L g = = < / l a t e x i t > v E < l a t e x i t s h a 1 _ b a s e 6 4 = " the disorder-averaged OTOC C(x, t) may be computed using the approach developed for RUC, mapping the problem to the partition function of an Ising-like model [57,58], see also [59][60][61][62][63][64].Additional technical complications arise due to the "staggered" structure of the dynamics, alternating physical and virtual Hilbert spaces [cf.Fig. 1(b)].Nevertheless, a fully analytic expression may be obtained [51], and in the hydrodynamic limit of large spacetime scales, it simply reads C(x, t) while v R and σ R are obtained by exchanging p ↔ q.We see that the coarse-grained OTOC has the same form as RUC [57,58], being characterized by two propagating fronts inside which C(x, t) 1, i.e. information is fully scrambled.However, the fronts propagate with asymmetric butterfly ve-locities and widths, with the faster front being the narrowest one (in contrast to the quantum circuit models constructed in Ref. [37]).Note that when ind = 0, i.e. p = q, we recover the result for RUC [57,58].In the other limit where ind = ln d, i.e. p = 1 and q = d 2 , the random QCA consists of a two-step evolution in which a right translation is followed by a layer of random unitary gates.Since the shift does not increase the operator support, σ L,R in (10) are then simply those of a RUC evolved up to a time t = t/2.Entanglement dynamics.-Next,we move on to compute the line-tension, cf.
) |I ∈ H, with the two Hilbert spaces associated with its input and output degrees of freedom.Considering a bipartition of the system with boundary at site x and y in the input and output, respectively, we define the operator entanglement [65-68] S (o) (y − x, t) as the associated entanglement entropy of |U(t) .The line-tension may then be computed via [25] Our random QCA has no conserved quantity and our unit length increment dx = 1 is defined to contain two sites, so that s eq = ln is the operator-entanglement velocity.
Averaged Rényi-n entropies are hard to compute, but the mapping to an Ising partition function gives access to the averaged purity and its logarithm − ln e −S (o) 2 (x,t) .Since this average is taken "inside the logarithm", it differs from the averaged Rényi-2 entropy S (o) 2 .However, the former is sufficient to see the key relationships obeyed by E. The line tension for S (o) 2 can be understood as a perturbatively "dressed" version of that for − ln e −S (o) 2 (x,t) [26]: the difference between the two vanishes as d → ∞, and in fact is numerically small even for finite d [69].So we approximate t where r p = p(q 2 − 1)/[q(p 2 − 1)].As a main difference from the case of RUC, the minimum is at a nonzero veloc- as confirmed by directly computing the growth of state Rényi-2 entropy following a quench from a product state [51].We stress that v L,R and v E do not depend on ind in a universal way.This could be expected from the study of RUC, where asymmetric butterfly velocities might be realized by specific arrangements of the local unitaries [37].
Crucially, we see that E 2 (±v ± ) = ±1, where v ± = v R,L , and that Eq. ( 7) is satisfied, after the identification This is our final main result: it states that the index, a microscopic dynamical topological invariant, appears at the hydrodynamic level as a constant background velocity for the coarse-grained entanglement dynamics.Based on this identification, the index also determines the qualitative features of the rate Γ(s), which is shown in Fig. 4. Further details on the random QCA, including a computation of the so-called tripartite mutual information [70] and its relation to the index, are reported in the SM [51]. Outlook.
-Our results open up several possibilities for future research.First, when viewing anomalous 1D dynamics as boundaries of 2D Floquet systems, it would be interesting to investigate the corrections to our theory when the assumption of ideal localization in the bulk is relaxed.In this case, we expect subleading effects emerging, due to a slow entropy flow from the boundary to the bulk, and vice versa.It would also be natural to apply our picture based on spacetime entropy currents to more general situations with inhomogeneous backgrounds, as models with genuinely spacetime-dependent entropy currents may be constructed by introducing additional structure.Next, it would be interesting to explore how chaotic anomalous dynamics is modified by local conservation laws, such as U (1) charges as done for RUC [59,60].Studies along this direction could reveal an intriguing effect of the index on the otherwise purely diffusive behavior of the charge.Finally, two natural generalizations of our study include adding randomized measurements [71][72][73][74][75][76][77][78] and higher dimensions, where the theory of QCA is much more open [79][80][81][82].

Supplemental Materials
We provide the detailed derivations of various dynamical quantities for the exactly solvable random QCA in the main text.We also discuss a different simplified random QCA and the parent Floquet systems of both.

MODELS OF RANDOM QCA
In the main text, we have considered a natural model of random QCA.It is defined on a chain with periodic boundary conditions, associated with the Hilbert space (C d ) ⊗2L (2L: system size).The unitary dynamics is discrete, and dictated by the evolution operator Here u : unitary operators, where p and q integers such that pq = d 2 .In our model, u and v are drawn randomly from the Haar distribution and independently for each point in space and time.If p = d, then the dynamics is anomalous, and is characterized by a nonzero index [10].Here, we omit its mathematical definition and simply recall that for the time-step evolution operator U s in (S1) it is easily computed and reads In fact, it is possible to construct other simple models of anomalous dynamics.Indeed it is known that any discrete unitary evolution with a given index may be obtained by combining layers of quantum circuits and shifts [10].When d = pq, where p and q are positive integers, the local Hilbert space factorizes as C d = C p ⊗ C q , and one may define two shift operators T p and T q translating to the right only one of the two types of local subspaces.Then, the operator T p,q = T −1 q T p has nonzero index ind = ln(p/q) [10].The dynamics generated by T p,q is clearly not generic.However, in order to obtain a simple chaotic model we may alternate this shift with a direct product of on-site Haar random unitaries u, yielding where u are drawn randomly from the Haar distribution and independently for each point in space and time.U(t) then admits the graphical representation where the dimension of thick (thin) legs is p (q).We see that blocks with even s + j are fully decoupled from those with odd s + j.Also, due to the properties of the Haar measure, we have at the level of random ensembles.Therefore, this construction leads to two identical copies of a simplified model of random QCA, reported in Fig. S1b), with index We will call this Model B, in order to distinguish it from the one defined by (S1) [and depicted in Fig. S1a)], which we will call Model A.
In order to verify the generality of the proposed EMT, we have performed analytic computations in both models.In fact, the calculations are very similar in the two cases.Therefore, we will detail them only for Model A, for which they are slightly more complicated, while for Model B we will merely report the final results.
p is a Haar-random bipartite unitary acting on two adjacent virtual sites.The independent for different [a, a + 1] and s.Schematically, we may represent Eq. (1) as where a thick (thin) leg is associated with C p (C q ).This model is a natural generalization of the Haar-random q [1,2] to QCAs with nontrivial topological indices ind = log p q [3,4].A central quantity we focus on in this note is the infinite-temperature OTOC: where where the dark unitary is the complex conjugate of the light one.For simplicity, we will drop the subscript in E Following the notations in Ref. [5], we may write down the explicit expression as .
Since U † t is equivalent to U t after taking the disorder average, by introducing F (x, t) = (pq) 2t−1 F (x, t), we grap 2 aphically represented as ults comparable with the related models (especially the simplified QCA), we bunch two adjacent sites 2j and gle one j.In this case, s eq = 2 ln d = ln(pq).ute the Rényi-2 operator entanglement using the standard method, i.e., mapping the disordered averaged exprestition of an Ising model with (anisotropic) three-body interactions.The configurations with nonzero contributions domain wall like v u wall contributions are given by e operator entanglement can be expressed as gest horizontal (i.e., u-axis in Eq. ( 7)) coordinate of the domain wall and a counts the total number of horizontal teps.Introducing Eq. ( 10) into to figure out the asymptotic expression for the operator entanglement, which determines the entanglement line e x, t with O(1) ↵ = x t , f (x, t, r) (11) has the following large deviation property (see Appendix A): 2 ) where u [a,a+1] s : C p ⊗ C q → C q ⊗ C p are unitaries drawn randomly from the Haar distribution, and independently for each point in space and time.For p = q the dynamics is anomalous, with index (S6).

DETAILS ON THE PARENT FLOQUET SYSTEMS
Here we briefly discuss the 2D parent Floquet dynamics for the random QCA studied in this work.They can be easily defined for both models A and B.
We begin by analyzing Model B, for which the construction can be straightforwardly obtained from the results presented in Ref. [4].We consider a rectangular system, in which we place a lattice Λ, rotated 45 degrees with respect to a horizontal axis.The local Hilbert space is C d = C p ⊗ C q , each factor forming a sublattice denoted by Λ p and Λ q , respectively.As shown in [4], for Λ p (Λ q ) it is easy to construct a quantum circuit of depth 4 made of swaps which implements a shift of a full period on each plaquette [cf.Fig. S2(b)].In turn, such quantum circuit can be clearly defined in terms of finite-time evolution of a local Floquet Hamiltonian.The ensuing dynamics is trivial in the bulk, while it implements a non-trivial shift at the boundaries.Choosing opposite directions for the shifts in Λ p and Λ q , the (bottom) boundary Floquet operator reads U F = T p,q = T −1 q T p , where T p and T q are translations to the right acting on the two types of local subspaces.We see that in order to obtain (S3), it is then enough to apply a unitary V = r u r , where u r is a random one-site operator acting on the local space C d at position r.Clearly, V may be obtained as finite-time evolution of a random local Hamiltonian.Furthermore, V does not spoil localization in the bulk.This completes the construction of the parent Floquet system for Model B.
Similarly, we can construct a parent Floquet evolution for Model A. This is done following a construction presented in [18], which we briefly sketch here for completeness.We again consider a square lattice with local Hilbert space C d .As shown in the left panel in Fig. S2, we perform the four-step-swap operation on the virtual level, which is pulled back to the physical level by a conjugation by disjoint random unitaries u : C d ⊗ C d → C p ⊗ C q .According to Ref. [18], if this is followed by a four-step-swap operation on the physical level, as is shown in the right panel, the edge dynamics will be a QCA in the Margolous form with v [2j+1,2j+2] = u [2j,2j+1] † S p,q , where S p,q swaps C q and C p .To make v's completely random, i.e., independent of u's, we only have to apply a layer of random two-site gates on the physical level, as shown in the middle panel in Fig. S2(a).Again, these disjoint local unitaries do not spoil the many-body localizability in the bulk.This completes the construction of the parent Floquet system for Model A.

MAPPING TO AN ISING PARTITION FUNCTION: THE OTOCS
Here we provide details on the calculations of the OTOC where O 0 , O 2x are arbitrary traceless operator satisfying . We used the same approach developed in the case of RUC [24,57,58], which is based on a mapping onto the partition function of a classical Ising-like model in 2D.In fact, our calculations follow closely those presented in these works.For this reason, we simply sketch the main steps and work out the parts of the analysis which are different from the case of RUC.We will focus on Model A, and only present the final results for Model B at the end.The first step is to represent the quantities of interest in terms of the "replica operator" ⊗N , where U * (t) denotes complex conjugation with respect to the computation basis.Clearly, U (N ) (t) has the same geometrical structure of U(t), but in a "replica space" where the local Hilbert-space dimension is d 2N .For instance, in the case U(t) is a RUC with brickwork structure, so is U (N ) (t), with local two-site gates u (N ) = (u ⊗ u * ) ⊗N .This representation is very convenient because physical random gates (9) physical right translation the unitaries u (N ) are distributed independently from one another, and averages over disorder then factorize [24,57,58].After averaging, one is left with a non-unitary 1D evolution, which in turn can be mapped onto a classical statistical-mechanics model of Ising spins in 2D.The quantity under investigation determines the boundary conditions of this partition function.The latter may finally be evaluated exactly by solving a problem of domain-wall counting.We refer to Ref. [57] for the details.

Model A
For the random QCA studied here, the very same approach can be applied.As the main difference, the QCA display staggered dynamics, which is reflected in the fact that the weights in the 2D Ising model are different for even and odd rows.Although this makes the enumeration of domain walls more involved, the procedure is straightforward, and therefore we omit it.As a final result we obtain the following exact expression, valid at finite distances and times where and As a first check, we see that for p = q = d we recover the result for the OTOC in RUC derived in Ref. [57], after replacing (x, t) = ( x 2 , t 2 ): this can be seen using Vandermonde's identity, which gives Let us now consider the OTOC in the limit of large spacetime distances.First, neglecting subleading terms, we may rewrite (S8) as The analysis of (S13) for large x and t is significantly more complicated than the case of RUC.This is because for p = q the functions f t−1 L (x), f t−1 R (x) can not be expressed as simple binomial coefficients.By inspection, we see that the hydrodynamic limit of C(x, t) is determined by the terms in the sums (S9), (S10) where u and n are of the same order.Hence, we need an asymptotic expansion for f 2t−1 L/R (u) when u = αt and t → ∞.A systematic analysis of this limit may be obtained by expressing f n L/R (x) in terms of the Gaussian hypergeometric function and exploiting known asymptotic expansion formulas [83,84].Since we are only interested in the leading behavior in t, we follow a simpler and more direct approach.We begin by considering the function which can be rewritten as where C z can be an arbitrary closed loop encircling the origin.Next, we parametrize z = ρe ik , with ρ, k ∈ R, and setting x = αt, we obtain For large t, we wish to evaluate this integral via the saddle-point method, which yields the saddle-point equation namely For r, α fixed, there is only one ρ * > 0 for which the solution to (S18) satisfies k ∈ [0, 2π].Therefore, we may choose ρ = ρ * and apply the saddle-point method.By examining the limit of r = 1, we know that the dominant contribution should arise choosing the "+" sign in (S18).Substituting the corresponding root into Eq.(S16) yields which finally determines the leading exponential behavior of the function (S14).Next, to the leading order in t, it is immediate to see that f 2t−1 L/R (αt) f (t, αt, r ±1 p ), where Therefore, we may exploit the asymptotic expansion (S19) and obtain the final result where A(t, α) is a function growing sub-exponentially in t, and Finally, let us consider the sums in Eqs.(S9), (S10).They are dominated by the terms where exp[tH(α, r ±1 p )]d −2u is largest, the others being negligibly small.This happens for values of u that are of the same order of t, namely u = αt.The leading terms are thus once again determined via the saddle-point condition which yields Recalling (S13), this gives us the left/right velocities Expanding up to the second order in α around the saddle points α * ± , we can also obtain the diffusive broadening of the light cones.Defining we get where .

(S29)
Assuming that near α * ± we can neglect the contributions coming from the subleading function A(t, α) (which varies slowly in t), we have that, in the region where g L (n, a) is not negligible, it is a discrete sum of Gaussian weights.Therefore, we can take the continuum limit, and obtain the final result where We have tested numerically this result, see Fig. S3 for an example.As an analytic check, we have verified that for p = q = d we obtain the same formula for σ in RUC obtained in [57], after replacing (x, t) with ( x 2 , t 2 ).It is interesting to consider the the limit p → 1.In this case, setting q = d 2 we find which is exactly the expression for a quantum circuit evolved up to time t = t/2.This is consistent, because a simple shift does not increase the operator support, and hence does not modify the broadening of the light cones.

Model B
Similar calculations may be performed for Model B. In fact, in this case the analysis is significantly simpler, and here only report the final result.Defining the OTOC for Model B reads where In the hydrodynamic limit, C(x, t) is still of the form (S30) where now Note that when p = 1 (q = 1), we have ) and σ L (t) = σ R (t) = 0, which is fully consistent with the left (right) shift.Furthermore, when p = q, we again reproduce the results in Ref. [57] after replacing (x, t) with ( x 2 , t 2 ).

THE OPERATOR ENTANGLEMENT AND THE TRIPARTITE INFORMATION
In this section we discuss different entanglement-related quantities.To this end, we define a doubled Hilbert space H = (C d ) ⊗2L ⊗ (C d ) ⊗2L along with the maximally entangled state |I = |φ + ⊗2L , where |φ + ≡ d n=1 |n ⊗ |n / √ d and n runs over a basis of C d .This allows us to vectorize the unitary operator as the so-called Choi state |U(t) = (U(t) ⊗ 1 1) |I ∈ H, with the two Hilbert spaces naturally associated with its input and output degrees of freedom.All the entanglement-quantities can be obtained by calculating the entropies of the reduced states on properly chosen subsystems.

Model A
Considering a bipartition of the system with boundary at site x and y in the input and output, respectively, our goal is to compute the bipartite entanglement S (o) (x, y, t) of |U(t) for an infinite system (L → ∞).Exploiting translation symmetry, we have S (o) (x, y, t) = S (o) (x−y, t).As mentioned, although computing Rényi-n entropies exactly is very challenging, our model makes it possible to obtain an analytic expression for the average of the purity, i.e. e −S (o) 2 (x,t) .Indeed, this quantity may be expressed in terms of two copies of U(t) ⊗ U * (t), and can be computed exploiting the very same mapping to an Isinig partition function in 2D used for the OTOCs [24,57].As before, we omit the details of this procedure which follow closely those for RUC.We obtain the following exact expression valid at finite x and t e which may be simplified as where f (t, x, r) and r p are defined in (S14), (S20) respectively, and we used the identity In order to extract information on the line tension, we need to analyze the hydrodynamic limit of (S39).This follows once again from (S19).Indeed, we first note that, within the light cone, we have the following leading behavior Therefore, the Rényi-2 entanglement line tension ts eq = log pq (pq + 1) The main qualitative features of E 2 (v) can be seen from Fig. S4, where we have plotted it for given values of d, p and q.
Eq. (S43) allows us to extract several quantities and check explicitly a few aspects of the generalized EMT.First al all, the value of E 2 (v) in v = 0 defines the operator Rényi-2 entanglement velocity Next, its minimum value determines the entanglement speed v E [24,25].To compute it, we need to solve E 2 (v) = 0.While this is a complicated transcendental equation, we can eliminate the intractable logarithm term by making the ansatz which leads to the solution It is immediate to verify that v m found in this way indeed satisfies (S45).Plugging (S46) into (S43) we obtain the entanglement velocity We can check explicitly that this is the growth rate of the state Rényi-2 entropy after a quench.Indeed, the latter is readily computed for an initial product state, and reads A r e r C f r x X q 3 P q a t B W s 2 s w / + l P X 5 A 3 O D m z c = < / l a t e x i t > ind s eq < l a t e x i t s h a 1 _ b a s e 6 4 = " e b O N D B q H 0 M m 0 J y < l a t e x i t s h a 1 _ b a s e 6 4 = " x y T P S 9 s y 5 d / 9 q G y 4 P w m D p 1

a t e x i t s h a 1 _ b a s e 6 4 = " J i 4 p Y I C j S h s I 3 Z o N t R t b / h x M p w o = " >
c 5 k 0 l q q C T r i 6 K U Y x P j M k 8 8 Z 4 o S w 1 c W g C h m 3 4 r J F d g g j U 3 d t S H 4 2 1 / e h e F x x z / t e J 9 O G t 2 z T R w 1 9 A a 9 R S 3 k o 3 e o i y 5 Q H w 0 Q Q d f o B / q F 7 p w b 5 6 d z 7 / x e W / e c z c x r 9 F 8 5 f / 4 C j H a t L g = = < / l a t e x i t > v E < l a t e x i t s h a 1 _ b a s e 6 4 = " < l a t e x i t s h a 1 _ b a s e 6 4 = " K 5 a 0 0 from which we see that − lim t→∞ ln e −S2(t) /ts eq equals v E in (S47).Finally, we can also easily check that the butterfly velocities (S31) satisfy the identities ) Eq. ( S51) is a natural generalization of E(v) = |v|, ∀|v| > v B for spatial-reflection symmetric quantum circuits [25].Importantly, although it was derived for random QCA, it holds more generally, and we can show that, for sufficiently large v, it is valid for arbitrary entanglement Rényi entropy and any QCA, even without spacetime disorder.To see this, we consider the case of v > 0. Suppose that vt is larger than the range of the QCA at time t, which has index ind t due to its additivity, we can take a specific bilayer form [10] such that a local operator (LO) with respect to the entanglement cut can transform the operator into an array of virtual identities: Obviously, the operator entanglement, irrespective of its Rényi order, is the logarithm of the dimension of the virtual Hilbert space across the cut, which reads Similarly, for v < 0 we have α (vt, t) = ln(s −vt eq e ind t ) = −s eq vt + ind t , (S54) as anticipated.
As discussed in the main text, given the line tension E 2 (v) we can define the entanglement production rate as From a straightforward calculation, we obtain Argument 3 and r p is given in (S20).The main qualitative features of Γ 2 (s) can be see from Fig. S4, where we have plotted it for given values of d, p and q.Note in particular, that it satisfies the following relations As anticipated, our model makes it possible to also compute another useful measure of scrambling, the so-called tripartite mutual information [70], see also [64,[85][86][87].We recall that, given two bipartitions A ∪ B and C ∪ D with A = (−∞, 0), B = (0, ∞), C = (−∞, x), D = (x, ∞) for the input and output, respectively, the tripartite information is given by I A:C:D (x, t) ≡ I A:C (x, t) + I A:D (x, t) − I A:CD (x, t).Using the unitarity of U(t), we can simplify I A:C (x, t) into S max − S AD (x, t) − S AC (x, t), where S max = 2L ln d.Taking the Rényi-2 version, S AC (x, t) is nothing but S (o) 2 (x, t), whose exponentiated average has been obtained previously.Therefore, all we have to do is to compute S max − S AD (x, t).
While both S max and S AD (x, t) diverge in the thermodynamic limit, their difference remains finite.Moreover, we can again make use of the map to the Ising model, obtaining After some straightforward calculations, we can relate the above quantity to g L,R in Eqs.(S9) and (S10): E is given in Eq. (S44).Noting that | ind | = ln max{p/d, q/d}, we have We conjecture that Eq. (S62) is valid universally for generic anomalous dynamics with translation invariance on the ensemble level.A heuristic argument is as follows: for a while let us take the periodic boundary condition so that S max − S AD can be decomposed into two terms arising from the left and right boundaries (cf.Fig. S5): Suppose that the segments across the left and right boundary are taken to be of the same length and using the translation invariance of the random QCA on the ensemble level, we know from the entropy formula of the index [16] that As long as we can justify we immediately obtain S max − S AD 2| ind |t, implying t −1 (S max − S AD ) | ind | for the half-infinite entanglement cut.For model A, (S65) can be established for the Rényi-2 entropy and so it should hold for arbitrary Rényi (including von Neumann) entropies, under the assumption of small fluctuations.Therefore, since it holds for Haar-random QCA, we expect it to hold also for sufficiently chaotic dynamics.As a consistency check, we will see that Eq. (S62) holds true also for Model B.
A r e r C f r x X q 3 P q a t B W s 2 s w / + l P X 5 A 3 O D m z c = < / l a t e x i t > ind s eq < l a t e x i t s h a 1 _ b a s e 6 4 = " e b O N D B q H 0 M m 0 J y 9 t c N 3 a 5 E 0 q E 0 4 7 k k I p e l i p / y u Q h y + y R 7 5 T c i j s t u A h e D i W Q V 8 1 3 v j q 9 C C e c C I 0 Z U q r t u b H u p k h q i h n J 7 E 6 i S I z w E P V J 2 6 B A n K h u O l 0 o g 2 W j 9 G A Y S X O E h l P 1 9 0 S K u F J j H p h O j v R A z X s T 8 T + v n e j w o p t S E S e a C D x 7 K E w Y 1 B G c p A N 7 V B K s 2 d g A w p K a v 0 I 8 Q C Y W b T K 0 T Q j e / M q L 0 D i p e G c V 9 / a 0 V L 3 M 4 y i C A 3 A I j o E H z k E V X I M a q A M M H s E z e A V v 1 p P 1 Y r 1 b H 7 P W g p X P 7 I M / Z X 3 + A P 6 B m / o = < / l a t e x i t > v < l a t e x i t s h a 1 _ b a s e 6 4 = " K / H j 4 4 G + 5 B l T x Z 0 Y M r C i n M I 5 n 6 8 = " > A < l a t e x i t s h a 1 _ b a s e 6 4 = " x y T P S 9 s y 5 d / 9 q G y 4 P w m D p 1 c 5 k 0 l q q C T r i 6 K U Y x P j M k 8 8 Z 4 o S w 1 c W g C h m 3 4 r J F d g g j U 3 d t S H 4 2 1 / e h e F x x z / t e J 9 O G t 2 z T R w 1 9 A a 9 R S 3 k o 3 e o i y 5 Q H w 0 Q Q d f o B / q F 7 p w b 5 6 d z 7 / x e W / e c z c x r 9 F 8 5 f / 4 C j H a t L g = = < / l a t e x i t > v E < l a t e x i t s h a 1 _ b a s e 6 4 = " For completeness, let us also briefly discuss how the previous results read for Model B. Since the calculations are analogous, and fact significantly simpler, we omit the derivation.First, for the Rényi-2 operator entanglement we find e −S (o) 2 (x,t) = p t−x q t+x g(2t, t + x − 1, 1 − p R ) + where p L , p R are defined in (S34), while g(n, m, p) is given in (S36).This allows us to obtain the line tension E = E 2 (0) = log pq (p 2 q 2 − 1) 2 4pq(p 2 − 1)(q 2 − 1) , (S68) where the minimum is attained at v m = (p−q)(pq+1) (p+q)(pq−1) .Once again, it is straightforward to check that v E coincides with the growth rate of the Rényi-2 entropy after a quench from an initial product state, and that the identities (S49), (S50) hold.Finally, from the previous results, we also obtain an explicit expression for the entanglement production rate (S56) Γ 2 (s) = 2 log pq pq + 1 p + q − 2 log pq p(q 2 − 1)e − 1 2 (1+ ind seq )s + q(p 2 − 1)e where we defined gL = g(2t, t + x − 1, p L ) , gR = g(2t, t − x, p R ) , (S72) with p L , p R and g(n, m, p) given in (S34) and (S36) respectively.We see that when the spacetime coordinate is outside the (strict) light cone, i.e., x > t (or x < −t), we have e I A:C:D 2 (x,t) 1, consistent with the intuition that the quantum information cannot propagate faster than "light".Furthermore, for the translation QCA with q = 1 (p = 1), we have again e I A:C:D 2 (x,t) 1.This result is consistent with the observation that no (quantum) scrambling takes place in swap (permutation) circuits [70].Finally, if we take p = q, we reproduce the results of Ref. [64] for RUC, as we should.Taking now x = 0, we obtain E is given in Eq. (S68), implying the validity of Eq. (S62).

FIG. 1 .
FIG. 1. (a): Pictorial representation of LP evolution as the edge dynamics of a Floquet qudit system.In the presence of an MBL phase in the bulk[27], it is always possible to decompose the one-period Floquet operator as[4] UF = U edge e −iH bulk T , where H bulk is an MBL Hamiltonian, and U edge is an effective 1D evolution acting on qudits within a few localization lengths from the boundary.(b): Any LP evolution U(t) may be approximated by a QCA in Margolous form, i.e. the single time-step unitary operator U admits a bilayer representation, where the local unitaries map a product of nearestneighbor Hilbert spaces into another, with possibly different individual input or output dimensions.

FIG. 2 .
FIG. 2. (a) Pictorial representation of a spacetime curve cutting through the unitary evolution operator U(t), in (1 + 1)D.For a given curve, the total line-tension is seqE(v)dt, where v(t) is the local velocity.(b) The line tension may be obtained by viewing U(t) as a state, and computing the corresponding bipartite entanglement.
t e x i t s h a 1 _ b a s e 6 4 = " p H h c V e S u 1 u r U t e x i t s h a 1 _ b a s e 6 4 = " 5 a T r P J S Z u 2 b s W d a F l 2 i 4 +

Fig. 2 .
Formally, we introduce a doubled Hilbert space H = (C d ) ⊗2L ⊗ (C d ) ⊗2L along with the maximally entangled state |I = |φ + ⊗2L , where |φ + ≡ d n=1 |n ⊗ |n / √ d and n runs over a basis of C d .This allows us to vectorize the evolution operator as FIG. S1.Models of random QCA.a) Model A: The unitary evolution operator is given in Eq. (S1).For p = d the dynamics is anomalous, with index (S2).b) Model B: The dynamics is dictated by the single-step evolution operator Us = (⊗ L−1 j=0 u FIG. S2.Parent Floquet systems for (a) Model A and (b) Model B. (a) Compared to Ref.there are two main differences: (i) In the first step, the bipartite unitaries that map (by conjugation) physical degrees of freedom to virtual ones are randomized; (ii) there is an additional middle step of applying randomized two-site gates on the physical level.(b) After an on-site randomization, sublattices Λp (black circles) and Λq (white circles) undergo clockwise and counterclockwise four-step swap operations (following the order red→blue→green→orange; same in (a)), respectively.The operations on the two sublattices can be performed simultaneously.
FIG.S3.OTOC in Model A. Exact numerical results, obtained by evaluating (S8), are compared against the hydrodynamic description (S30) for different values of the index.Vertical lines correspond to the left/right light-cone velocities.The curves are in excellent agreement and we verified that the difference decreases uniformly in R as t increases.
t e x i t s h a 1 _ b a s e 6 4 = " 5 a T r P J S Z u 2 b s W d a F l 2 i 4 + G P 8 Y O o = " > A A A B 6 n i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K q M e g F 4 / x k Q c k S 5 i d 9 C Z D Z m e X m d l A W P I J X j w o 4 t U v 8 u b 0 g 3 o 3 S N p i 7 B n l 9 5 k b R P 6 v Z Z 3 b o 9 r T Y u Z 3 W U w Q E 4 B E f A B u e g A a 5 B E 7 Q A B o / g G b y C N + P J e D H e j Y 9 p t G T M Z v b B H x i f P 4 1 p n L o = < / l a t e x i t > s < l a t e x i t s h a 1 _ b a s e 6 4 = " n T l t u f U f Y m + b I p P O 4 i i 0 l C T 8 5 J s = " > A A A C I 3 i c b V D L S s N A F J 3 4 r P E V d e l m s B T q e e z 8 9 0 5 d H 6 7 m + 4 v 9 8 g 9 e Z S 6 z r T n E 3 o W b u 8 B a Z v I q A = = < / l a t e x i t > s eq < l a t e x i t s h a 1 _ b a s e 6 4 = " c G G b

4 2 1
/ e h e F x x z / t e J 9 O G t 2 z T R w 1 9 A a 9 R S 3 k o 3 e o i y 5 Q H w 0 Q Q d f o B / q F 7 p w b 5 6 d z 7 / x e W / e c z c x r 9 F 8 5 f / 4 C j H a t L g = = < / l a t e x i t > v E FIG. S4.(a) Entanglement line tension E(v) and (b) production Γ(s) (in terms of Rényi-2 entropy) for Model A with d = 6 and p 3, q = 12.The red dashed lines refer to the tangents at (vL,R, E(∓vL,R) = vL,R ± ind /seq); (b) (±seq, 0).The blue dashed line in (b) corresponds to Γ(s) − (ind /seq)(s/seq), whose maximum gives v (o) E .

2 (
S60) Since g L,R (almost) saturates 1 inside the light cone, the dynamics of the Rényi-2 tripartite information should be governed by e I A:C:D

8 <
l a t e x i t s h a 1 _ b a s e 6 4 = " 5 a T r P J S Z u 2 b s W d a F l 2 i 4+ G P 8 Y O o = " > A A A B 6 n i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K q M e g F 4 / x k Q c k S 5 i d 9 C Z D Z m e X m d l A W P I J X j w o 4 t U v 8 u b f O E n 2 o N G C h q K q m + 6 u I B F c G 9 f 9 c g o r q 2 vr G 8 X N 0 t b 2 z u 5 e e f + g q e N U M W y w W M S q H d e a V 2 n c d R h C M 4 h l P w 4 B J q c A t 1 a A C D A T z B C 7 w 6 w n l 2 3 p z 3 R W v B y W c O 4 R e c j 2 8 8 T o 3 E < / l a t e x i t > v R < l a t e x i t s h a 1 _ b a s e 6 4 = " 6 3 H D T J c C c e n f w 7 0 g 3 o 3 S N p i 7 B n l 9 5 k b R P 6 v Z Z 3 b o 9 r T Y u Z 3 W U w Q E 4 B E f A B u e g A a 5 B E 7 Q A B o / g G b y C N + P J e D H e j Y 9 p t G T M Z v b B H x i f P 4 1 p n L o = < / l a t e x i t > s < l a t e x i t s h a 1 _ b a s e 6 4 = " n T l t u f U f Y m + b I p P O 4 i i 0 l C T 8 5 J s = " > A A A C I 3 i c b V D L S s N A F J 3 4 r P E V d e l m s B T q

a 5 z 6 B
FIG. S6.(a) Entanglement line tension E(v) and (b) production Γ(s) (in terms of Rényi-2 entropy) for Model B with p = 3 q = 2.The red and blue dashed lines share the same meaning as those in Fig. S4.