Entanglement growth in diffusive systems

Entanglement helps in understanding diverse phenomena, going from quantifying complexity to classifying phases of matter. Here we study the influence of conservation laws on entanglement growth. Focusing on systems with U(1) symmetry, i.e., conservation of charge or magnetization, that exhibits diffusive dynamics, we theoretically predict the growth of entanglement, as quantified by the Rényi entropy, in lattice systems in any spatial dimension d and for any local Hilbert space dimension q (qudits). We find that the growth depends both on d and q, and is in generic case first linear in time, similarly as for systems without any conservation laws. Exception to this rule are chains of 2-level systems where the dependence is a square-root of time at all times. Predictions are numerically verified by simulations of diffusive Clifford circuits with upto ~ 105 qubits. Such efficiently simulable circuits should be a useful tool for other many-body problems. Entanglement is one of the most fascinating quantum properties, being a resource in quantum computation and an obstacle in classical simulations. The author analyses entanglement growth in terms of purity, and verifies the results by introducing a class of random Clifford circuits that conserve magnetization, finding that for systems with diffusive dynamics that entanglement growth differs depending on the dimension as well as between qubits or qutrits.

E ntanglement is one of crucial quantum resources responsible for the emerging 2nd quantum revolution-exploiting quantumness to perform tasks not possible by classical means, for instance, quantum computation, teleportation, or secure communication 1 . Even if not easily measurable 2 , it is an extremely powerful theoretical concept. This was further underlined by another discovery from the '80, from a seemingly unrelated field, namely the quantum Hall effect 3 . It gradually brought to light the fact that there can be phases of matter that have topological order which goes beyond the Landau's paradigm of classifying all phases of matter just by local order parameters. Today we understand that such topological order is connected to certain patterns of entanglement 4 . A modern view in fact uses entanglement to distinguishing different phases of matter 5,6 . Entanglement though plays a role also beyond the equilibrium phases. An example is for instance a putative non-thermal manybody-localized phase 7 , one of the distinguishing features of which is slow logarithmic-in-time growth of entanglement 8 .
Conservation laws and the associated symmetries are one of the most important properties of laws of physics. On the smallest scale, the elementary particles differ by their symmetries, and on the large scale, as well, the most violent objects we know-black holes-are believed to be defined only by their conserved quantities, charge, mass and angular momentum 9 . Furthermore, the symmetry to translations in time and its associated generator is the very object that governs dynamics. In short, symmetries are crucially responsible for the simplicity of nature at its core.
An important question is what role do conservation laws play on the dynamics of entanglement? Its growth with time is important also from a practical point of view. Namely, if it is small then efficient classical simulation of such systems is possible 10 . For generic local systems and initial states one expects that dynamics explores the whole available Hilbert space and therefore entanglement grows linearly with time. This holds true even for integrable systems, see e.g. refs. 11,12 . Because symmetries are about constraints, and because entanglement is given essentially by the number of degrees of freedom (two-level systems) involved, one might argue that symmetries will certainly affect entanglement growth. On the other hand, however, in the thermodynamic limit (TDL) one could also argue that conservation of a single charge should not matter much in a large Hilbert space.
Therefore it was surprising and interesting when it was shown 13 (focusing on diffusive 1D systems with conserved charge) that the entanglement, as quantified by the Rényi entropy S r ðtÞ :¼ grows in fact as S 2 $ ffiffi t p starting from a generic separable initial state, instead of the "expected" S 2~t , see also refs. 14,15 . This finding, if holding for generic systems, would have many consequences. For instance, one could argue that simple charge conservation causes the "Rényi complexity" $ 2 S 2 to grow only as $ b ffi ffi t p , i.e. slower than exponentially (though still super-polynomially). A system with diffusive conserved charge would seem to be a less powerful quantum information resource than a one without it.
We address the question of the Rényi entropy growth in local systems in any spatial dimension d and for any local Hilbert space dimension q. Theoretical predictions for a bipartition (L is the linear system size, being also proportional to the subsystem size), are numerically verified on large systems, with the total number of qubits up-to e.g. 252 × 252 ≈ 6 × 10 4 in 2D, and 48 × 48 × 48 ≈ 10 5 in 3D. While we confirm 13 , the main finding is that in higher d and q the asymptotic ffiffi t p growth is in fact not what one will typically observe. While for diffusive qubit systems (spin-1/2, i.e. q = 2) the asymptotic growth is still $ ffiffi t p , it starts in d > 1 only at a time when the entropy already becomes extensively large, S 2L d−1 , in other words, in the TDL the S 2 grows linearly with time at any finite value of the entropy. For qudit systems (q > 2), even in d = 1, one expects instead a linear asymptotic growth, except in cases where the dynamics of all diagonal operators is diffusive. In this respect the often studied qubit systems in 1D are rather special-diffusive growth is there observed already at early times and for all diffusive systems because a single diffusive charge, together with an identity operator, already exhausts the algebra of diagonal operators. As a side result, the presented new class of efficiently simulable systems with nontrivial dynamics could be useful in addressing other questions of many-body physics.

