Entanglement growth in diffusive systems

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 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 linear in time, similarly as for generic systems without any conservation laws. Exception to this rule are chains of 2-level systems where the dependence is a square-root of time. All predictions are numerically verified by simulations of diffusive Clifford circuits with upto $10^5$ qubits. Such efficiently simulable circuits should be a useful tool for other many-body problems.

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. However, in the thermodynamic limit (TDL) conservation of a single e.g. charge should not matter much in a large Hilbert space. Therefore it was rather surprising when it was presented [13] (focusing on diffusive 1D systems with conserved charge) that the entanglement, as quantified by the Rény entropy S r (t) := log 2 (tr A ρ r A )/(1 − r), grows in fact as S 2 ∼ √ t instead of the "expected" S 2 ∼ t, see also [14,15]. This interesting finding would have many important consequences. For instance, due to simple charge conservation the complexity ∼ 2 S2 would grow with time slower than exponentially (though still superpolynomially)! Computationally a system with diffusive conserved charge would be much less powerful than one without it.
We address the question of the Rény entropy growth in local systems in any spatial dimension d and for any local Hilbert space dimension q. Theoretical predictions, summarized in Table I, 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. Perhaps the most interesting finding is that the entanglement growth is generically in fact linear in time, except arXiv:1912.03645v1 [cond-mat.str-el] 8 Dec 2019 for diffusive qubit (spin-1/2) chains. As a side results, the presented new class of efficiently simulable systems with nontrivial dynamics could be useful in addressing other many-body questions.
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 qdimensional local Hilbert space, whose dynamics has a nontrivial conservation law, e.g. a 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 in particular want to exclude integrable systems that can be special). Linear dimension is denoted by L with the total number of qudits n := L d .
We shall discuss the entropy growth as quantified by the Rény entropy S r (integer 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. We remark that in generic systems all S r , including S 1 , behave in the same way. Whether S 1 grows in the same way also in our case is not completely clear and we leave it to future work. We will mostly focus on S 2 as a representative case of S r =1 , and will use a bipartition to regions A and B. We remark that sometimes S 2 , rather than e.g. S 1 , is the relevant quantity [16], and is also easier to measure [17]. The reduced density operator is ρ A (t) = tr B |ψ(t) ψ(t)|. The size of region A will be extensive, |A| ∼ L d , as otherwise in the TDL one can get strong effects simply due to typicality -tracing a random state over |B| |A| results in the spectrum of ρ A whose relative deviation from a flat one is negligible ∼ q −(|B|−|A|)/2 [18]. To ease comparison of different d and L time will be measured in such units that one will generate a unit of entanglement in a unit of time, S 2 (1) ∼ O(1).
Let us first argue why and how conservation of magnetization (charge) matters for the long-time behavior of S r . As we will see, in the TDL S r is self-averaging and we will for simplicity focus on the purity I(t) := 2 −S2(t) = tr(ρ 2 A (t)). An intuitive meaning of I or S 2 is that is measures the effective number of the explored degrees of freedom needed to describe ρ A (t). In particular, purity is I ∼ 1 N eff , where N eff ∼ 2 L eff is the effective Hilbert space size on which ρ A is supported, while S 2 ∼ log 2 N eff ∼ L eff . Average I over all computational initial states isĪ where D Because of self-averaging, even I(t) for a particular initial state will at large times depend only on dynamics of diagonal operators. Therefore, provided dynamics of all diagonal operators is diffusive, one expects L eff ∼ √ t, resulting in S 2 ∼ √ t. For qubits, q = 2, there are just two local diagonal operators, 1 k and σ z k . If σ z k is diffusive one therefore expects a longtime asymptotic growth S 2 ∼ √ 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 (magnetization). While one might think that diffusive modes that contribute to purity decay as e − √ t will due to their slow decay still dominate over non-diffusive ones which decay as e −t , this is not the case. In a system of n qubits the number of diagonal operators that are products of only diffusive magnetization σ z k and identity 1 k is 2 n , while the number of all other non-diffusive diagonal ones is (q n − 2 n ). The average purity will behave as I ∼ 2 n e − √ t + (q n − 2 n )e −t , so that in the TDL only the exponential term survives. For q ≥ 3 one will have 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 always has linear growth S 2 ∼ t, even at long times. For short-time behavior it is crucial to distinguish correlations spreading in a direction transversal to the boundary of dimension d − 1 between regions A and B, and the tangential direction. 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 nearest-neighbor gates per unit of time. Such scaling is in-line with the mentioned units of time -probability for a random n.-n. gate to connect 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, and A (d−1) is the area of the boundary between A in B (e.g., in 2D this will be a circumference l, in 3D a true two-dimensional area A, in 1D the number of boundaries c). The initial growth of entanglement is therefore expected to be We expect this linear growth to hold for any S r , including r = 1. Linear growth of S 2 will continue until time . After that one will crossover into the asymptotic diffusive growth S 2 ∼ √ t, until at t ∞ ∼ L d+1 a finite-size saturation is reached, 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 infinity to large values of entropy L d−1 < ∼ S 2 < ∼ L d and will be hard if not impossible to experimentally observe. Generic behavior for d > 1 is linear growth, which is in fact as fast as allowed by the Lieb-Robinson bound [19]. Qubits in d = 1 are rather special because the linear growth ends at t 1 ∼ 1 and one gets S 2 ∼ √ t for all not too short times. As we move away from qubits one has generic linear growth for all times. Table I 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 entanglement dynamics are socalled random quantum circuits [20] composed of a series of (random) 2-site unitaries. Random circuits can be thought of as handy toy models of many-body physics but also as a useful theoretical concept being unitary designs [21]. One of the first exact results was obtained by rewriting the dynamics of purity on average as a classical Markov process [22], mapping it to a solvable quantum spin chain and getting an exact expression for the gap ∆ or the decay rate [23], i.e., entanglement speed [24] 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 [23] v E = (1 − 4 5 cos π L ) 1 5 . If one would instead take a regular brick-wall pattern of applied gates, like in [25], one instead has to calculate the gap of a product of Markovian matrices, obtaining a "multiplicative" form v E = 2 ln 5 4 cos π L , going in the TDL to v E 2 ln 5 4 , as also calculated in [25,29]. Studies of random circuits have expanded in recent years with many nice exact results, see e.g. [24][25][26][27][28][29][30]. They are also instrumental in a race for quantum supremacy [31].
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 [32,33] as where ω := e 2πi/q , and all additions are modulo q (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 simple commutation of X and Z, 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 x1 Z z1 ⊗ · · · ⊗ X xn Z zn |a = ω z·x |a ⊕ x .
Evolution of states is however not done by updating each computational basis state -that would be ineffi- cient for highly entangled states -but rather by a stabilizer formalism [34]. A state |ψ on n qubits is called a stabilizer state if it is a joint eigenstate with eigenvalue 1 of n independent stabilizer generators g j from the GPG. Such a state is uniquely determined by n stabilizers g j . For qubits one can obtain it as a product of projectors, |ψ ψ| = Π j (1 + g j )/2, for qutrits one has |ψ ψ| = Π j (1 + g j + g 2 j )/3. A Clifford circuit is a series of unitaries U (called Clifford gates), 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 it would be the case for a generic U ). And here lies the advantage of Clifford circuits. Instead of updating the state |ψ one instead updates each generator g j , whose number is always n and which will remain elements of the GPG [32,34]. Performing one gate, i.e., updating all stabilizers, takes O(n 2 ) operations. Entanglement and state overlaps can also be calculated efficiently [35,36].
A common choice of Clifford gates are: the phase gate P|j = ω j(j−1)/2 |j , the Hadamard gate H|j = 1 √ q k ω kj |k , and a 2-qudit controlled-NOT gate CNOT 12 |j, k = |j, k ⊕ j . The dynamics of Clifford circuits therefore boils down to modular arithmetic [37]. So far such 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 [38] (fluctuations though can exhibit interesting behavior [24]). We shall study a new class of random Clifford circuits that conserve magnetization and whose dynamics is diffusive.
Numerical verification.-We shall first focus on qubits. For qubits the elements of the local GPG are just the ordinary Pauli matrices {σ x , σ y , σ z , 1}. 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 )) on a randomly selected n.-n. pair of sites on a d dimensional square lattice. It is easy to verify that U XY 1 j σ z k U † XY = σ z j 1 k and U 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 | ↑↓ = | ↓↑ . 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 | ↓ . . . ↓↑ . . . ↑ 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).
Starting with an initial state |ψ ∼ (| ↑ + | ↓ ) ⊗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 at late times. It is one of (rare) examples where highly entangling dynamics can nevertheless be simulated. The calculation of entanglement is simplified by the fact that at any time it is composed of an integer number M of generalized EPR pairs, 1 √ q j |j 1 ⊗ |j 2 , 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 entanglement will be generic and consistent with the presented theory. For another solvable evolution that also results in a flat spectrum of ρ A see Ref. [39]. In Fig. 2 we show S 2 for 1D, 2D, and 3D lattice, and for different bipartite splitting of n spins into region A and B. For 1D we see that the asymptotic growth is S 2 = c 2t/π 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, see e.g. [13,14]) it would be very difficult to see the true asymptotic growth. 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. Numerics confirms the short-time growth given by Eq.(2) without any additional prefactors  Fig. 2(a)). In (a) we also show standard deviation (grey shading) and one realization at two selected times.
(middle-1 3 bipartition has A (1) = 4L/3). We also note that due to unfavorable prefactors, to see the asymptotic growth ∼ √ Lt one needs fairly large systems. Even for L = 252 one can see only about one decade in time of S 2 ∼ √ t, 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 ∼ √ L 2 t. Again, for t < t 1 the linear growth (2) 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 middle-1 3 cut, despite a large number of qubits, n = 48 3 = 110592, one can barely hint the eventual ∼ √ t growth. Finally, we show entanglement profiles (Fig. 3), i.e. S 2 for a bipartite cut with region A being the first k spins. Also shown are fluctuations of S 2 between different circuit realizations, showing that the relative fluctuations scale as σ(S 2 )/S 2 ∼ 1/S 2 , and therefore in the TDL at large times dynamics is selfaveraging. It suffices to look at a single random circuit realization.
We also check the case of qudits with q > 2. To that end we simulate a qutrit chain (q = 3, i.e. spin-1 particles) where the local diagonal basis Clifford circuit with nearest-neighbor gate being U = H 2 CNOT 21 CNOT 12 CNOT 12 H 1 , which gives rise to diffusive conservative dynamics of both diagonal matrices, k , we observe the expected diffusive S 2 ∼ √ t (Fig. 4). To get a generic case in which not all diagonal operators are conserved we use local dimension q = 4 visualized as a qubit ladder with the local space corresponding to one rung. Per unit of time we apply L steps, each of which consists 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 magnetizationpreserving U XY on a random bond in the upper leg, or a non-conserving U g = CNOT 12 H 1 Z 2 2 P 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 ∼ √ t. 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.
Conclusion.-We have presented a theory of Rény entanglement entropy growth in lattice systems with conservation of total magnetization due to U (1) symmetry.
We show that in general the entanglement grows linearly in time until a late-time crossover to a slower, parametrically short, square-root growth. In 1D qubit system the diffusive √ t 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 pushed to infinity in the thermodynamic limit. For lattice systems with more than 2 local levels (e.g., spin-1 particles) one will generically observe only linear growth, irrespective of diffusion, the only exception being a case where dynamics of all diagonal operators is diffusive. The case of qubits in 1D displaying a square-root growth is therefore very special: first, for qubits a single U (1) exhausts all non-trivial diagonal operators and second, correlations can propagate only transversely to the boundary accross a single bond, and as such diffusion poses a "bottleneck" already at short times.
Diffusion therefore places only relatively weak influence on entanglement growth, nevertheless, it can distinguish both the spatial dimensionality as well as the size of the local Hilbert space. For other transport types of a conserved quantity and qubits one can conjecture that asymptotically S r ∼ t 1/z , with z being the dynamical exponent. An interesting question is the influence of symmetries other than magnetization conservation. A promising direction is also using presented nontrivial Clifford circuits to further explore the many-body physics.
M.Ž. is supported by Grants No. J1-1698 and No. P1-0402 from the Slovenian Research Agency.