Entanglement growth in diffusive systems

Entanglement 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¹. Even if not easily measurable², 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³. 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⁴. 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 many-body-localized phase⁷, one of the distinguishing features of which is slow logarithmic-in-time growth of entanglement⁸.

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⁹. 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¹⁰. 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¹³ (focusing on diffusive 1D systems with conserved charge) that the entanglement, as quantified by the Rényi  entropy

grows in fact as \({S}_{2} \sim \sqrt{t}\) starting from a generic separable initial state, instead of the “expected” S₂ ~ 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” \(\sim {2}^{{S}_{2}}\) to grow only as \(\sim {b}^{\sqrt{t}}\), 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⁴ in 2D, and 48 × 48 × 48 ≈ 10⁵ in 3D. While we confirm¹³, the main finding is that in higher d and q the asymptotic \(\sqrt{t}\) 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 \(\sim \sqrt{t}\), it starts in d > 1 only at a time when the entropy already becomes extensively large, S₂ ~ L^d−1, in other words, in the TDL the S₂ 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.

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₁, are expected to behave in the same way, whereas for diffusive systems the S₁ (which we don’t discuss) can behave differently¹³. We will mostly focus on S₂ as a representative case of S_r>1. We remark that sometimes S₂, rather than e.g. S₁, is a more relevant quantity¹⁶, and is furthermore also easier to measure¹⁷. Using a bipartition to regions A and B the reduced density operator is \({\rho }_{{\rm{A}}}(t)={{\rm{t}}r}_{{\rm{B}}}\left|\psi (t)\right\rangle \left\langle \psi (t)\right|\). 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 \(\left|\psi (t)\right\rangle\) becomes increasingly closer to ~ \({\mathbb{1}}\)¹⁸. 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₂, and not simple kinematic effects of Hilbert space sizes, we require ∣A∣ → ∞ in the TDL. For finite ∣A∣ the saturation value of S₂ would also be finite, so that one could not unambiguously differentiate between different powers in S₂ ~ 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) \sim {\mathcal{O}}(1)\), i.e., in the language of quantum circuits \(\sim {\mathcal{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₂. As we shall see, in the TDL S₂ 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)}={\rm{tr}}({\rho }_{{\rm{A}}}^{2}(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 \sim \frac{1}{{N}_{{\rm{eff}}}}\), where \({N}_{{\rm{eff}}} \sim {2}^{{L}_{{\rm{eff}}}}\) is the effective Hilbert space size on which ρ_A is supported, resulting in \({S}_{2} \sim {\mathrm{log}\,}_{2}{N}_{{\rm{eff}}} \sim {L}_{{\rm{eff}}}\).

More quantitatively, the average purity I over all computational initial states is

where \({D}_{k}^{({c}_{k})}:=\left|{c}_{k}\right\rangle {\left\langle {c}_{k}\right|}_{k}\) is a basis of diagonal matrices (projectors) with \({c}_{k}\in {{\mathbb{Z}}}^{q}\) labeling the local computational state, and \({D}_{{\rm{A}}}^{(\overrightarrow{{\bf{c}}})}(t)={{\rm{t}}r}_{{\rm{B}}}[{U}^{t}{D}_{1}^{({c}_{1})}\otimes \cdots \otimes {D}_{n}^{({c}_{n})}{U}^{-t}]\). We see that at large times what matters is the spreading of the reduced diagonal operators \({D}_{{\rm{A}}}^{(\scriptstyle\overrightarrow{{\bf{c}}})}(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²⁰ 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 \(\bar{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}_{{\rm{A}}}^{(\overrightarrow{{\bf{c}}})}(t)\), regardless of details one can say that if the dynamics of all diagonal operators is diffusive, one expects diffusive \(\bar{I}\) and \({S}_{2} \sim \sqrt{t}\). For qubits, q = 2, there are just two local diagonal operators, 1_k and \({\sigma }_{k}^{{\rm{z}}}\), and therefore if \({\sigma }_{k}^{{\rm{z}}}\) is diffusive one expects a long-time asymptotic growth \({S}_{2} \sim \sqrt{t}\)¹³. 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 \({{\rm{e}}}^{-\sqrt{t}}\) 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 \({\sigma }_{k}^{{\rm{z}}}\) and the identity 1_k is 2ⁿ, while the number of all other diagonal ones, that will in general be non-diffusive, is (qⁿ − 2ⁿ). The diffusive contribution to I will then be \(\sim {2}^{n}{{\rm{e}}}^{-\sqrt{t}}\) while a non-diffusive one is ~(qⁿ − 2ⁿ)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₂ ~ t.

How about the short-time behavior of S₂(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₂ ~ 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 \(\sim \frac{{L}^{d-1}}{{L}^{d}}=\frac{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 \(\frac{{A}^{(d-1)}}{d{L}^{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₁ ~ L^d−1 at which S₂(t₁) ~ A^(d−1). After that one will crossover into the asymptotic diffusive growth \({S}_{2} \sim \sqrt{t}\), until at t_∞ ~ L^d+1 a finite-size saturation value S₂(t_∞) ~ L^d is reached.

We see that in higher spatial dimensions the region of diffusive growth is parametrically small, it lasts from t₁ ~ 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₂ ≲ L^d and will be hard to observe. Qubits in d = 1 are rather special because the linear growth ends at S₂ ~ L⁰ = 1 (i.e., at short time t₁ ~ 1) and one gets \({S}_{2} \sim \sqrt{t}\) in the whole range of S₂ (and t). In short, in d = 1 the asymptotic \(\sim \sqrt{t}\) 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₂. 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²¹). 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²² 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²³. One of the first exact results was obtained by rewriting the dynamics of purity on average as a classical Markov process²⁴, mapping it to a solvable quantum spin chain and getting an exact expression for the gap Δ or the decay rate²⁵, i.e., entanglement speed²⁶ 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²⁵ \({v}_{{\rm{E}}}=(1-\frac{4}{5}\cos \frac{\pi }{L})\asymp \frac{1}{5}\). If one would instead take a regular brick-wall pattern of applied gates, like in²⁷, one instead has to calculate the gap of a product of Markovian matrices, obtaining a “multiplicative” form \({v}_{{\rm{E}}}=2\mathrm{ln}\,\frac{5}{4\cos \frac{\pi }{L}}\), going in the TDL to \({v}_{{\rm{E}}}\asymp 2\mathrm{ln}\,\frac{5}{4}\), as also calculated in^27,28. Studies of random circuits have expanded in recent years, including U(1) conserving ones²⁰, 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³².

Let us check the above predictions for S₂ 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² products X^vZ^w, with v, w = 0, …, q − 1. The generalized Pauli group (GPG) on n sites is then formed by the tensor product of q² 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, \({X}^{{x}_{1}}{Z}^{{z}_{1}}\otimes \cdots \otimes {X}^{{x}_{n}}{Z}^{{z}_{n}}\left|{\bf{a}}\right\rangle ={\omega }^{{\bf{z}}\cdot {\bf{x}}}\left|{\bf{a}}\oplus \bf{x}\right\rangle\).

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 \({\rm{P}}\left|j\right\rangle ={\omega }^{j(j-1)/2}\left|j\right\rangle\), the Hadamard gate \({\rm{H}}\left|j\right\rangle =\frac{1}{\sqrt{q}}{\sum }_{k}{\omega }^{kj}\left|k\right\rangle\), and a 2-qudit controlled-NOT gate \({{\rm{C}}NOT}_{12}\left|j,k\right\rangle =\left|j,k\oplus j\right\rangle\). The dynamics of Clifford circuits therefore boils down to modular arithmetic³⁸. They also form a unitary 2-design³⁹ (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₂(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⁴⁰ (fluctuations though can exhibit interesting behavior²⁶). 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 \(\{{\sigma }_{,}^{{\rm{x}}}{\sigma }_{,}^{{\rm{y}}}{\sigma }_{,}^{{\rm{z}}}{\mathbb{1}}\}\). To preserve the total magnetization our Clifford circuit consists of applying the XY gate \({U}_{{\rm{XY}}}:=\exp (-{\rm{i}}\frac{\pi }{4}({\sigma }_{j}^{{\rm{x}}}{\sigma }_{k}^{{\rm{x}}}+{\sigma }_{j}^{{\rm{y}}}{\sigma }_{k}^{{\rm{y}}}))\) to a randomly selected n.-n. pair of sites on a d dimensional square lattice. It is easy to verify that \({U}_{{\rm{XY}}}^{\dagger }{{\mathbb{1}}}_{j}{\sigma }_{k}^{{\rm{z}}}{U}_{{\rm{XY}}}={\sigma }_{j}^{{\rm{z}}}{{\mathbb{1}}}_{k}\) and \({U}_{{\rm{XY}}}^{\dagger }{\sigma }_{j}^{{\rm{z}}}{{\mathbb{1}}}_{k}{U}_{{\rm{XY}}}={{\mathbb{1}}}_{j}{\sigma }_{k}^{{\rm{z}}}\), and therefore the total magnetization \({\sigma }_{j}^{{\rm{z}}}+{\sigma }_{k}^{{\rm{z}}}\) is conserved. It also implies that a pair of oppositely polarized spins is exchanged, \({U}_{{\rm{XY}}}\left|\uparrow \downarrow \right\rangle =\left|\downarrow \uparrow \right\rangle\). 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 \(\left|\,\downarrow \ldots \downarrow \uparrow \ldots \uparrow \right\rangle\) the 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).

The inset shows the domain wall profile z_k, while the main plot show a diffusive growth of transferred magnetization ΔZ across the domain wall (averaged over 10³ circuit realizations and for L = 1000).

Full size image

Starting with the initial state \(\left|\psi \right\rangle \sim {(\left|\uparrow \right\rangle +\left|\downarrow \right\rangle )}^{\otimes n}\) stabilized by \({g}_{j}={\sigma }_{j}^{{\rm{x}}}\), we can simulate our Clifford circuit for thousands of qubits up-to very long times, despite the entanglement eventually being a volume-law S₂ ~ 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, \(\frac{1}{\sqrt{q}}{\sum }_{j}{\left|j\right\rangle }_{1}\otimes {\left|j\right\rangle }_{2}\), stabilized by two generators X₁X₂ and \({Z}_{1}{Z}_{2}^{-1}\). 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₂ 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₁ can behave differently), or use any finer measures of complexity that involve the individual spectral components, e.g. ref. ⁴¹, we argue that they do result in generic behavior of S₂. 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. ⁴⁴.

In Fig. 2a–c we show S₂ 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\sqrt{2t/\pi }\) with c being the number of boundaries between A and B (c  = 1 for a half-cut, and c = 2 for the middle-\(\frac{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-\(\frac{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⁽¹⁾ = 4L/3). We also note that to see the asymptotic growth \(\sim \sqrt{Lt}\) one needs fairly large systems; even for L = 252 one can see only about one decade in time of \({S}_{2} \sim \sqrt{t}\), while on the other hand three decades of S₂ ~ 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}=\frac{{A}^{(2)}}{3{L}^{2}}t\) to the long-time \({S}_{2} \sim \sqrt{{L}^{2}t}\). Again, for t < t₁ the linear growth Eq. (3) has no additional prefactors (for a half-cut A⁽²⁾ = L², for the middle-\(\frac{1}{3}\) cut \({A}^{(2)}=\frac{2}{3}{L}^{2}\)). For the middle-\(\frac{1}{3}\) cut, despite a large number of qubits, n = 48³ = 110,592, one can barely hint the eventual \(\sim\sqrt{t}\) growth. Finally, we show entanglement profiles (Fig. 3a), i.e. S₂ for a bipartite cut with region A being the first k spins. Also shown are the fluctuations of S₂ between different circuit realizations (Fig 3b), showing that the relative fluctuations scale as \(\sigma ({S}_{2})/{S}_{2} \sim 1/\sqrt{{S}_{2}}\), 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.

a In 1D (solid curves are for the middle-\(\frac{1}{3}\) bipartition), 2D in b, and 3D in c (solid blue, red and green curves are for the half-cut bipartition).

Full size image

a Shows profiles at different times, and b fluctuations, all in diffusive 1D Clifford circuit with L = 500. In a we also show standard deviation (gray shading) and one realization at two selected times.

Full size image

We also check the case of qudits with q > 2, also studied in¹³. To that end we simulate a qutrit chain (q = 3, i.e. spin-1 particles) where the local diagonal basis is spanned by \(\{{Z}_{k},{Z}_{k}^{2}={Z}_{k}^{-1},{Z}_{k}^{3}={\mathbb{1}}\}\), 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₂ CNOT₂₁CNOT₁₂CNOT₁₂ H₁, which gives rise to diffusive conservative dynamics of both diagonal matrices, \({U}_{{\rm{D}}}^{\dagger }{{\mathbb{1}}}_{j}{Z}_{k}{U}_{{\rm{D}}}={Z}_{j}{{\mathbb{1}}}_{k}\), \({U}_{{\rm{D}}}^{\dagger }{Z}_{j}{{\mathbb{1}}}_{k}{U}_{{\rm{D}}}={{\mathbb{1}}}_{j}{Z}_{k}\), and \({U}_{{\rm{D}}}^{\dagger }{{\mathbb{1}}}_{j}{Z}_{k}^{2}{U}_{{\rm{D}}}={Z}_{j}^{2}{{\mathbb{1}}}_{k}\), \({U}_{{\rm{D}}}^{\dagger }{Z}_{j}^{2}{{\mathbb{1}}}_{k}{U}_{{\rm{D}}}={{\mathbb{1}}}_{j}{Z}_{k}^{2}\), we observe the expected diffusive \({S}_{2} \sim \sqrt{t}\) (Fig. 4). For further qutrit numerics see “Methods”, including circuits with non-Clifford gates.

Solid curves (blue, green, red) are for the ladder (q = 4) with diffusive dynamics only on the upper leg, while the dashed orange curve (labeled “full L = 4000”) is for the ladder with U(1) conservation on both legs. Solid olive curve is for the qutrit (q = 3) chain with diffusive dynamics of all diagonal operators.

Full size image

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}_{{\rm{Z}}Z}={\rm{i}}\exp (-{\rm{i}}\frac{\pi }{2}{\sigma }_{1}^{{\rm{z}}}{\tau }_{2}^{{\rm{z}}})\) 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}_{{\rm{G}}}={{\rm{C}}NOT}_{12}\ {{\rm{H}}}_{1}\ \exp ({\rm{i}}\frac{\pi }{4}{\tau }_{2}^{{\rm{z}}})\) on a random bond of the lower leg. Magnetization is conserved only in the upper leg and one therefore expects generic behavior with S₂ ~ t. This is indeed observed in Fig. 4. We can also see that for times less than ≈ 10² 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} \sim \sqrt{t}\) (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.

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₂, a crossover to slower square-root growth happens. In 1D qubit systems the diffusive \(\sqrt{t}\) growth is generic because the crossover happens already at small values of entanglement S₂ ~ 1, while in higher dimensions the regime of such growth is in the thermodynamic limit pushed to infinitely large values of S₂ (at any finite value of S₂ 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₂, 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₁ is special¹³, and one could more generally study how the whole average eigenvalue spectrum converges to that of a random state¹⁹.

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.

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³⁵. A state \(\left|\psi \right\rangle\) 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, \(\left|\psi \right\rangle \left\langle \psi \right|={\Pi }_{j}({\mathbb{1}}+{g}_{j})/2\), for qutrits one has \(\left|\psi \right\rangle \left\langle \psi \right|={\Pi }_{j}({\mathbb{1}}+{g}_{j}+{g}_{j}^{2})/3\). Here lies the advantage of Clifford circuits. Instead of updating the state \(\left|\psi \right\rangle\) 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 \({\mathcal{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₂. 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₂ 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}_{{\rm{X}}YP}={U}_{{\rm{XY}}}\exp [-{\rm{i}}\frac{\pi }{4}({\sigma }_{1}^{{\rm{z}}}+{\sigma }_{2}^{{\rm{z}}})]\), and one with a gate \({U}_{{\rm{G}}}={{\rm{C}}NOT}_{12} {{\rm{H}}}_{1}\exp ({\rm{i}}\frac{\pi }{4}{\sigma }_{2}^{{\rm{z}}})\) that does not conserve the magnetization.

Data for magnetization conserving U_XY and U_XYP almost overlap, both growing diffusively. All is for L = 1002 and averaged over 100 realizations.

Full size image

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₂ ≈ 0.5t, where \(0.5=\frac{l}{2L}\) with l = L, whereas the asymptotic growth goes into \({S}_{2}\asymp \sqrt{2Lt}\) (no fitting parameters).

Initial linear growth transitions at t_1 into asymptotic diffusive growth.

Full size image

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₁ + Z₂ and \({Z}_{1}^{2}+{Z}_{2}^{2}\), results in a diffusive asymptotic growth of S₂ (Fig. 4). Instead of the two "Clifford”-basis diagonal operators Z_j and \({Z}_{j}^{2}\) we can also use the language of spin S = 1 particles: there the two non-trivial diagonal operators are \({S}_{j}^{{\rm{z}}}={\rm{d}}{\mathrm{iag}}(1,0,-1)\) and \({\tilde{S}}_{j}^{{\rm{z}}}=3{({S}_{j}^{{\rm{z}}})}^{2}-2\cdot {{\mathbb{1}}}_{j}={\rm{d}}{\mathrm{iag}}(1,-2,1)\). The gate U_D of course also conserves those two, \({U}_{{\rm{D}}}^{\dagger }({S}_{1}^{{\rm{z}}}+{S}_{2}^{{\rm{z}}}){U}_{{\rm{D}}}={S}_{1}^{{\rm{z}}}+{S}_{2}^{{\rm{z}}}\), \({U}_{{\rm{D}}}^{\dagger }({\tilde{S}}_{1}^{{\rm{z}}}+{\tilde{S}}_{2}^{{\rm{z}}}){U}_{{\rm{D}}}={\tilde{S}}_{1}^{{\rm{z}}}+{\tilde{S}}_{2}^{{\rm{z}}}\). Note that U_D is up-to phases equal to the spin-1 SWAP gate \({U}_{{\rm{S}}WAP}=\exp (-{\rm{i}}\frac{\pi }{2}[{({{\bf{S}}}_{1}\cdot {{\bf{S}}}_{2})}^{2}+({{\bf{S}}}_{1}\cdot {{\bf{S}}}_{2})+2\cdot {\mathbb{1}}])\), \({U}_{{\rm{D}}}{U}_{{\rm{S}}WAP}^{\dagger }={\rm{d}}{\mathrm{iag}}(1,1,1,1,\omega ,{\omega }^{2},1,{\omega }^{2},\omega )\).

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, \(\psi (0) \sim {(\left|1\right\rangle +\left|0\right\rangle +\left|\text{-}1\right\rangle )}^{\otimes L}\). The shown S₂ is the average one over between 10³ (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} \sim \sqrt{t}\) growth only if all diagonal operators are conserved and diffusive. For q = 3 this means that both S^z and \({\tilde{S}}^{{\rm{z}}}\) should be conserved (alternatively, both Z and Z²).

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 U_I = U_ISOU_D, where \({U}_{{\rm{ISO}}}=\exp (-{\rm{i}}\frac{\pi }{\sqrt{2}}{{\bf{S}}}_{1}\cdot {{\bf{S}}}_{2})\). The gate U_ISO is not a Clifford gate, and conserves S^z, \({U}_{{\rm{I}}SO}^{\dagger }({S}_{1}^{{\rm{z}}}+{S}_{2}^{{\rm{z}}}){U}_{{\rm{ISO}}}={S}_{1}^{{\rm{z}}}+{S}_{2}^{{\rm{z}}}\), but not \({\tilde{S}}_{1}^{{\rm{z}}}+{\tilde{S}}_{2}^{{\rm{z}}}\). The gate U_I therefore conserves only \({S}_{1}^{{\rm{z}}}+{S}_{2}^{{\rm{z}}}\), 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₂ 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.

a S₂(t) for two spin-1 random circuits: a Clifford system with U_D conserving both diagonal operators and resulting in diffusive growth (dotted curves for L = 10, 20, 2000 sites; times are multiplied by 1.4), and a non-Clifford U_I (full curves for L = 10, 14, 18, 20) that conserves only one diagonal operator and is conjectured to lead to the asymptotic linear growth. The logarithm in the definition of S₂ is here base-3. In b we plot a finite-time scaling exponent α (from S₂ ~ t^α) for the same data (solid lines for U_I, dashed for U_D; line-style is changed into a dotted curves for times when S₂ starts to approach a finite-L saturation value). Chain lines suggest convergence with L.

Full size image

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 up-to times t ≈ 6 the two evolutions result in the same growth of S₂ (one could fit S₂ ≈ 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}\asymp \sqrt{t}\) 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₂ ~ 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₂ ~ 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} \sim \sqrt{t}\) only at t ≈ 10³ where S₂ ≈ 10²).

In Fig. 8 we show data for further non-Clifford random circuits. We show results for \({U}_{{\rm{XX}}2}=\exp (-{\rm{i}}\frac{\pi }{2\sqrt{2}}[{S}_{1}^{{\rm{x}}}{S}_{2}^{{\rm{x}}}+{S}_{1}^{{\rm{y}}}{S}_{2}^{{\rm{y}}}])\) that conserves only \({S}_{1}^{{\rm{z}}}+{S}_{2}^{{\rm{z}}}\), but not \({\tilde{S}}_{1}^{{\rm{z}}}+{\tilde{S}}_{2}^{{\rm{z}}}\). Data in frame (a) is compatible with the linear asymptotic growth. If we on the other hand change the gate to \({U}_{{\rm{XX}}}=\exp (-{\rm{i}}\frac{\pi }{\sqrt{2}}[{S}_{1}^{{\rm{x}}}{S}_{2}^{{\rm{x}}}+{S}_{1}^{{\rm{y}}}{S}_{2}^{{\rm{y}}}])\), which conserves both \({S}_{1}^{{\rm{z}}}+{S}_{2}^{{\rm{z}}}\) and \({\tilde{S}}_{1}^{{\rm{z}}}+{\tilde{S}}_{2}^{{\rm{z}}}\), we see that the growth is much slower, like S₂ ~ t^0.7 at short times. While it is hard to make any asymptotic claims about \({S}_{2} \sim \sqrt{t}\) 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 of \({\tilde{S}}_{1}^{{\rm{z}}}+{\tilde{S}}_{2}^{{\rm{z}}}\). Finally, as a third example we show the anisotropic XXZ-like gate \({U}_{{\rm{XXZ}}}=\exp (-{\rm{i}}\frac{\pi }{3}[{S}_{1}^{{\rm{x}}}{S}_{2}^{{\rm{x}}}+{S}_{1}^{{\rm{y}}}{S}_{2}^{{\rm{y}}}+1.5{S}_{1}^{{\rm{z}}}{S}_{2}^{{\rm{z}}}])\) that again conserves only the total magnetization, and therefore one has S₂ ~ t visible already at short times.

Top frame a is for U_XX2, frame b for U_XX, and frame c for U_XXZ (see text for details). The logarithm in the definition of S₂ is always base-3.

Full size image

Data for ladder systems (q = 4)

Here we show further data supporting the claim that for ladders the growth of S₂ 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}_{{\rm{XY}}}=\exp (-{\rm{i}}\frac{\pi }{4}({\sigma }_{j}^{{\rm{x}}}{\sigma }_{k}^{{\rm{x}}}+{\sigma }_{j}^{{\rm{y}}}{\sigma }_{k}^{{\rm{y}}}))\), or \({U}_{{\rm{G}}}={{\rm{C}}NOT}_{12}\ {{\rm{H}}}_{1}\exp ({\rm{i}}\frac{\pi }{4}{\sigma }_{2}^{{\rm{z}}})\). On the rungs we use either \({U}_{{\rm{Z}}Z}={\rm{i}}\exp (-{\rm{i}}\frac{\pi }{2}{\sigma }_{1}^{{\rm{z}}}{\tau }_{2}^{{\rm{z}}})\), or \({U}_{{\rm{S}}xy}=\exp (-{\rm{i}}\frac{\pi }{4}{\sigma }_{1}^{{\rm{z}}})\exp (-{\rm{i}}\frac{\pi }{4}{\tau }_{2}^{{\rm{y}}})\exp (-{\rm{i}}\frac{\pi }{4}{\tau }_{2}^{{\rm{z}}})\). 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}_{{\rm{G}}}^{\dagger }{{\mathbb{1}}}_{j}{\sigma }_{k}^{{\rm{z}}}{U}_{{\rm{G}}}={\sigma }_{j}^{{\rm{x}}}{\sigma }_{k}^{{\rm{z}}}\), and \({U}_{{\rm{G}}}^{\dagger }{\sigma }_{j}^{{\rm{z}}}{{\mathbb{1}}}_{k}{U}_{{\rm{G}}}={\sigma }_{j}^{{\rm{x}}}{{\mathbb{1}}}_{k}\), so that one has \({U}_{{\rm{G}}}^{\dagger }({\sigma }_{j}^{{\rm{z}}}+{\sigma }_{k}^{{\rm{z}}}){U}_{{\rm{G}}}={\sigma }_{j}^{{\rm{x}}}({{\mathbb{1}}}_{k}+{\sigma }_{k}^{{\rm{z}}})\). The rung gate U_ZZ does not break conservation of \({\sigma }_{1}^{{\rm{z}}}+{\sigma }_{2}^{{\rm{z}}}\), nor of \({\tau }_{1}^{{\rm{z}}}+{\tau }_{2}^{{\rm{z}}}\) because one has \({U}_{{\rm{Z}}Z}^{\dagger }{\sigma }_{1}^{{\rm{z}}}{{\mathbb{1}}}_{2}{U}_{{\rm{Z}}Z}={\sigma }_{1}^{{\rm{z}}}{{\mathbb{1}}}_{2}\), and \({U}_{{\rm{Z}}Z}^{\dagger }{{\mathbb{1}}}_{1}{\tau }_{2}^{{\rm{z}}}{U}_{{\rm{Z}}Z}={{\mathbb{1}}}_{1}{\tau }_{2}^{{\rm{z}}}\) (as well as \({U}_{{\rm{Z}}Z}^{\dagger }{\sigma }_{1}^{{\rm{x}}}{{\mathbb{1}}}_{2}{U}_{{\rm{Z}}Z}=-{\sigma }_{1}^{{\rm{x}}}{{\mathbb{1}}}_{2}\), \({U}_{{\rm{Z}}Z}^{\dagger }{{\mathbb{1}}}_{1}{\tau }_{2}^{{\rm{x}}}{U}_{{\rm{Z}}Z}=-{{\mathbb{1}}}_{1}{\tau }_{2}^{{\rm{x}}}\)). 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}_{{\rm{S}}xy}^{\dagger }{\sigma }_{1}^{{\rm{z}}}{{\mathbb{1}}}_{2}{U}_{{\rm{S}}xy}={\sigma }_{1}^{{\rm{z}}}{{\mathbb{1}}}_{2}\), \({U}_{{\rm{S}}xy}^{\dagger }{\sigma }_{1}^{{\rm{x}}}{{\mathbb{1}}}_{2}{U}_{{\rm{S}}xy}= -{\sigma }_{1}^{{\rm{y}}}{{\mathbb{1}}}_{2}\), while it breaks conservation on the lower leg, \({U}_{{\rm{S}}xy}^{\dagger }{{\mathbb{1}}}_{1}{\tau }_{2}^{{\rm{z}}}{U}_{{\rm{S}}xy}={{\mathbb{1}}}_{1}{\tau }_{2}^{{\rm{y}}}\), \({U}_{{\rm{S}}xy}^{\dagger } {{\mathbb{1}}}_{1}{\tau }_{2}^{{\rm{x}}}{U}_{{\rm{S}}xy}={{\mathbb{1}}}_{1}{\tau }_{2}^{{\rm{z}}}\).

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}\asymp \sqrt{t}\). 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₂ ≍ 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₂ ≍ 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₂ (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₂ ≍ t.

Notation (A) − (B) − (C) denotes Clifford evolution with the gate U_A on the upper leg, U_B on the rung, and U_C on the lower leg.

Full size image