Results
Theoretical prediction. A class of systems that we study are lattice systems with local nearest-neighbor (n.-n.) interactions in d spatial dimensions and with q-dimensional local Hilbert space, whose dynamics has a nontrivial conservation of the total particle number or the total spin in z-direction (i.e., a U(1) symmetry). The dynamics of that conserved degree of freedom is assumed to be diffusive, while the rest of dynamics is generic (we exclude integrable systems). Specifically, the influence of possible non-U(1) symmetries is left for future work. In numerical demonstrations we also focus on Floquet systems in order to avoid having to deal with an additional conserved quantity (the energy). Linear dimension is denoted by L, and the total number of qudits by n: = L d .
We shall discuss the entropy growth as quantified by the Rényi entropy S r (integer index r > 1) starting from a pure product initial state. We prefer S r over von Neumann entropy S r→1 due to its analytical simplicity. In generic systems all S r , including S 1 , are expected to behave in the same way, whereas for diffusive systems the S 1 (which we don't discuss) can behave differently 13 . We will mostly focus on S 2 as a representative case of S r>1 . We remark that sometimes S 2 , rather than e.g. S 1 , is a more relevant quantity 16 , and is furthermore also easier to measure 17 . Using a bipartition to regions A and B the reduced density operator is ρ A ðtÞ ¼ tr B ψðtÞ j i ψðtÞ h j. The size of region A will be extensive, |A|~L d , in order to avoid the effects of measure concentration that becomes prominent when the ratio of subsystem Hilbert space sizes |B|/ |A| → ∞. Specifically, when that ratio grows the reduced density operator ρ A (t) for a typical state ψðtÞ j ibecomes increasingly closer to~1 18 . More precisely, tracing a random state over |B| ≫ |A| results in a spectrum of ρ A whose relative deviation from a flat one is~q −(|B|−|A|)/2 19 , and become negligible. Therefore, because we are interested in the influence of dynamics on S 2 , and not simple kinematic effects of Hilbert space sizes, we require |A| → ∞ in the TDL. For finite |A| the saturation value of S 2 would also be finite, so that one could not unambiguously differentiate between different powers in S 2~t α . To facilitate comparison of different d and L we will measure t in such units that one will generate a unit of entanglement in a unit of time, S 2 ð1Þ $ Oð1Þ, i.e., in the language of quantum circuits $ Oð1Þ gates connecting regions A and B are applied per unit of time. Compared to local Hamiltonian evolution this means a rescaling of time by L d−1 . None of our conclusions depends on the chosen time-units, i.e., on the values of potential crossover times. In all our numerical demonstrations shown in figures time is a dimensionless integer denoting the number of applied steps of a given quantum circuit.
Let us first argue why and how conservation of magnetization (charge) matters for the long-time behavior of S 2 . As we shall see, in the TDL S 2 is self-averaging (which is expected for generic, i.e., chaotic systems) and we will for simplicity focus on the purity IðtÞ :¼ 2 ÀS 2 ðtÞ ¼ trðρ 2 A ðtÞÞ. A non-rigorous intuitive meaning of the entropy is that it measures the effective number of the explored degrees of freedom needed to "describe" ρ A (t). For purity one can write I $ 1 N eff , where N eff $ 2 L eff is the effective Hilbert space size on which ρ A is supported, resulting in S 2 $ log 2 N eff $ L eff .
More quantitatively, the average purity I over all computational initial states is where D ðc k Þ k :¼ c k j i c k h j k is a basis of diagonal matrices (projectors) with c k 2 Z q labeling the local computational state, and We see that at large times what matters is the spreading of the reduced diagonal operators D ð c ! Þ A ðtÞ, specifically their Hilbert-Schmidt norm. In particular, the average purity gets contribution from all possible products of initial projectors, where at each site we have q different ones. The operator spreading in diffusive systems has a rich structure, having in general diffusive and ballistic features, see 20 for details. Operators with a large initial overlap with the conserved charge will have a hydrodynamic (power-law) tail, as well as some other operators (a trivial example is the associated conserved current). Such operators will tend to cause diffusive behavior of I. On the other hand, regardless of the diffusion, most operators will not exhibit any diffusive tails at long times. While in Eq. (2) one actually needs the reduced D ð c ! Þ A ðtÞ, regardless of details one can say that if the dynamics of all diagonal operators is diffusive, one expects diffusive I and S 2 $ ffiffi t p . For qubits, q = 2, there are just two local diagonal operators, 1 k and σ z k , and therefore if σ z k is diffusive one expects a long-time asymptotic growth S 2 $ ffiffi t p 13 . However, for higher dimensional qudits, q ≥ 3, the diagonal basis is spanned by q linearly independent diagonal operators, only one of which is the conserved operator (the local magnetization). While one might think that diffusive modes that contribute to purity decay as e À ffi ffi t p will due to their slow decay still dominate over non-diffusive ones, which decay as e −t , a simple counting argument shows that this is not to be expected. Namely, in a system of n qubits the number of diagonal operators that are products of only diffusive magnetization σ z k and the identity 1 k is 2 n , while the number of all other diagonal ones, that will in general be non-diffusive, is (q n − 2 n ). The diffusive contribution to I will then be $ 2 n e À ffi ffi t p while a non-diffusive one is~(q n − 2 n )e −t , so that in the TDL the non-diffusive contribution wins. Simply put, for higher q there are exponentially more non-diffusive operators than diffusive ones. For generic q ≥ 3 with only one conserved charge one therefore expects the asymptotic linear growth S 2~t .
How about the short-time behavior of S 2 (t)? We shall argue that it is, instead, always linear in time, even for qubits. Let us limit our discussion to qubits, q = 2, as for q ≥ 3 one anyway has linear growth S 2~t even at long times. For short-time behavior it is crucial to account for correlations spreading in a direction transversal to the boundary of dimension d − 1 and area A (d − 1) between regions A and B (in 2D A (d−1) is a circumference l, in 3D a true two-dimensional area A, in 1D the number of boundaries c, Table 1). Starting from a product initial state the dynamics tries to generate entanglement across the boundary. For local (n.-n.) interaction the natural first candidate sites to be entangled are all~L d − 1 n.-n. pairs lying on the boundary between A and B. Only after all those qubits are entangled can a slowing down due to diffusion in a transversal direction kick in. Let us be more specific, with a view on numerical demonstration. In our random quantum circuits we will apply L gates between random n.-n. qubits per unit of time. Such scaling is in-line with the mentioned units of time-probability that such a random n.-n. gate connects A and B is $ L dÀ1 L d ¼ 1 L , and therefore applying L of them means we will have~1 gates connecting regions A and B, and therefore, at least initially, generate one bit of entanglement in a unit of time. More precisely, the probability that a random gate connects A and B is A ðdÀ1Þ dL d , where the denominator dL d is the number of all nearest-neighbor bonds on a d dimensional square lattice. The initial growth of entanglement is therefore expected to be We expect this linear growth to hold for any S r , including r = 1. Such linear growth will continue until the time t 1~L d−1 at which S 2 (t 1 )~A (d−1) . After that one will crossover into the asymptotic diffusive growth We see that in higher spatial dimensions the region of diffusive growth is parametrically small, it lasts from t 1~L d − 1 till t ∞~L d+1 . Furthermore, in the TDL it is pushed to infinitely large values of entropy L d−1 ≲ S 2 ≲ L d and will be hard to observe. Qubits in d = 1 are rather special because the linear growth ends at S 2~L 0 = 1 (i.e., at short time t 1~1 ) and one gets S 2 $ ffiffi t p in the whole range of S 2 (and t). In short, in d = 1 the asymptotic $ ffiffi t p growth is "easy" to observe, while in d > 1 it is hard because it appears in the TDL only at infinitely large values of S 2 . Therefore the generic behavior after a quench from a product state is in d > 1 the linear growth (which is as fast as allowed by the Lieb-Robinson bound 21 ). Table 1 summarizes these findings.
Clifford circuits. It is always useful to take the simplest model, analytically or numerically, that displays the physics one wants to explore. A setting for which one can get exact results for the entanglement dynamics are so-called random quantum circuits 22 composed of a series of (random) local unitaries. Random circuits can be thought of as handy toy models of many-body physics but also as a useful theoretical concept called a unitary designs 23 . One of the first exact results was obtained by rewriting the dynamics of purity on average as a classical Markov process 24 , mapping it to a solvable quantum spin chain and getting an exact expression for the gap Δ or the decay rate 25 , i.e., entanglement speed 26 v E in modern language. For instance, for a circuit composed of a random 2-site unitaries applied to a random n.-n. pair of qubits in a chain with L sites, one gets 25 v E ¼ ð1 À 4 5 cos π L Þ 1 5 . If one would instead take a regular brick-wall pattern of applied gates, like in 27 , one instead has to calculate the gap of a product of Markovian matrices, obtaining a "multiplicative" form v E ¼ 2ln 5 4 cos π L , going in the TDL to v E 2ln 5 4 , as also calculated Table 1 Entanglement growth in diffusive lattice systems of linear size L in d spatial dimensions with q-level local Hilbert space (subsystem size is also ∝ L).
For qubits one has two regimes: linear growth for t < t 1 , and a slow square-root growth at later times, before the finite-size saturation S 2 (t∞)~L d is reached at t∞. Time units are such that dS 2 (0)/dt ≈ 1 for all L and d, A (d−1) is the area of a (d − 1)-dimensional boundary between the two subsystems (in dimension d = 3, 2, 1 this goes into area A, circumference l, and the number of boundaries c).
in 27,28 . Studies of random circuits have expanded in recent years, including U(1) conserving ones 20 , with many nice exact results, see e.g. refs. 20,[26][27][28][29][30][31] . They have been also notably used in a race toward quantum supremacy 32 . Let us check the above predictions for S 2 by numerical simulations of random circuits. In order to be able to simulate large systems we resort to the so-called Clifford circuits. For a q level system the local generalized Pauli operators X and Z are defined 33,34 as where ω: = e 2πi/q , and all additions are modulo q (the sign ⊕).
Generators of the local Pauli group are all q 2 products X v Z w , with v, w = 0, …, q − 1. The generalized Pauli group (GPG) on n sites is then formed by the tensor product of q 2 local Paulis, allowing also for all overall phases ω j . Due to ZX = ωXZ, a product of two members of the GPG is again in the GPG. The action of such Pauli operators on the computational basis states is simple, for instance, A Clifford circuit is a series of Clifford gates U, each of which preserves the GPG. That is, U j,k acting nontrivially on sites j and k maps a member of the GPG to another member of the GPG (instead of to a superposition of GPGs as for generic U). A sequence of such Clifford gates acting on a stabilizer state can be efficiently simulated (see Methods for details). A common choice of Clifford gates are the phase gate P j The dynamics of Clifford circuits therefore boils down to modular arithmetic 38 . They also form a unitary 2-design 39 (correctly reproduce Haar averages over all 2nd order polynomials in ρ A (t → ∞), e.g., a purity), with the same convergence behavior as generic random circuits. Therefore, in absence of conservation laws S 2 (t) behaves similarly for Clifford and for generic random circuits. So far Clifford circuits have been extensively studied in quantum information, but not so much in condensed matter or statistical physics. The reason being that their dynamics is typically either ballistic or localized 40 (fluctuations though can exhibit interesting behavior 26 ). We shall study a new class of random Clifford circuits that conserve magnetization and whose dynamics is diffusive. By looking at random circuits we are also able to focus exclusively on the role of the U(1) symmetry without any stray effects caused by other conservation laws (e.g., conservation of energy).
Numerical verification. Let us first focus on qubits. For qubits the elements of the local GPG are just the ordinary Pauli matrices fσ x ; σ y ; σ z ; 1g. To preserve the total magnetization our Clifford circuit consists of applying the XY gate U XY :¼ expðÀi π 4 ðσ x j σ x k þ σ y j σ y k ÞÞ to a randomly selected n.-n. pair of sites on a d dimensional square lattice. It is easy to verify that U y XY 1 j σ z k U XY ¼ σ z j 1 k and U y XY σ z j 1 k U XY ¼ 1 j σ z k , and therefore the total magnetization σ z j þ σ z k is conserved. It also implies that a pair of oppositely polarized spins is exchanged, U XY "# j i ¼ #" j i. Because the pair (j, k) is chosen at random, it is also immediately clear that the dynamics of magnetization is diffusive, e.g., starting from a domain wall initial state # #" " j ithe average profile at time t can be expressed exactly in terms of binomial probabilities, that can be approximated in the large-t limit by the error function (see Fig. 1 for an explicit numerical demonstration).
Starting with the initial state ψ j i $ ð " j i þ # j iÞ n stabilized by g j ¼ σ x j , we can simulate our Clifford circuit for thousands of qubits up-to very long times, despite the entanglement eventually being a volume-law S 2~L d . Entanglement calculation is simplified by the fact that the state at any time is composed of an integer number M of generalized EPR pairs, 1 stabilized by two generators X 1 X 2 and Z 1 Z À1 2 . Eigenvalues of the reduced density operator ρ A (t) are all equal, and using the base-q logarithm one has S r = M for all r. In this respect Clifford circuits are special, however, their dynamics of S 2 will be generic and consistent with the presented theory. In addition, we will also numerically demonstrate that a similar behavior is obtained also for non-Clifford circuits where the spectrum of the reduced density operator is not flat. Therefore, while one can not extract the difference between different S r from Clifford circuit simulations (in particular that S 1 can behave differently), or use any finer measures of complexity that involve the individual spectral components, e.g. ref. 41 , we argue that they do result in generic behavior of S 2 . This is in-line with the fact that while Clifford circuits are not universal, already very small modifications, see e.g. refs. 42,43 (that might not influence many quantities), do result in universal behavior, e.g. universal quantum computation. For another solvable evolution that also results in a flat spectrum of ρ A see ref. 44 .
In Fig. 2a-c we show S 2 for 1D, 2D, and 3D lattice, and for different bipartite splitting of n spins into regions A and B. For 1D we see that the asymptotic growth is S 2 ¼ c ffiffiffiffiffiffiffiffiffi ffi 2t=π p with c being the number of boundaries between A and B (c = 1 for a half-cut, and c = 2 for the middle-1 3 cut). We also observe that at short times t ≲ 10 the growth is a bit faster than diffusive. This means that in small systems L~30 (being a typical maximal size amenable to other methods) it would be very difficult to see the true asymptotic growth over a significant range of times. In 2D we show data only for the case where the region A is the middle-1 3 part of the full square lattice with n = L × L qubits, as this bipartition gives a clearer transition between linear and diffusive growth (see "Methods" for the half-cut data). Numerics confirms the short-time growth given by Eq. (3) without any additional prefactors (A (1) = 4L/3). We also note that to see the asymptotic growth $ ffiffiffiffi ffi Lt p one needs fairly large systems; even for L = 252 one can see only about one decade in time of S 2 $ ffiffi t p , while on the other hand three decades of S 2~t . In 3D the situation is even less favorable for slow asymptotic diffusive growth. Nevertheless, in the more favorable half-cut bipartition we can see a transition from the short-time S 2 ¼ A ð2Þ 3L 2 t to the long-time S 2 $ ffiffiffiffiffiffi ffi L 2 t p . Again, for t < t 1 the linear growth Eq. (3) has no additional prefactors (for a half-cut A (2) = L 2 , for the middle-1 3 cut A ð2Þ ¼ 2 3 L 2 ). For the growth. Finally, we show entanglement profiles (Fig. 3a), i.e. S 2 for a bipartite cut with region A being the first k spins. Also shown are the fluctuations of S 2 between different circuit realizations (Fig 3b), showing that the relative fluctuations scale as σðS 2 Þ=S 2 $ 1= ffiffiffiffi ffi S 2 p , and therefore in the TDL at large times dynamics is self-averaging. It suffices to look at a single random circuit realization. In "Methods" we also show data for more complicated Clifford gates than U XY , leading to similar results.
We also check the case of qudits with q > 2, also studied in 13 . To that end we simulate a qutrit chain (q = 3, i.e. spin-1 particles) where the local diagonal basis is spanned by fZ k ; Z 2 k ¼ Z À1 k ; Z 3 k ¼ 1g, and we take the initial state stabilized by generators g j = X j . Taking a Clifford circuit with the n.-n. gate being U D = H 2 CNOT 21 CNOT 12 CNOT 12 H 1 , which gives rise to diffusive conservative dynamics of both diagonal matrices, Fig. 4). For further qutrit numerics see "Methods", including circuits with non-Clifford gates.
To get a generic case in which not all diagonal operators are conserved we use local dimension q = 4 visualized as a qubit ladder (n = 2L) with the local space corresponding to one rung. We apply L steps per unit of time, each step consisting of: the gate U ZZ ¼ i expðÀi π 2 σ z 1 τ z 2 Þ applied to a random rung (σ α are Pauli matrices on the upper and τ α on the lower leg), and either magnetization-preserving U XY on a random bond in the upper leg, or a non-conserving U G ¼ CNOT 12 H 1 expði π 4 τ z 2 Þ on a random bond of the lower leg. Magnetization is conserved only in the upper leg and one therefore expects generic behavior with S 2~t . This is indeed observed in Fig. 4. We can also see that for times less than ≈ 10 2 slower growth is observed in which diffusive dynamics competes with increasingly dominating non-diffusive dynamics of other operators, and therefore, once again, one needs large systems with L ≳ 100 in order to see the true linear asymptotic growth. We contrast this linear growth with a special case: if we instead of U G apply U XY also on the lower leg, such that magnetization on both legs is conserved, one again gets a non-generic S 2 $ ffiffi t p (orange dashed curve in Fig. 4). We remark that this latter case of using U XY on both legs corresponds to a Trotterized dynamics of the Hubbard chain (using Jordan-Wigner transformation the upper leg represents spin-up fermions, the lower spin-down, U XY is hopping, while U ZZ is the on-site interaction). In "Methods" we show further ladder examples.

