Power-law entanglement growth from typical product states

Generic quantum many-body systems typically show a linear growth of the entanglement entropy after a quench from a product state. While entanglement is a property of the wave function, it is generated by the unitary time evolution operator and is therefore reflected in its increasing complexity as quantified by the operator entanglement entropy. Using numerical simulations of a static and a periodically driven quantum spin chain, we show that there is a robust correspondence between the entanglement entropy growth of typical product states with the operator entanglement entropy of the unitary evolution operator, while special product states, e.g. $\sigma_z$ basis states, can exhibit faster entanglement production. In the presence of a disordered magnetic field in our spin chains, we show that both the wave function and operator entanglement entropies exhibit a power-law growth with the same disorder-dependent exponent, and clarify the apparent discrepancy in previous results. These systems, in the absence of conserved densities, provide further evidence for slow information spreading on the ergodic side of the many-body localization transition.


I. INTRODUCTION
Entanglement is an intrinsic property of quantum many-body wave functions and describes the amount of information a subsystem A contains about its complement B. The entanglement between a bipartition A :: B of a closed quantum system is quantified by the entanglement entropy, which distinguishes separable product states with zero entanglement from states which are entangled to various degrees. For example, it was found that ground states of gapped quantum many-body systems exhibit an entanglement entropy which scales as the surface area of the subsystem A [1][2][3]. In generic quantum many-body systems that thermalize and in which the eigenstate thermalization hypothesis [4][5][6][7][8][9] is valid, the entanglement entropy of all eigenstates of the Hamiltonian is furthermore identified with the thermodynamic entropy of the subsystem and is therefore proportional to the volume of the subsystem for eigenstates corresponding to finite temperatures [10][11][12][13][14][15].
The phenomenology becomes much richer in a nonequilibrium setting after a quantum quench: even if the wave function at early times is prepared as a lowly entangled state, the complex quantum many-body dynamics leads to a rapid production of entanglement and typically the entanglement entropy is found to grow linearly in time under the evolution with short range Hamiltonians [16][17][18][19]. This was also observed in experiments with ultracold atomic systems [20].
The production of entanglement is typically observed in the quench dynamics of wave functions |ψ(t) . However, it is uniquely due to the action of the time evolution operator U (t) on the wave function and should therefore be reflected in the growing complexity of U (t). Consequently, a generalization of the entanglement entropy to operators was introduced in [40,41], which allows for a state-independent quantification of the complexity of operators across a bipartition of the system. In Refs. [42][43][44], this concept was directly applied to the time evolution operator U (t), and a generic linear growth of the operator entanglement entropy (opEE) in ergodic quantum systems was found, while MBL systems are characterized by a logarithmic growth in time [42].
In the vicinity of the MBL transition, the situation is more complicated: At intermediate disorder, where the system still thermalizes, an anomalously slow thermalization [14,[45][46][47][48][49] is found, which was connected to slow, subdiffusive transport [50][51][52][53][54][55][56][57][58][59][60][61][62][63][64]. The anomalous thermalization is also reflected in a sublinear powerlaw growth ∝ t α of the wave function entanglement entropy [65], even in systems which do not have globally conserved densities [66], suggesting that the generic slow dynamics is a universal precursor of MBL. In such pre-MBL systems, the entanglement production exponent α varies continuously with disorder and vanishes at the MBL transition, where the logarithmic growth takes over. The same phenomenology was found for the growth of the operator entanglement entropy of the time evolution operator U (t), with a disorder-dependent exponent. Comparing the exponents found in the disordered Heisenberg model [35,51,[67][68][69][70], of the entanglement entropy growth after a quench from a product state [65] and of the operator entanglement entropy growth of U (t) [42], it turns out that these exponents do not agree.
In the present article, we address the above disagreement and show that it is due to the fact that completely unentangled σ z product states are not typical separable states, and that the growth of the operator entanglement entropy agrees perfectly with the growth of the wave function entanglement entropy if typical product states are considered. We provide further evidence that the phenomenon of generally slow dynamics is universally found in one-dimensional systems at disorder strengths weak enough so that the system still thermalizes and does not require any conserved densities with associated transport.
Our article is organized as follows: In Sec. II we introduce a static and a periodically driven model representing generic disordered quantum XYZ spin chains, as well as the methods used to characterize them. We show that both models exhibit an MBL transition. The static model possesses no symmetries nor extensive conservation laws in addition to energy, while also the energy conservation is broken in the driven model. In Sec. III we present the framework of the entanglement measures considered for both wave functions and operators. In Sec. IV we compare the growth of the wave function entanglement entropy after a global quench, S ψ (t), with the growth of the operator entanglement entropy of evolution operators, S U (t). For both models, we find the same universal power-law growth t α of S ψ (t) and S U (t) within the ergodic regime-provided that S ψ (t) is obtained after a quench from typical initial product states. If the quench comes from either σ z product states or a class of intermediate states, S ψ (t) and S U (t) differ. Such a class of intermediate states is characterized by their maximal bond dimension in a matrix product state representation and a precise definition is given in Sec. IV C. We further elucidate the influence of the initial states on the wave function entanglement production and argue that the underlying mechanism is that of monogamy of entanglement. Finally, in Sec. V we conclude by summarizing and discussing our main results.

A. Static and Driven XYZ chain
We study the generic XYZ chain with open boundary conditions in the presence of a disordered tilted field. Its Hamiltonian is given by where σ γ i denote Pauli matrices, J γ coupling constants and h γ ≡h γ /| h| fixed amplitudes of the tilted field h; γ ≡ x, y, z. The field is fixed in the direction of h but its magnitude is disordered and taken from a box distribution h i ∈ [−W, W ]. We set (J x, J y , J z ) = (0.5, 0.7, 1.0) and (h x ,h y ,h z ) = (0.95, 1.0, 1.1) to remove all the extensive conservation laws except energy conservation. This means that we can not further reduce the manybody Hilbert space and there is only one block of size dim(H) = 2 L .
To also break the energy conservation, we subject the system (1) to monochromatic driving. The resulting Floquet model is defined by with driving period T ≡ 2π/ω = 2 and driving amplitude A = 0.5. The Floquet operator over one driving period is defined aŝ where T denotes the time-ordering operator. Sincê U F (T ) is a unitary operator, its eigenvalues ω n lie on the complex unit circle.

B. Method
For the characterization of the models, we fully diagonalize the Hamiltonian of the static model (1) up to system size L = 14 with dim(H) = 16384, and consider the statistics of adjacent energy spacings. In the case of the Floquet model (2), we have to generate the timeevolution operatorÛ F (T ). Since our driving protocol is monochromatic, we use small time steps dt = 0.02 and a second-order Trotter decomposition to separate the constant part exp(−idtH 0 ) from the (diagonal) driven part exp(−idtH D (t)), requiring only the diagonalization of H 0 to calculate the corresponding matrix exponential and repeated matrix products to step through the period T . We then fully diagonalizeÛ F (T ) and consider the statistics of adjacent eigenphases to establish the ergodic regime of the model.
In the main part of this work, we consider both the production of wave function entanglement when starting from a product state |ψ(t = 0) as well as the growth of operator entanglement directly in the time evolution operator U .
For the wave function dynamics, we use exact-time evolution with a Krylov-space method [60,65,[71][72][73], which relies on the sparse structure of both H 0 and H D . Again, to faithfully describe the monochromatic driving, small timesteps dt = 0.02 inside the period are used and the results analyzed including intraperiod values. This method allows for the calculation of entanglement entropies in systems up to L = 26 with dim(H) = 2 26 ≈ 6.7 · 10 7 .
For the operator entanglement entropy, we calculate the time evolution operator U (t) as described above with the same limitation to system sizes up to L = 14.
Our results are averaged over 50-100 disorder realizations, for several values of disorder strength within the ergodic regime at intermediate disorder W

C. Characterization of the Model
The static model is similar to other disordered spin chains, in particular, the Heisenberg spin chain given by J x = J y = J z = 1, h x = h y = 0 and h z = 1, which has in these units an MBL transition at W c ≈ 7.5 ± 2.0 [51,[67][68][69][70], therefore we expect to find an MBL transition at similar values of disorder.
Using level spacing statistics as diagnosis, we corroborate that for the chosen parameters both models (1) and (2) undergo an MBL transition at critical values of disorder W 0 c ≈ 6 and W D c ≈ 12, respectively. The adjacent level spacing ratio is defined as [74] where δ n is defined in terms of consecutive energy levels δ n ≡ E n+1 − E n of (1) or consecutive phase spacings δ n ≡ θ n+1 −θ n extracted from the eigenvalues ω n = e −iθn of (3), depending on whether the static or the Floquet model is considered. For the static model (1) the energy is conserved and different energy densities represent different temperatures of the system. Since in the Heisenberg model numerical evidence for a mobility edge (i.e. an energy density-dependent critical disorder W c ) was found, we study the average level spacing ratio [r] for different energy densities = (E n − E min )/(E max − E min ) defined as in Ref. [69]. We include the eigenvalues E n of the Hamiltonian on a vicinity of radius δE = 0.1 around = 0.5 and = 0.25 and average over 100-10000 realizations of the disorder. The results of this analysis for different system sizes L are shown as a function of disorder strength W in Fig. 1(a),(b), revealing a transition from a mean level spacing expected from the Gaussian unitary ensemble (GUE) at weak disorder with [r] GUE ≈ 0.59982(8) [75] to the value predicted by a Poisson distribution [r] POI ≈ 0.38629 [76] at strong disorder. Results for different system sizes show a crossing at intermediate disorder which, although the precise location still drifts with system size seems compatible with the existence of a mobility edge, as the crossing appears at significantly lower disorder for = 0.25 compared to = 0.5.
In contrast, Floquet models which have an MBL transition do not exhibit mobility edges [29,46]. Therefore, we include eigenvalues from the full unit circle to calculate the level spacing statistics. The result is shown in Fig. 1(c), where the crossing of [r] with L is now enclosed between the limit values corresponding to the Poisson and the circular unitary ensemble (CUE) distributions, [r] POI ≈ 0.38629 and [r] CUE ≈ 0.59982(8) [77], respectively. As expected [29], the critical disorder of the Floquet model is larger than in the static case. In the remainder of this paper, we will focus on the region of weak enough disorder such that both models are well located in the ergodic phase (W ≤ 4).

III. WAVE FUNCTION AND OPERATOR ENTANGLEMENT ENTROPY
In this section, we provide a technical discussion pertinent to the concepts of the entanglement entropy of both quantum wave functions |ψ and quantum evolution op-eratorsÛ for the case of a complementary real-space bipartition in terms of two subsystems A and B ≡ A, with a tensor product Hilbert space H = H A ⊗H B . For a more detailed discussion we refer the reader to [40,42,43].
The system in a pure state |ψ ∈ H has a density matrix ρ = |ψ ψ| with unit purity Tr ρ 2 = 1. The state of the subsystem A is described by the reduced density matrix ρ A , given by the partial trace of ρ over the subsystem B, Even if ρ is a pure state, the reduced density matrix is in general a mixed state with purity Tr ρ 2 A ≤ 1. The von Neumann entropy of the reduced density matrix ρ A is called the entanglement entropy and defined as Any pure state |ψ ∈ H in turn can be expressed in terms of its Schmidt decomposition where form a complete orthonormal basis of H A (H B ) and are obtained from the left and right singular vectors of the matrix Ψ. Algorithmically, (when using a complete computational basis ordered such that the tensor product structure is preserved), the entanglement spectrum { λ i } of a wave function |ψ is obtained by reshaping the wave function into a matrix Ψ such that the row indices correspond to basis states of subsystem A and the column indices correspond to basis states of subsystem B. Then, a singular value decomposition (SVD) of Ψ yields the singular values { √ λ i }. The von Neumann entanglement entropy (6) is readily generalized to the α Rényi entropy in terms of the squared singular values (the entanglement spectrum) {λ i }; from which (6) is recovered in the limit α → 1.
From the foregoing concepts, the generalization of the entanglement entropy to the space of linear operators is straightforward. Linear operatorsÔ : H → H form a basis of the Hilbert spaceH : H → H, endowed with the inner product ·, · :H ×H → C, which in turn is inherited from H. Given two linear operatorsÔ,Ô ∈H, their inner product is defined as The same extension can be done for the aforementioned bipartition (A :: B) in terms of the Hilbert spaces associated to each subsystem,H A : H A → H A andH B : H B → H B , respectively. Given any linear operator-in particular-a unitary quantum evolution operator which obeys the orthonormality (unitarity) condition Û ,Û = 1, the above definitions (5)-(8) can be extended using Eq.(9). Again, using a computational basis order which respects the tensor-product structure of the Hilbert space, we can write any basis state |i = |i A , i B = |i A ⊗ |i B in terms of basis states of the subsystems A and B. Then, the time evolution operatorÛ has a matrix representation in the form Similarly to the case of wave functions, but now dealing with two pairs of indices, we can calculate the operator entanglement spectrum by first vectorizing the matrix U (interpreting it as a vector inH ⊗H) and then writing it as a matrix u with all the indices corresponding to the A subsystem as row indices and the indices corresponding to the B subsystem as column indices: Algorithmically, this corresponds to interpreting the unitary matrix U as a tensor of rank 4, then performing a tensor transposition to sort the A and B indices, followed by a reinterpretation as a (rectangular) matrix u. The operator entanglement spectrum { λ op i } is readily obtained by an SVD of u.
In the following, we will concern ourselves with comparing the entanglement dynamics of wave functions and time evolution operatorsÛ in the static (1) and the Floquet model (2), taking into account the following consideration. Since there is a Hilbert space isomorphism between the space of states and the space of linear operators H⊗H ∼ =H, the comparison between S ψ (t) and S U (t) should be done with respect to system sizes L and L/2, respectively, which correspond to the same Hilbert space dimension. The Hilbert space dimension determines the maximal entanglement entropy in a finite system, which is given by S ψ max = ln min [dim(H A ), dim(H B )] for wave functions and S U max = ln min dim(H A ), dim(H B ) = 2S ψ max for operators, where the dynamics in ergodic systems is expected to reach values close to the maximum at late times [78].
While in the case of the time evolution operator, there is no ambiguity with respect to the initial state (Û (t = 0) = 1 and therefore S U (t = 0) = 0), we are free to choose any wave function |ψ(t = 0) as an initial state for the time evolution. Since we are interested in the production of entanglement, it is natural to require that the initial state be a product state |ψ(t = 0) = |ψ A ⊗ |ψ B which has minimal entanglement S ψ (t = 0) = 0.
In this work, we mostly focus on typical product states, which are defined by random Haar measure states |ψ A and |ψ B on each subsystem. These states are the most general product states and have zero entanglement between A and B, while being maximally entangled inside each subsystem. We will discuss the dependence on initial states for different classes of product states in Sec. IV C, including the much simpler (and less typical) σ z basis states which are widely used in related numerical studies.

A. Static model
In Fig. 2 we show the growth of the entanglement entropy of |ψ(t) (solid lines) following a global quench from a typical product state of the form (12), and compare it with the growth of the operator entanglement entropy of U (t) (dashed lines) for disorder strengths W = 1, 2, 3, 4, all well in the ergodic regime, as discussed in Sec. II C. All results are averaged over ≈ 100 realizations of the disorder and for each disorder realization a different initial product state is used.
In the upper panels of Fig. 2 we observe that for all values of disorder (starting from zero entanglement) there is a rapid growth of both entropies at short times and a saturation at late times. The saturation values are close to the Page value (indicated by dotted-dashed lines), as expected [42]. Interestingly, the growth of the operator entanglement entropy is almost identical to the wave function entanglement entropy growth, clearly showing that they encode the same information.
The doubly logarithmic scale reveals that the growth of both entropies follows a power law t α in time until saturation, and the domain of the power law extends to later times for larger systems due to the larger saturation value (proportional to L).
To analyze the value of the dynamical exponent α, we show the discretized logarithmic derivative d ln S(t)/d ln t of both entropies as a function of time in the lower panels of Fig. 2, where the derivative is taken over small time windows of size δt = 0.1 (results are essentially independent of the choice of δt). This analysis clearly shows that both the operator entanglement entropy of U (t) and the wave function entanglement entropy of |ψ(t) grow like a power law with the same exponent α. The domain of the power law grows with system size, stabilizing the plateau in the logarithmic derivative. And, most importantly, the value of the exponent decreases continuously as a function of disorder, confirming that there is slow dynamics in the system before it undergoes an MBL transition. These results appear to be converged with system size at short enough times, and due to the slower dynamics the domain of the power law is larger for stronger disorder W .

B. Floquet model
Analogously to the case of the static model, in Fig. 3 we show the comparison between the growth of the wave function entanglement entropy of typical product states |ψ(t) and the operator entanglement entropy of U (t), generated by the monochromatic drive (2).
The results are very similar to the static case: We find a power law growth of both entropies at short times until it saturates to a value close to the Page value at late times. The domain of the power law grows with system size L and the operator entanglement entropy of U (t) follows very closely the wave function entanglement entropy (comparing the entanglement entropy of a system of size L to the opEE of a system of size L/2).
We find again that at stronger disorder (still well in the ergodic regime), the exponent α of this power law is strongly suppressed, indicating slow dynamics in this system prior to the MBL transition.

C. Initial state dependence
In Figs. 2 and 3, we have shown results for the entanglement entropy production starting from a typical product state, which is maximally entangled inside the subsystems A and B.
It is interesting to ask the question whether other classes of product states (which are less strongly entangled inside the subsystems) yield the same results. We therefore introduce the following general matrix product state (MPS) ansatz for any product state with respect to the A :: B bipartition: Here, A is the length of the subsystem A and M σ k ∈ C χ×χ are independent random matrices with i.i.d Gaussian elements associated to each site and spin polarization of the system. At the edges k = 1, L of the system and at the boundary k = A , A + 1, M σ k ∈ C χ are vectors instead, making the wave function a product state if cut at the boundary between A and B. The wave function |ψ χ is normalized by a proper choice of the constant N .
The typical product states used in Figs. 2 and 3 correspond to the maximal bond dimension χ = 2 A /2 (for The legends in (c) apply to all panels and legends in (d) apply to the lower panel.
If we choose instead χ = 1, we recover product states with random phases on each site but zero entanglement, independently of the bipartition of the system. The typical product states were also considered in [80], and the special case of χ = 1 with and without random phases in Refs. [18,19,65,66].
The general ansatz in Eq. (13) allows for a smooth interpolation between these extreme cases by choosing different intermediate bond dimensions χ.
In Fig. 4, we present the growth of the entanglement entropy in both the static and Floquet models for W=2.0 for different classes of initial product states defined by Eq. 13, labelled by their bond dimensions χ. In addition to χ = 1, we show data for the χ = 1 case without random phases, i.e. for pure σ z basis states.
Somewhat surprisingly, different classes of product states show a remarkably different behavior in terms of their entanglement growth even though we consider only the ergodic phase of our models. We note in passing that different (logarithmic) entanglement production rates were found in the MBL phase when considering different types of χ = 1 product states [81].
The entanglement in completely unentangled product states grows the fastest, while it grows more slowly in typical product states which are maximally entangled within the subsystems. This different growth rate seems to be reflected in different exponents α of the power law in time.
We argue here that typical product states with maxi-mal χ are the most representative for the overall behavior of the system in the sense that they contain the largest number of degrees of freedom. This means that if one generates a random product state, the likelihood of finding a χ = 1 state vanishes compared to a maximal χ state.
Why is the growth of entanglement slower if we start from a product state which is initially already entangled inside each subsystem, compared to an unentangled χ = 1 product state? The different behavior in the two cases is illustrated in Fig. 5, where the top panels exhibit the entanglement entropy growth as a function of time for different sizes of the subsystem. Panel (a) shows the case of a typical product state, which has zero entanglement for A = L/2 by construction and maximal entanglement for all other values of A given this constraint. Panel (c) shows the same data as line plots for different A . It is clear that the entanglement entropy for a bipartition which is already close to maximally entangled can only grow very little. However, due to the monogamy of entanglement [82,83], in order to generate entanglement across the cut at A = L/2, highly nontrivial processes have to occur to "free" some degrees of freedom before they can entangle with the other subsystem. This is reflected in the flat behavior (no growth) of cuts at e.g. A = 5 at short times and slows down the entanglement growth across the cut at the centre.
Conversely, the case of a χ = 1 σ z product state does not impose such constraints and the entanglement for all cuts builds up homogeneously.
The entanglement production at very early times t 1 is different from the situation described above, the typical product states produce entanglement much faster, compared to the χ = 1 product states, and it is after this fast start that the entanglement production becomes slower. Interestingly, this t 1 behavior was observed for similar typical product states [80], where the faster entanglement production was attributed to the positive curvature ∂ 2 S( A , t)/∂ 2 A at the central cut A = L/2, compared to the flat structure produced by the χ = 1 initial condition. This situation is reflected in Fig. 5 (a),(b) at t 1 around the central cut. The aforementioned second partial derivative is the subleading correction to the coarsegrained behavior of the local entanglement entropy increase rate and it is explained by the phenomenology of entanglement production [80,84,85], which is conjectured to apply to generic non-integrable systems.

V. DISCUSSION
We have systematically compared the operator entanglement entropy of the time evolution operator U of disordered static and driven quantum spin chains, well in the ergodic regime of the phase diagram. Our models are chosen such that they exhibit a many-body localization transition at strong disorder.
In the case of the disordered quantum Heisenberg chain, the operator entanglement entropy of the time evolution operator was calculated in Ref. [42], and a slow, sublinear power-law growth was found. However, comparing these exponents to the exponents of the power-law growth of the wave function entanglement entropy after a quench from a χ = 1 σ z product state revealed a discrepancy of these exponents.
Here, we clarify this discrepancy: Our argument given above based on the concept of monogamy of entanglement explains why χ = 1 product states can exhibit genuinely faster entanglement production compared to typical product states. However, since the time evolution operator needs to encode the entanglement production for any initial state, its entanglement entropy can only be expected to grow with a rate similar to that seen in typcial wave functions, which are given by initial product states with maximal bond dimension χ and represent an overwhelming majority in the class of all A :: B separable pure states.
This becomes even more apparent if we consider the time evolution operator U in the computational σ z basis: Take an σ z product state |φ , which is given by a single basis state |j : i|φ = φ i = δ i,j . Then, the time evolution of this state is given by U ik φ k = U ij , which is the j-th column of U . On the other hand, when considering a product state with maximal χ, we will get a random average over all columns of U , which will therefore reflect the typical behavior of all columns.
At times t 1 the opposite situation happens, the typical product states exhibit faster entanglement production compared to the χ = 1 product states.
In summary, we have clarified the apparent discrepancy between the growth of entanglement entropies of wave functions and time evolution operators in the slow dynamical regime prior to the MBL transition, showing that typical product states display identical power law entanglement entropy growth exponents compared with the operator entanglement entropy of the evolution operator. This is in agreement with the same correspondence recently observed in the case of linear growth in [80]. Furthermore, our results show that this correspondence is valid across the full ergodic phase and holds both for a static model with a mobility edge as well as for a periodically driven system.
Our results furthermore provide additional evidence that the slow dynamics in the ergodic phase is visible even in general situations of systems without conservation laws and can be observed in power-law growths of quantum information measures such as the operator entanglement entropy of the evolution operator and the entanglement production in wave functions.

VI. ACKNOWLEDGEMENTS
We thank David A. Huse for useful comments on the manuscript and in particular for pointing us to the early time behavior.
Appendix A: Rényi operator entanglement entropies In this appendix we provide further results on Rényi operator entanglement entropies. Rényi entropies are a generalization of the von Neumann-Shannon entropy. They highlight the behavior of different scales in the entanglement spectrum via the Rényi index α. Here, we show results for the growth of the α Rényi operator entanglement entropy of the time evolution operator U (t) of the static model in Eq. (1). In Fig. 6, the top panel shows that the Rényi entropies grow for all values of α. Large α highlights the behavior of the largest singular value, while smaller α > 0 represent their average behavior. The lower panel shows the discretized logarithmic derivative, which reveals that all Rényi entropies with 0.4 ≤ α ≤ 2 essentially have the same power law growth exponent in time. For very large Rényi index α = 100, the behavior of the largest singular value is recovered and the entropy is identical to the limit α → ∞. Interestingly, in this case we observe a slightly slower entropy growth with a smaller exponent.