Operator Entanglement in Local Quantum Circuits I: Chaotic Dual-Unitary Circuits

The entanglement in operator space is a well established measure for the complexity of the quantum many-body dynamics. In particular, that of local operators has recently been proposed as dynamical chaos indicator, i.e. as a quantity able to discriminate between quantum systems with integrable and chaotic dynamics. For chaotic systems the local-operator entanglement is expected to grow linearly in time, while it is expected to grow at most logarithmically in the integrable case. Here we study local-operator entanglement in dual-unitary quantum circuits, a class of"statistically solvable"quantum circuits that we recently introduced. We identify a class of"completely chaotic"dual-unitary circuits where the local-operator entanglement grows linearly and we provide a conjecture for its asymptotic behaviour which is in excellent agreement with the numerical results. Interestingly, our conjecture also predicts a"phase transition"in the slope of the local-operator entanglement when varying the parameters of the circuits.

for some manifestly non-local objects as well, such as, e.g., the time-dependent many-body propagator [38][39][40][41][42]. For this reason, in this work we will always refer to this quantity as local-operator entanglement.
The main problem is that, to date, there are essentially no exact benchmarks for the dynamics of the local-operator entanglement except for certain results in integrable models [30][31][32][33] and in random models in a particular asymptotic limit [37]. Providing such exact benchmarks for non-integrable many-body quantum dynamical systems is our main objective. We consider systems represented as local quantum circuits [43][44][45][46][47][48][49][50][51][52][53][54][55], i.e. qudit chains (quantum spin-(d − 1)/2 chains with arbitrary integer d ≥ 2) where the time evolution is generated by the discrete application of unitary operators coupling neighbouring sites. We measure the local-operator entanglement through Rényi entropies at integer Rényi order. Our strategy is to write them in terms of partition functions on a non-trivial space-time domain that we contract in terms of row and corner transfer matrices. We divide this endeavour into two separate works investigating two conceptually distinct classes of quantum circuits that are generically non-integrable. Using the special properties of these classes we prove and conjecture some exact statements about dynamics of local-operator entanglement.
In the present paper, we study the dynamics of local-operator entanglement for dualunitary local quantum circuits. This is a class of local quantum circuits where the dynamics remains unitary also when the roles of space and time are swapped [56]. As we showed in a recent series of works, dual-unitarity is an extremely powerful property and enables the exact calculation of many statistical and dynamical properties. These include spectral statistics [57], (state) entanglement spreading [58], and dynamical correlations [56] (see also [54] and [59] for other useful features of dual-unitarity). Focussing on dual-unitary quantum circuits with no local conservation laws -the chaotic subclass -we conjecture a general formula for the dynamics of the local-operator entanglement. The idea is to compute the local-operator entanglement by considering separately the entanglement produced by the two edges of the spreading operator -as if the opposite edge were effectively sent to infinity. Dual-unitarity allows us to evaluate these contributions exactly revealing a simple and remarkable prediction, which is in excellent agreement with exact short-time numerical results. First, we find that in chaotic dual-unitary circuits the local-operator entanglement always grows linearly with time. Second, the slope of growth displays an abrupt transition when varying the parameters of the circuits. Third, the slope is maximal on one side of the transition. This has to be contrasted with the linear local-operator entanglement growth in Haar-random noisy circuits [37], where the slope is around half of the maximal one. These results once again put forward dual-unitary circuits (in appropriate parameter ranges) as minimal models -with fixed local Hilbert space dimension and local interactions -for the maximally-chaotic dynamics.
In the companion paper [60] (Paper II) we consider the dynamics of local-operator entanglement in local quantum circuits exhibiting local dynamical conservation laws, i.e. solitons. These conservation laws are generically not enough to generate an integrable structureà la Yang-Baxter. Limiting to the circuits of qubits (d = 2), we classify all instances of circuits with solitons and show that if a spreading operator crosses some soliton, the dynamics of its local-operator entanglement can be computed explicitly and exhibits saturation. Interestingly, we show that all circuits admitting moving solitons are dual-unitary. Importantly, since they have conservation laws, those dual-unitary circuits are not chaotic as the ones studied here.
The rest of this paper is laid out as follows. In Section 2 we give a detailed definition of the quantum many-body systems of interest for this work -local quantum circuits -and introduce a useful diagrammatic representation to study their dynamics. In Section 3 we introduce the local operator entanglement entropies and write them in terms of partition functions on appropriate space-time surfaces. In Section 4 we specialise the treatment to dual-unitary local quantum circuits, recalling their main defining features and characterising the "completely chaotic" class of interest in this paper. In Section 5 we formulate our conjecture and use it to explicitly compute the local-operator entanglement dynamics (explicitly comparing it with the numerical evaluation of the space-time partition functions). Finally, Section 6 contains our conclusions. Five appendices complement the main text with a number of minor technical points.