Discussion
We have presented a theory of the Rényi entropy growth in lattice systems that conserve the total magnetization due to U(1) symmetry. We show that in general qubit systems the entanglement grows linearly in time until, at an area-law value of S 2 , a crossover to slower square-root growth happens. In 1D qubit systems the diffusive ffiffi t p growth is generic because the crossover happens already at small values of entanglement S 2~1 , while in higher dimensions the regime of such growth is in the thermodynamic limit pushed to infinitely large values of S 2 (at any finite value of S 2 the growth is linear). For lattice systems with more than 2 local levels (spin-s particles with s ≥ 1) and a single conserved charge, non-diffusive degrees dominate and one expects  the linear growth. An exception is a situation where the dynamics of all diagonal operators is diffusive. Two-level systems in 1D are special because a single diffusive charge exhausts all non-trivial local projectors. Entanglement growth can therefore distinguish both the spatial dimensionality as well as the size of the local Hilbert space, with the influence of a diffusive charge diminishing when either of the two increases. It would be interesting to generalize our results to other transport types beyond diffusion, an obvious conjecture in 1D is that asymptotically one will have S r > 1~t 1/z , with z being a dynamical transport exponent. An interesting question is the influence of more complicated symmetries than U(1), as well as the presence of multiple symmetries. Energy conservation (i.e., time-translation symmetry) which comes automatically in autonomous Hamiltonian systems, and which was not discussed specifically, is an important example. While diffusive energy transport likely plays a similar role, it carries few technical complications, for instance, it is a nontrivial problem to establish diffusion of energy in the first place. A class of systems not touched upon are integrable systems where one has an extensive set of local conserved quantities. Those will typically be ballistic, however it needs not be so. The question is, can such non-ballistic modes influence the entropy growth in some generic fashion? In the present work we focus on S 2 , and conjecture that all integer S r > 1 behave similarly; there could however be non-trivial time-scales connected with the index r. von Neumann entropy S 1 is special 13 , and one could more generally study how the whole average eigenvalue spectrum converges to that of a random state 19 .
A promising direction is also employing introduced nontrivial Clifford circuits to further explore the many-body physics, and, more generally, to understand their dynamical properties and in which aspects are they different than those of the many-body Hamiltonian systems.