Local Quantum Circuits
In this work we consider periodically-driven quantum many-body systems represented as local quantum circuits. These systems consist of a periodic chain of 2L sites, where at each site is embedded a d-dimensional local Hilbert space H 1 = C d so that the total Hilbert space is (1) The time evolution in the system is discrete and each time-step is divided into two halves. In the first half the time evolving operator is where Z L = Z ∩ (−L/2, L/2], while U x,y ∈ U (H 1 ⊗ H 1 ) is the unitary "local gate" connecting the qdits at sites x and y and encoding all physical properties of a given quantum circuit. In the second half, instead, the system is evolved by where T is a −periodic translation by one site Here {o j } are generic operators in H 1 . In summary, the "Floquet operator" -the time evolution operator for one period of the drive (one time-step) -is given by Note that, since the local gate is unitary, the Floquet operator U is also unitary. Moreover, from the definition (5) it immediately follows that U is invariant under two-site shifts Note that in this work we consider translationally invariant quantum circuits which are specified by a single 2-qudit gate U x,x+1/2 = U for all x ∈ 1 2 Z 2L , while we expect that the formalism we develop here should be useful also for generalizations to disordered and/or noisy quantum circuits.
Local quantum circuits admit a convenient diagrammatic representation. One depicts states as boxes with 2L outgoing legs (or wires) representing the local sites and operators as boxes with a number of incoming and outgoing legs corresponding to the number of local sites they act on. Each leg carries a Hilbert space H 1 . For instance, the identity operator on a single site, 1, is represented as while a generic single-site operator a is represented as The local gate and its Hermitian conjugate are instead represented as where we added a mark to stress that U and U † are generically not symmetric under space reflection (left to right flip) and time reversal (up to down flip, transposition of U ). The time direction runs from bottom to top, hence lower legs correspond to incoming indices (matrix row) and upper legs to outgoing indices (matrix column). With these conventions, the diagrammatic representation of U reads as where we labelled sites x by half integers, x ∈ 1 2 Z 2L , and boundary conditions in space (horizontal direction) are periodic. This means that the ultralocal operator evolved up to time t is represented as Before concluding we note that time-evolving operators transform covariantly under the following gauge transformation in the space of local gates Specifically, we have

Operator-to-state mapping
The time evolution of operators in H can be mapped into that of states in the "doubled" Hilbert space H ⊗ H by performing an operator-to-state (or vectorization) mapping Choosing any basis {|n } of H we completely specify the mapping by defining so that the time evolution maps to a y (t) −→ |a y (t) ≡ n,m n|a y (t)|m |n ⊗ |m * = (U † ⊗ U † * ) t |a y .
The complex conjugation (·) * is defined such that * n|O * |m * = ( n|O|m ) * , (18) meaning that the vectorization mapping is linear (and not antilinear!) with respect to both, the ket and the bra parts 1 . For convenience we arrange the states |n ⊗ |m * in H ⊗ H in such a way that the time evolution generated by U † ⊗ U † * is "local in space". Specifically, where {|i ; i = 1, 2, . . . , d} is a real, orthonormal basis of H 1 . In general, for any set of states |a , |b · · · ∈ H 1 , we use a compact notation |a b . . . = |a ⊗ |b ⊗ · · · . The mapping defined in this way is directly represented by folding the circuit where each thick wire carries a d 2 dimensional local Hilbert space and we introduced the "double gate" Note that the red gate is upside down, meaning that U is transposed (c.f. (U † ) * = U T on the r.h.s. of (17)). Finally, we also introduced the (normalised) local states associated to to the identity operator and to the operator a a −→ a = a ≡ |a .
We stress that in this paper we always consider local operators that are Hilbert-Schmidt normalised tr[aa † ] = 1 .
For non-normalised operators one should include the appropriate normalisation factors in (20) and (24). Finally, we remark that from the unitarity of U it follows where we introduced

Local Operator Entanglement
The entanglement of a time-evolving operator O(t) is defined as the entanglement of the state |O(t) corresponding to it under the state-to-operator mapping. Specifically, here we are interested in the entanglement of a connected real space region A with respect to the rest of the system. Since the state corresponding to a time-evolving operator is pure, this quantity is conveniently measured by the Rényi entanglement entropies [61] S (n) where ρ A (t) is the density matrix at time t reduced to the region A. Specifically, here we consider the evolution of the entanglement of the ultralocal operator a y and select half of the chain A = [0, L/2). Moreover, here and in the following we will always take a to be Hilbert-Schmidt orthogonal to the identity operator, i.e. traceless, to project out its trivial component.
With our choices of operator and subsystem the graphical representation for the reduced density matrix reads as In the representation (29) we took y < t ≤ L. We considered the right inequality, because we are interested in the thermodynamic limit and the results no longer depend on L for L ≥ t, while we take the left inequality, t > y, because in the opposite case the reduced density matrix is pure and hence the entanglement vanishes. This is due to the fact that in quantum circuits there is a strict lightcone for the propagation of information: nothing can propagate faster than a given maximal velocity (this is stricter than the Lieb-Robinson bound which allows for exponentially small corrections). In particular, in our units (see Eq. (20)) the maximal velocity is 1. Finally, we assumed y to be an integer. The case of half-integer y can be recovered by the reflection R of the chain around the bond between 0 and 1/2. This results in where S is the "swap-gate" and we designate explicitly the dependence on the local gate. From now on we always take y to be integer. Using the representation (29) we see that the calculation of tr A [ρ n A (t, y; a)] is reduced to that of a partition function of a vertex model (generically with complex weights). For instance, in the simplest nontrivial case n = 2 we have Using the graphical rules (26) this equation is reduced to This expression can be rewritten as where we introduced the "light-cone coordinates" and the row/column transfer matrices Computing higher moments (i.e. higher Rényi orders) requires n > 2 replicas and involves partition functions on more complicated surfaces. In order to represent them compactly it is convenient to introduce the d 2x − × d 2x + "corner transfer matrix" [62,63], defined by the following matrix elements where {|I ; I = 1, 2, . . . , d 2 } is an orthonormal basis of H 1 ⊗H 1 . Note that the corner transfer matrix is related to the row transfer matrices H In terms of C[a] we can easily express the Rényi entropies as follows The problem of computing operator entanglement is then reduced to that of computing the Before concluding this section we note that under the gauge transformation (13) the traces of the reduced transfer matrix transform as follows This means that the gauge transformation only causes a rotation in the space of ultralocal operators.