Methods
Simulating Clifford circuits. Evolution of states under conjugation by Clifford circuits is not done by updating each computational basis state-that would be inefficient for highly entangled states-but rather by a stabilizer formalism 35 . A state ψ j i on n qubits is called a stabilizer state if it is a unique joint eigenstate with eigenvalue 1 of n independent stabilizer generators g j from the GPG. For qubits one can obtain it as a product of projectors, ψ j i ψ h j ¼ Π j ð1 þ g j Þ=2, for qutrits one has ψ j i ψ h j ¼ Π j ð1 þ g j þ g 2 j Þ=3. Here lies the advantage of Clifford circuits. Instead of updating the state ψ j i one instead updates each generator g j , whose number is always n and which will remain elements of the GPG 33,35 . Performing one gate, i.e., updating all stabilizers, takes Oðn 2 Þ operations. Entanglement and state overlaps can also be calculated efficiently 36,37 .
Additional data for 1D and 2D Clifford systems. In the main text we used the XY gate in qubit (q = 2) systems to demonstrate diffusive growth of S 2 . Such gate is quadratic in fermionic operators. Here we show that using more complicated Clifford gates (which in particular are not quadratic), similarly as has been done for q > 2, leads to similar results. In Fig. 5 we show data for S 2 for three different types of Clifford circuits. One that uses only the XY gate (similar data as in Fig. 2 a), one with a gate U XYP ¼ U XY exp½Ài π 4 ðσ z 1 þ σ z 2 Þ, and one with a gate U G ¼ CNOT 12 H 1 expði π 4 σ z 2 Þ that does not conserve the magnetization. In Fig. 6 we show data for the same 2D diffusive qubit Clifford circuit utilizing U XY gate as in Fig. 2b, but for a half-cut bipartition. We can see that the agreement with theoretical short-time as well as long-time prediction (Table 1) is good. The short-time growth is S 2 ≈ 0.5t, where 0:5 ¼ l 2L with l = L, whereas the asymptotic growth goes into S 2 ffiffiffiffiffiffiffi 2Lt p (no fitting parameters).
Additional data for non-Clifford qutrit systems (spin S = 1). In the main text we demonstrated that the dynamics given by the Clifford qutrit gate U D , which conserves both the non-trivial diagonal operators Z 1 + Z 2 and Z 2 1 þ Z 2 2 , results in a diffusive asymptotic growth of S 2 (Fig. 4). Instead of the two "Clifford"-basis diagonal operators Z j and Z 2 j we can also use the language of spin S = 1 particles: there the two non-trivial diagonal operators are S z j ¼ diagð1; 0; À1Þ andS Note that U D is up-to phases equal to the spin-1 SWAP gate U SWAP ¼ expðÀi π 2 ½ðS 1 Á S 2 Þ 2 þ ðS 1 Á S 2 Þ þ 2 Á 1Þ, U D U y SWAP ¼ diagð1; 1; 1; 1; ω; ω 2 ; 1; ω 2 ; ωÞ.
We are going to show additional data for a number of spin-1 quantum circuits that are not of the Clifford type. We will always use a half-cut bipartition and the same initial state as used for Clifford circuits, that is a uniform superposition of all computational states, ψð0Þ $ ð 1 j i þ 0 j i þ À1 j iÞ L . The shown S 2 is the average one over between 10 3 (small L) and 20 (for largest L = 21) circuit realizations. The aim is to further shed light on the fact that we expect the asymptotic S 2 $ ffiffi t p growth only if all diagonal operators are conserved and diffusive. For q = 3 this means that both S z andS z should be conserved (alternatively, both Z and Z 2 ).
We first check the evolution using a 2-site gate U I that is a concatenation of the diffusive U D we already used in the main text and the isotropic gate, that is The gate U I therefore conserves only S z 1 þ S z 2 , the dynamics of which is diffusive due to the spatial randomness (a gate is applied to a random n.n. bond). We therefore expect the asymptotic growth of S 2 to be linear for U I despite diffusive total magnetization. The results are shown in Fig. 7a, b. We immediately have to remark that a drawback of non-Clifford circuits is that only very small systems can be simulated, and correspondingly the reachable times are far from the asymptotic ones.
For comparison we also show data for the Clifford circuit with U D and L = 2000 (the same data as in Fig. 4), as well as for L = 10 and L = 20. Because the initial rates of the entropy production are different for U D and U I we multiply the times for the U D data by 1.4 so that the curves overlap at short times. We can see that upto times t ≈ 6 the two evolutions result in the same growth of S 2 (one could fit S 2 ≈ 0.59t 0.76 ); for instance, for L = 10 qutrits it is hard to claim any difference between U D and U I . After t ≈ 10 though deviations start to appear: the Clifford case U D that conserves both diagonal operators starts to converge to slower S 2 ffiffi t p growth, whereas U I that conserves only the total S z starts to grow faster. This is furthermore seen also in the dependence of the scaling exponent in S 2~t α on time. Due to taking a numerical derivative the data for α is much more noisy (particularly for L = 20 where the ensemble size is 100), however one can nevertheless see a clear difference between the two cases (in-line with observation in Fig. 7 a). For the  non-Clifford circuit α increases with time, while for the Clifford case it decreases. While from such short-time data it is impossible to make a definite claim about the asymptotic growth, what we observe is compatible with the asymptotics S 2~t for U I . If the Clifford case, where we can simulate large systems, is any indication of the required sizes necessary to reach the asymptotics, we can say that likely about five times larger systems would be required to really see the asymptotic linear growth for U I (for the Clifford data in Fig. 4, we can see that one converges to S 2 $ ffiffi t p only at t ≈ 10 3 where S 2 ≈ 10 2 ).
In Fig. 8 we show data for further non-Clifford random circuits. We show results , we see that the growth is much slower, like S 2~t 0.7 at short times. While it is hard to make any asymptotic claims about S 2 $ ffiffi t p based on such numerics (exact numerics for larger systems gets hampered by memory requirements; the Hilbert space size for L = 21 qutrits is about the same as for ≈33 qubits), what is very distinct is that the exponent is very different in (a) and (b) despite a very similar 2-site gate; the only difference between the two is the conservation ofS z 1 þS z 2 . Finally, as a third example we show the anisotropic XXZlike gate U XXZ ¼ expðÀi π 3 ½S x 1 S x 2 þ S y 1 S y 2 þ 1:5S z 1 S z 2 Þ that again conserves only the total magnetization, and therefore one has S 2~t visible already at short times.
Data for ladder systems (q = 4). Here we show further data supporting the claim that for ladders the growth of S 2 is generically linear in time. We simulate Clifford ladders with a large number of rungs L = 10000 using different 2-site gates. On legs, upper or lower, we apply either the already seen U XY ¼ expðÀi π 4 ðσ x j σ x k þ σ y j σ y k ÞÞ, or U G ¼ CNOT 12 H 1 expði π 4 σ z 2 Þ. On the rungs we use either The protocol is always the same: at each step we apply one of the leg gates on a random bond on either the upper or the lower leg, and a rung gate on an independent random rung. The type of the gate applied on rungs as well as on the upper and lower leg is held fixed, so that one can get 8 different protocols out of the 4 mentioned gates. The gate U XY conserves the total magnetization on the respective leg on which it acts, while U G does not. Namely, U y The rung gate U ZZ does not break conservation of σ z 1 þ σ z 2 , nor of τ z 1 þ τ z 2 because one has U y ZZ σ z 1 1 2 U ZZ ¼ σ z 1 1 2 , and U y ZZ 1 1 τ z 2 U ZZ ¼ 1 1 τ z 2 (as well as U y ZZ σ x 1 1 2 U ZZ ¼ Àσ x 1 1 2 , U y ZZ 1 1 τ x 2 U ZZ ¼ À1 1 τ x 2 ). The gate U ZZ though does introduces non-trivial phases ±1 in the dynamics of Pauli x and y matrices. The gate U Sxy on the other hand preserves conservation of magnetization only on the upper leg, U y Sxy σ z 1 1 2 U Sxy ¼ σ z 1 1 2 , U y Sxy σ x 1 1 2 U Sxy ¼ Àσ y 1 1 2 , while it breaks conservation on the lower leg, U y Sxy 1 1 τ z 2 U Sxy ¼ 1 1 τ y 2 , U y Sxy 1 1 τ x 2 U Sxy ¼ 1 1 τ z 2 . In Fig. 9 we show results of numerical simulation for different protocols. Taking the (XY) − (ZZ) − (XY) protocol where the leg gates U XY as well as the rung gates U ZZ conserve magnetization on the upper and the lower leg, dynamics of all diagonal operators is diffusive and one has S 2 ffiffi t p . The same data for L = 4000 has been already shown in Fig. 4. We can break conservation of magnetization on the lower leg by using U G , which as we can see results in the asymptotic growth S 2 ≍ t (red curve in Fig. 9, the same data as in Fig. 4). We can however break the conservation on the lower leg also by changing the rung gate to U Sxy . This is illustrated by the protocol (XY) − (Sxy) − (XY), which again results in S 2 ≍ t. Note  that here the dynamics along the two legs is purely diffusive-the gate U XY is used on both legs-it is only the non-trivial rung dynamics that breaks one U(1) symmetry and causes the asymptotic linear growth of S 2 (similar result would be obtained also if at each step of the protocol the U XY gate would be applied simultaneously to a pair of upper-and lower-leg bonds forming a local plaquette with two rungs). Finally, using conservation breaking U G on the lower leg as well as the rung U Sxy , one again has S 2 ≍ t.

Data availability
Data is available upon reasonable request from the author.