Completely Chaotic Dual-Unitary Circuits
In this paper we consider dual-unitary circuits, i.e. local quantum circuits where the evolution remains unitary upon switching space and time directions. This means that the local twoqudit gate U remains unitary if we consider left pair of wires as incoming states and right pair of wires as outgoing states. More formally, defining the "dual" (space) propagatorŨ by means of the relation the circuit is dual-unitary, if both, U andŨ are unitary [56]. Dual unitarity can be expressed explicitly as or diagrammatically as = , = . (46a) Considering the double gate (27), these relations lead to We have shown in [56] that the dual-unitarity condition is not as stringent as one might think. For instance, in the case of qubits (d = 2) it only fixes two parameters of the fifteen specifying a generic matrix in U ∈ U(4) and allows for a rich variety of dynamical behaviours [56]. Here, in particular, we focus on a specific class of dual-unitary circuits, which we term the completely chaotic class. To define it, we consider the transfer matrices V x [1] and H x [1]. Since any such transfer matrix is a contracting operator, i.e.
T |v ≤ |v , their eigenvalues are contained in the unit circle of the complex plane (see Appendix A for a proof of (48)). Using only the relations (47a) and (47b), we can find x + 1 independent simultaneous eigenvectors of V x [1] and H x [1] associated with eigenvalue one. They read as where we introduced the "rainbow" states |r l and their orthonormal counterparts |r l , satisfying r k |r l = δ k,l . Note that the hermitian conjugates of these vectors are always left (49) are right eigenvectors if the circuit is dual-unitary.
We are now in a position to introduce the following We stress that (49) are in general only a subset of the eigenvectors of V x [1] and H x [1] associated with eigenvalue one. For instance, integrable dual-unitary circuits (e.g. the oneparameter dual-unitary line of the integrable trotterised XXZ model [64], or the self-dual Kicked Ising model at the non-interacting point) have much more such eigenvectors (see Paper II for additional examples of such circuits). A thorough numerical analysis, however, proves that (i) the completely chaotic class exists; (ii) it is the generic case. In other words, generating a dual-unitary gate at random we will find with probability 1 that there are no eigenvectors of V x [1] and H x [1] with unit magnitude eigenvalues other than (49). The rest of the spectrum is gapped within a circle of radius strictly smaller than one.
Before moving on to the calculation of the local operator entanglement dynamics, it is interesting to investigate the relation between the definition 4.1 of completely chaotic circuits and the intuitive definition of chaos based on absence of local conservation laws. We will show that the class of completely chaotic dual-unitary circuits is in general more restrictive than that of chaotic ones. Namely, if a dual-unitary circuit has some non-trivial local conservation law V x [1] and H x [1] acquire some additional eigenvectors corresponding to the eigenvalue 1. In our discussion we will focus on circuits admitting conservation laws with local density which can be written either as or as where the local densities q ± x act non-trivially (have support) on r sites. More precisely, these densities act non-trivially on the intervals [x, x + (r − 1)/2] ∩ Z L /2 and [x − (r − 1)/2, x] ∩ Z L /2 respectively. Moreover, we choose the densities such that tr x [q ± x ] = 0 (here the trace is over the local Hilbert space at the x-th site). Note that this can be done without loss of generality: all charges can be written as combinations of Q + and Q − .
Due to the two-site shift symmetry of the time evolution in the system we considered local conservation laws obtained by summing only on a sub-lattice (say even sites). In order for Q ± to be conserved their local densities must satisfy continuity equations of the form for some "currents" J ± x supported on r + 1 sites (for concreteness in writing (54) and (55) we assumed r odd). As shown in Appendix B in dual-unitary circuits the relations (54)  . This means that conserved-charge densities in dual-unitary circuits satisfy either or Let us show that these relations imply that V x [1] and H x [1] have additional eigenvectors corresponding to the eigenvalue 1. Focussing on the first and writing it in diagrammatic form for r = 5 we have Tracing out the identities in the last three sites and repeatedly multiplying by the double gate W we find Finally, contracting the last two sites with |•• we find We then have that the vector is an eigenvector of (V Note that |v cannot be zero. Indeed, if this were to be the case also the l.h.s. of (59) would vanish leading to an absurd: the r.h.s. of that equation features the non-zero operator q − x conjugated by unitary matrices. To conclude or argument we note that is an eigenvector of V •• x corresponding to the eigenvalue 1. Then, we can construct many additional eigenvectors of V x≥r [1] corresponding to eigenvalue 1, for example Finally we note that, if (56) holds, the charges are also conserved. Considering the densities of such charges and proceeding as before we can construct exponentially many (in x) eigenvectors of V x≥r [1] with eigenvalue 1. An analogous reasoning considering a conserved density q + x would instead produce additional eigenvectors of H x≥r [1] corresponding to eigenvalue 1.

Dynamics of Local Operator Entanglement
It is generically very difficult to say much about the dynamics of local operator entanglement in interacting systems (in fact, this quantity is generically out of reach even in the presence of integrability). Here we show that in completely chaotic dual-unitary circuits one can make some quantitative progress.
In the first part, we prove that in the two limits x ± → ∞ the local operator entanglement can be determined exactly. These two limits correspond to varying the initial position of the operator in order to measure the entanglement generated at the edges of the light-cone (x − → ∞ gives the entanglement generated by the right edge and x + → ∞ that generated by the left: see Fig. 1 for a pictorial representation). Note that, the operator "breaks" the left-right symmetry of the problem and one should not expect the results of the two limits to coincide. Indeed, we find that they are physically very different. In particular, while the entanglement generated by the right edge has flat spectrum and grows at the maximal speed, the one generated by the left edge is much richer. First it has a non-trivial spectrum and second, while the von Neumann entropy always grows at the maximal speed, higher Rényi entropies show a phase transition in the speed of the entanglement growth when varying the parameters of the gate. Specifically the growth depends on the largest eigenvalue λ governing the decay of the dynamical correlations (cf. [56]).
In the second part of the section we show that the "local operator n-purities" for any x + and x − are well described (even at short times) by summing the two limits x ± → ∞, namely This indicates that in completely chaotic dual unitary circuits the bulk of the light-cone region rapidly becomes highly entangled (and hence it does not contribute to the purities) and the leading contributions to the purities arises from the edges. Interestingly, if the dynamical correlations of a circuit decay fast enough we observe the local operator entanglement growing at the maximal speed (log d 2 ); otherwise the growth is slower and depends on λ. In the latter case the entanglement spectrum is non-trivial. Finally we extend the above results to a class of chaotic but not completely chaotic dual unitary circuits including the self dual Kicked Ising model. In particular, we compute exactly the operator entanglement in the two limits x ± → ∞ and, by comparing with numerical simulations, we show that the property (67) continues to hold far enough from integrable points.

The Two Limits
Let us start by considering the special limits described above, where one focusses on the entanglement generated by the edge of the light-cone produced by the spreading operator a.

The Limit x − → ∞
Let us first consider the entanglement generated by the right light-cone edge, namely we consider the limit x − = t − y → ∞ while keeping x + = t + y fixed. In this limit it is convenient to use the representation (40). We start by nothing that, since the operator a is traceless (cf. Sec. 3), the only eigenvector of H x [1] with non-zero overlap with the "initial state" |a † , •, ..., •, a is |e x + = |r x + . This is because the scalar product e x + |a † • · · · •a is the only one among { e i |a † • · · · •a } i=1,...,x + which does not produce the trace of a. In particular Figure 1: Pictorial illustration of the operatorial evolution, depicting the two limits (73) and (76). The two limits are taken along the dashed arrows. The first limit x − → ∞ has x + constant, which means that as time increases we move the operator to the left getting contribution from the red region. The second limit x + → ∞ has x − constant, meaning that, as time increase we move the operator to the right and get contribution from the blue region.
Where we used that a is Hilbert-Schimdt normalised. Plugging in the definition of C[a] † C[a] we then have where we introduced the d 2x + × d 2x + matrix M k (k ∈ {0, 1, . . . , x + }) defined as M k can be expressed in matrix form as where 1 denotes the identity element in End(H 1 ⊗ H 1 ) and is a projector to the state • on site k. Plugging (69) into (42) we find that the end result for a Rényi entropy of generic order n reads lim x − →∞ S (n) (y, t) = lim where, to take the trace, we used Eq. (73) gives linear growth of the operator entanglement entropy with the maximal slope, and holds in the absence of "non-generic" eigenvectors of eigenvalue 1.

The Limit x + → ∞
Let us now consider the limit x + = t + y → ∞ while keeping x − = t − y fixed. This limit can be evaluated using (41) but we immediately see that it is more complicated than the previous one: we need to deal with the operator-dependent transfer matrix V x − [a] and all of x − + 1 eigenvectors. The calculation yields Therefore we obtain Once again, this result holds in the absence of additional eigenvectors of V x − [1] with unit magnitude, i.e. for completely chaotic dual-unitary circuits. The missing information in (76) is the value of e k |V x − [a]| • · · · • , which can be expressed in terms of This expression can be evaluated by writing the elements of the sum using the single qudit map introduced in [56] for calculating the dynamical correlation functions. The central part of (77), which results from the contraction with the rainbow state |r l , simplifies due to the unitarity of the gate, and produces a factor d −l . The rest of the expression can be written in terms of the maps as The maps can be expressed using a d 2 ×d 2 matrix and the expressions are then easily evaluated numerically. Moreover, the maps M −,U † and M +,U have the same eigenvalues and their respective eigenvectors |e −,U † and |e +,U are connected via the relation |e +, The leading asymptotic behaviour is governed by the leading eigenvalue 2 λ (|λ| ≤ 1) of the map M +,U and can be determined analytically by posing in (79) and (80). Here c l is bounded in l, i.e.
lim sup Plugging in (76) we find the following asymptotic result ∆S (n) | asy,x + = lim The result is intriguing, we see a transition between maximal and a sub-maximal growth, governed by the slowest decay of the two point dynamical correlation functions. Moreover, we see that the entanglement spectrum is not flat in this limit, but the result encodes a non-trivial n-dependence, see Fig. 2. This is very different from the limit x − → ∞, where all entropies experience maximal growth. Furthermore, there is another interesting observation to make.
Performing an analytical continuation of the result in n and taking the limit n → 1 + we find that the the growth of von-Neumann entropy (n → 1 + ) is always maximal.

The Conjecture
Let us now consider the local operator entanglement for generic x − and x + . To describe its leading in time behaviour we propose the following conjecture Conjecture 5.1. For chaotic dual-unitary local circuits, at long times the operator entanglement entropies for n > 1 are well described by the sum of the two limits (73) and (76). Namely As pictorially represented in Fig. 3, the conjecture consists in replacing the trace of the n-th power of the reduced density matrix of the "vectorised" spreading operator a y (t) with the sum of two terms. These terms are the trace of n-th powers of density matrices corresponding to operators obtained from a y (t) by sending to infinity respectively its left (−x − ) or right Figure 2: The asymptotic slope ∆S (n) | asy,x + (83) as a function of the gate parameter J (see Appendix C for details on the parametrisation) for different values of n (different colors). The slope is n-independent in the maximal-growth region but both the size of this region and the slope away from it depend on n. (x + ) edges. Note that the conjecture cannot hold for the von-Neumann entropy, as the limit n → 1 + of the r.h.s. of (84) is singular (the argument of the logarithm goes to 2). Conjecture 5.1 yields the following form for the entanglement entropies where the "slope" ∆S (n) (y, t) and the "offset" µ n (y, t) are bounded in t. We evaluated Eq. (85) using Eq. (73) and Eq. (76) and compared it to the results of exact short-time numerical simulations -obtained by direct diagonalisation of the corner transfer matrix (see Appendix E for details). The comparison, for n = 2 and y = 0, is reported in Fig 4. The figure presents results for both the slope ∆S (2) (0, t) and the constant shift µ 2 (0, t), which is very sensitive to small errors in the slope. The agreement observed is remarkable, even for the short times accessible by the numerics. A similar level of agreement is observed also for n ≥ 2.
The asymptotic value of the slopes in the limit t → ∞ with fixed "ray" ζ = y/t are given by lim t→∞ y/t=ζ Equation (86) predicts a "phase transition" between a region with slopes that are symmetric in ζ (as it happens for random unitary circuits, see Ref. [37]) and a region where instead they show an interesting asymmetry in ζ, and, moreover, they become n-dependent. In particular, we see that for ζ = 0 the slopes coincide with those given in (83) and the slope in the symmetric region is the maximal one (log d 2 is the maximal entanglement growth attainable in a circuit with d-dimensional local Hilbert space). For comparison, we computed numerically the dynamics of the local operator Rényi-2 entropy in Haar-random non-dual-unitary local qubit circuits for ζ = 0. We considered two cases: (i) we chose the same (constant) gate for all the space-time points (clean case); (ii) we took different i.i.d. U x,t for each space-time point in the circuit (noisy case). In both cases we obtained roughly half of the maximal slope of the entanglement entropy growth, see Fig. 5. This is in accordance with the predictions of Ref. [37] and proves that, in some parameter ranges, dual-unitary circuits are "more chaotic" than the average. The idea behind Conjecture 5.1 is most easily explained considering the purity e −S (2) (y,t) . Looking at the representation (33) we see that this quantity can be written as the partition function of a statistical mechanical model (with complex weights) on a rectangle of dimensions 2x + and 2x − . The conjecture corresponds to restricting this partition function to the sum over configurations spanning eigenvectors of eigenvalue 1 of both row and column transfer matrices. The same idea applies to n-purities with n > 2. Physically, this corresponds to assume that the bulk of the light-cone is highly scrambled (i.e. it gives a very small contribution to the purity), while the regions at a finite distance from the light-cone edges present the minimal scrambling (i.e. give the leading contribution to the purity). This is justified by noting that close to the light-cone edge the operator retains the maximal amount of information on the initial condition. We expect this picture to hold true for more general, non-dual unitary, chaotic systems if one replaces the light cone spreading at the maximal speed (1 in our units) with an effective one spreading at the "butterfly velocity" v B [6] of the system. Indeed, v B is by definition the velocity at which the scrambled region spreads in time.

Self-Dual Kicked Ising Model (d = 2)
Conjecture 5.1 is assumed to describe the asymptotic dynamics of the local operator entanglement in any chaotic dual-unitary circuit. In order for it to have any predictive power, however, one must be able to compute the limits x ± → ∞. While in the previous subsections we showed that this can be done for the completely chaotic subclass, here we show that the limits can be computed exactly also in the paradigmatic example of dual-unitary circuits in d = 2: the self-dual kicked Ising model [57,58]. This model is not completely chaotic according to Definition 4.1 because it possesses additional structure. Specifically, its local gate fulfils where w, α are some hermitian and traceless matrices in SU(2) 3 . This condition leads to x additional eigenvectors of the transfer matrices V x [1] and H x [1] with eigenvalue one. In fact, Figure 4: The slope ∆S (2) (y, t − 1/2) = S (2) (y, t) − S (2) (y, t − 1) and the constant offset factor µ 2 (y, t) (much more sensitive) versus the parameter J for a dual unitary gate with r = 0.5, φ = 0.7, θ = 0 (see Appendix C for the definition of the gate). We show the results for operators a 1 = σ 3 (left panel) and a 2 = α 1 σ 1 + α 2 σ 2 + α 3 σ 3 a fixed random operator with α 1 = 0.3289, α 2 = 0.0696, α 3 = 0.6221. The points correspond to exact numerical results, and the lines are the predictions using the conjecture (84). The operator is initialised at y = 0, and we set t = 7 for the right panel. Figure 5: The slope of Rényi n = 2 entanglement entropy for the operator σ 3 evolving according to (non-dual-unitary) U (4) Haar-random gates. In the clean case we average over 10 realisations, in the noisy case over 20 (100 for t ≤ 6). The results suggest that the slope is close to log 2, which is half of the maximal slope. Note that this agrees well with large d result from [37], where we get the slope 2 6 5 ( √ 2 − 1) log 2 ≈ 0.9941 log 2, if we use the parameters s spread , v b for d = 2 (cf. Ref. [37] for a definition of these parameters).
as shown in Appendix D, all reflection symmetric dual-unitary circuits fulfilling (87) are gauge equivalent to the self-dual kicked Ising.
We can use the gauge transformation (13) to set α = σ 3 , which holds in the standard formulation of self-dual kicked Ising model. The additional eigenvectors with eigenvalue one are then given by where 3 stands for the operator σ 3 . To construct an orthonormal basis, we consider the following linear combination Having the eigenvectors, we evaluate the limits: where we parametrised the initial local operator as The last equality in Eq. (91) follows from e j | V x + [a] |• · · · • = 0. Therefore, the additional eigenvectors change only the constant prefactors. Equations (90) and (91) show that, in the long time limit, the offset constant µ 2 is different with respect to the result in completely chaotic dual-unitary circuits (cf. Eqs. (73) and (76)), but the slope is the same. With the limits x ± → ∞ at hand we are now in a position to compare the prediction of Conjecture 5.1 with (short-time) numerics. A comparison is shown in Fig. 6. The figure shows that, far enough from some special points in parameter space (see the caption), there is good agreement even for numerically accessible times (t ≤ 8).

Conclusions
In this paper we studied the local operator entanglement growth in dual-unitary circuits. We identified a completely chaotic class for which the local operator entanglement always grows linearly in time. For this class we provided a quantitative description of the localoperator-entanglement dynamics based on a simple conjecture, which is strongly supported by numerical results. We postulated that, at late enough times, the local operator purities (traces of powers of the reduced density matrix) can be determined by considering separately the entanglement produced by the two edges of the spreading operator and then summing them together. In other words, we wrote the exponentials of the operator entanglement entropies as sums of two contributions, respectively obtained by sending the right edge of the spreading operator to infinity and the left edge to minus infinity. Our conjecture, together with the dual-unitarity property, allows to evaluate analytically the local operator entanglement of generic operators initially localised on a single site. These results have been extended to Figure 6: Prediction and numerical data for self-dual kicked Ising model. The points correspond to exact results, and the lines are the predictions using the conjecture (84). In the left panel we show the Rényi n = 2 entropy at times 6, 7, 8 versus the magnetic field parameter h. The deviations from the prediction observed away from the central region are due to the vicinity of solvable points h = k π 4 , k ∈ Z, where the model is a Clifford circuit. There the prediction fails, because the transfer matrix has additional eigenvectors of eigenvalue 1. But as time increases, the region where the prediction holds grows. In the right panel we show the more sensitive constant offset factor µ 2 versus the parameter h at t = 7. The operators a 1 , a 2 are the same as those used in Fig. 4 and they are initialised in y = 0. the self-dual kicked Ising model (which does not fall into the completely chaotic class). We argued that a modified form of our conjecture should hold in generic chaotic systems, i.e. for non-dual-unitary circuits, however, without dual-unitarity it does not directly yield analytical predictions.
Interestingly, our conjecture predicts that the slope of the local operator entanglement displays an abrupt transition when varying the parameters of the circuits. Moreover, the point in which the transition occurs depends on the Rényi index. On one side of the transition the slope of growth is the maximal allowed by the geometry of the circuit (log d 2 ), which is approximately twice as large as that observed in Haar-random circuits [37]. This indicates that a subset of our chaotic dual-unitary circuits can be regarded as minimal solvable models for the maximally chaotic dynamics. On the contrary, on the other side of the transition the slope is not maximal, depends on the Rényi index, and approaches 0 when the dual unitary gate approaches the SWAP gate.
Our work raises a number of questions that can guide future research. First, our conjecture seems to describe the numerics even at small times, suggesting that it holds up to very small corrections. It would be interesting to investigate this aspect further and, possibly, rigorously prove the conjecture. Second, the class of systems that we introduced here (see also [56]) can be used to study exactly many aspects of non-equilibrium dynamics in chaotic systems, from relaxation of local observables to the behaviour of out-of-time-ordered correlations.
which generate a two-parameter r ± family of models (the phase ψ ± is irrelevant). The reflection symmetric case r + = r − = cos h is therefore gauge equivalent to the self-dual kicked Ising model, with h being the magnetic field in the z direction. This is also seen in the eigenvalues of the maps M ±,U (cf. (78)), which exactly match.

E Numerical methods
Calculating the operator entanglement entropy numerically is computationally expensive with resources scaling exponentially with t. In our case, we iteratively constructed the corner transfer matrix C[a], as defined in (39). First we construct the doubled gate W , from which we build the first row of C[a]. Then we add additional precomputed rows via matrix computations until we end up with the final corner transfer matrix. In the last and by far the most expensive step we calculate d 2(x + +x − ) matrix elements, with each costing d 2x + operations. At y = 0, d = 2 the total cost scales as 2 6t , which is still much better than using the row/column transfer matrices H x [a] and V x [a], where the cost scales as 2 8t .