On entropy growth and the hardness of simulating time evolution
Norbert Schuch, Michael M. Wolf, Karl Gerd H. Vollbrecht, J. Ignacio Cirac
aa r X i v : . [ qu a n t - ph ] A p r On entropy growth and the hardnessof simulating time evolution
N. Schuch, M.M. Wolf, K.G.H. Vollbrecht, and J.I. Cirac
Max-Planck-Institut f¨ur Quantenoptik,Hans-Kopfermann-Str. 1, D-85748 Garching, Germany.
Abstract
The simulation of quantum systems is a task for which quantum com-puters are believed to give an exponential speedup as compared toclassical ones. While ground states of one-dimensional systems canbe efficiently approximated using Matrix Product States (MPS), theirtime evolution can encode quantum computations, so that simulatingthe latter should be hard classically. However, one might believe thatfor systems with high enough symmetry, and thus insufficient param-eters to encode a quantum computation, efficient classical simulationis possible. We discuss supporting evidence to the contrary: We pro-vide a rigorous proof of the observation that a time independent localHamiltonian can yield a linear increase of the entropy when actingon a product state in a translational invariant framework. This cri-terion has to be met by any classical simulation method, which inparticular implies that every global approximation of the evolutionrequires exponential resources for any MPS based method.
One of the main challenges in the field of Quantum Computation is to deter-mine the boundaries between classical and quantum computers and pinpointthe computational problems for which quantum devices will be more powerfulthan classical ones. The research on this topic goes in two directions: On theone hand, several efficient quantum algorithms for mathematical problemswhich are believed to be hard classically have been found, the most prominentbeing Shor’s celebrated algorithm for factoring and discrete logarithms [1].On the other hand, the idea to use quantum simulators to simulate quan-tum systems has gained growing attention, as simulating quantum systems isarguably the most natural thing to do with a quantum computer [2]. More-over, in the near future quantum devices are far more likely to be used forsimulating quantum systems than for factoring large numbers, both due to1he much smaller number of qubits required and the weaker demands on thecontrol of the system. In particular, there is no need to have a universalquantum computer available [3–5].In fact, while universal quantum computation requires a dynamics whichis either inhomogeneous in time [6, 7] or in space [8, 9] (or both, as in thestandard network model), physically interesting dynamics to simulate canbe homogeneous in time, i.e., generated by a time-independent Hamiltonian,and translational invariant.Regarding static properties translational invariance appears to consid-erably simplify the problem and enables us in many cases to find efficientclassical approximation algorithms. For instance ground states of local one-dimensional Hamiltonians, in particular in the presence of an energy gap,have an efficient classical description [10, 11] in terms of so-called MatrixProduct States (MPS) [12], which are used as a variational ansatz in the theDensity Matrix Renormalization Group (DMRG) algorithm [13, 14].If static properties of a quantum system are efficiently simulatable on aclassical computer in the presence of translational symmetry, why not alsotime evolution governed by a time-independent local Hamiltonian?The present paper is devoted to providing evidence for the hardness ofclassically simulating quantum mechanical time evolution for systems whichare homogeneous in time and space. Hence, quantum computers (or simula-tors) may indeed yield an exponential speed-up compared to their classicalcounterparts.Clearly, we cannot hope for a complexity theoretic separation betweenclassical an quantum due to the lack of parameters needed to encode a com-putational task. We will rather show first that a sufficiently fast (e.g. linear)increase of the entanglement entropy makes it hard – and for MPS basedapproaches [15, 16] impossible – to even approximate the state of the sys-tem efficiently. That such a linear increase can occur even in the simplestone-dimensional systems was observed in [17] based on arguments of confor-mal field theory together with numerical diagonalization of chains up to 100sites. The main part of the present paper is devoted to a rigorous proof ofthis observation.
We show the hardness of simulating time evolution using MPS in two steps:In a first stage, we take a particular system initially in a translational invari-ant product state of spin- particles, and let it evolve under a translationalinvariant nearest neighbor Hamiltonian. For this system, we prove a linear2ower bound to the entropy of any contiguous block as a function of time,as observed by Calabrese and Cardy in [17]. This evolution has also beenstudied in [18], where similar effects were observed. We then combine thislinear bound with the continuity inequality for the von Neumann entropywhich leads to a linear lower bound on the block entropy of any approxi-mation of the state [19]. As the computational effort needed to deal witha general MPS grows exponentially in the entropy, it would thus increaseexponentially with the time t to be simulated.In the following, logarithms are understood to the base 2. In particular,we define the von Neumann entropy S ( ρ ) = − tr[ ρ log ρ ] with the base 2logarithm. The natural logarithm will be denoted by ln. Theorem 1
Consider a chain with an odd number N of spins with periodicboundary conditions, corresponding to a Hilbert space ( C ) ⊗ N , which at time t = 0 is in the state | ψ (0) i = | ↑ . . . ↑i ≡ | . . . i and evolves under the IsingHamiltonian H = − X j (cid:2) σ xj σ xj +1 + σ zj (cid:3) , (1) so that | ψ ( t ) i = e − i H t | ψ (0) i . Define S L ( t ) = − tr[ ρ L ( t ) log ρ L ( t )] with ρ L ( t ) = tr L +1 ,...,N | ψ ( t ) ih ψ ( t ) | . Then for N ≥ , L ≥ , ≤ t ≤ eL/ ,and t ≤ N/ , S L ( t ) ≥ π t − ln t − . We postpone the proof of the Theorem to Sec. 3.Note that although the lower bound is independent on L – the intuitionbehind being that entanglement arises at the boundaries – the maximallyattainable entanglement is proportional to the length of the block due to theconstraint t ≤ eL/ Let us now see what the entropy scaling implies for an ansatz state | φ i ≡| φ ( t ) i which we use to approximate | ψ i ≡ | ψ ( t ) i . Let ρ L ≡ ρ L ( t ) = tr B | ψ ih ψ | , We choose σ x , σ z the Pauli matrices with eigenvalues ± L ≡ σ L ( t ) = tr B | φ ih φ | be the reduced states on a block A of length L ( B isthe complement of A ). The error we make in this approximation is ǫ ≥ k| ψ ih ψ | − | φ ih φ |k ≥ k ρ L − σ L k =: 2 T , where we have used that the trace norm k η k = tr | η | is contractive underthe partial trace. We now exploit Fannes’ continuity inequality for the vonNeumann entropy in its improved version by Audenaert [20], | S ( ρ L ) − S ( σ L ) | ≤ T log(2 L −
1) + H ( T, − T ) ≤ T L + 1 , where H ( T, − T ) = − T log T − (1 − T ) log(1 − T ) ≤ S ( σ L ) ≥ S ( ρ L ) − T L − ≥ S ( ρ L ) − ǫL − . Using the lower bound on S ( ρ L ( t )) ≡ S L ( t ) of Theorem 1, we have (underthe corresponding conditions on N , L , and t ) that S ( σ L ( t )) ≥ π t − ǫL − ln t − . To get an optimal bound, we choose L as small as possible for the given t , L = 4 t/e , which implies the new constraint t ≥ e/
2, and obain S ( σ L ( t )) ≥ (cid:18) π − ǫe (cid:19) t − ln t − . (2)In order to turn this into a lower bound on the size of an MPS | φ ( t ) i used to approximate the state at time t , note that the entropy of any block isnaturally upper bounded by the rank of the reduced state ρ L , which in turnis at most D . Here, D is the so-called bond dimension of the MPS, whichis polynomially related to the complexity of the MPS [14]. Then, (2) giveslog D ≥ (cid:18) π − ǫe (cid:19) t − ln t − , (3)i.e. the bond dimension – and thus the resources for a straight forward MPSsimulation of the time evolution – will scale exponentially with time as soonas the accuracy ǫ < ǫ = 2 e/ π ≈ . S ( ρ A ( t )) of a region A has a linearleading term, there is an accuracy ǫ beyond which S ( σ A ( t )) has to growlinearly as well. Hence, approaches based on projected entangled pair states(PEPS) [21] or variants thereof are burdened with the same problem.4 Proof of Theorem 1
We now prove Theorem 1. Let us briefly sketch the main steps: First, wemap the spin chain to a one-dimensional system of free fermions, which canbe solved exactly, and bound the error made by going to the thermodymamiclimit N → ∞ . The entropy of a contiguous subblock of length L is equalto the entropy of the corresponding fermionic modes; using a parabola asa bound on the binary entropy we can lower bound the block entropy by aquadratic function in the correlation matrix. Combining this with the exactsolution for the thermodynamic limit, we obtain a lower bound on the blockentropy which is essentially of the form L X n =1 n J n (2 t )(2 t ) , with J n Bessel functions of the first kind. In the appendix we bound thissum by deriving a lower bound for L = ∞ and an upper bound on the errormade by extending the sum. Taking all this together then proves Theorem 1.Before starting with the proof, let us recall the conditions imposed, namely N ≥ L ≥
10, 4 ≤ t ≤ eL/
4, and t ≤ N/
5. In addition we can assume L ≤ N/ We start by mapping the spin operators to fermionic Majorana operators viathe Jordan-Wigner transformation c n = n − Y l =1 σ ( l ) z σ ( n ) x , c n + N = n − Y l =1 σ ( l ) z σ ( n ) y , n = 0 , . . . , N − . (4)We will refer to c , . . . , c N − as the position and to c N , . . . , c N − as themomentum operators; each pair c k , c N + k defines a fermionic mode. For theinitial state | ψ (0) i = | · · · i we have | ψ (0) i = N Y k =1 12 ( ic k c k + N ) = lim β →∞ e + iβ P k c k c k + N β ) , and the Hamiltonian H = i "X j ( c j − c j +1 mod N ) c j + N (5)5s equal to the spin Hamiltonian (1) of Theorem 1 on the − σ (1) z · · · σ ( N ) z , while for the +1 eigenspace, the coupling between spins N and 1 has opposite sign. Since Π | ψ (0) i = −| ψ (0) i for odd N and [ H , Π] = 0,(5) indeed describes the evolution of the Hamiltonian of Theorem 1. Notethat the following results equally hold for even N where they describe anIsing system with a flipped coupling across the boundary.Since the initial state is a Gaussian and the Hamiltonian a quadraticfunction in the Majorana operators, the system remains in a Gaussian statethroughout its evolution and is thus completely characterized by its secondmoments, the antisymmetric correlation matrix [22, 23]Γ kl = − i tr ( ρ [ c k , c l ])which for the initial state is Γ = (cid:18) − (cid:19) . It evolves as Γ t = e − Ht Γ e Ht , (6)where we defined the Hamiltonian matrix H by H = i N − X k,l =0 H kl [ c k , c l ] . We now apply Fourier transformations F kl = e πi kl/N / √ N to both the posi-tion and the momentum subspace. Thereby, Γ and H become block diagonalwith blocks ˆΓ ( φ k ) = (cid:18) −
11 0 (cid:19) (7)and ˆ H ( φ k ) = 12 (cid:18) − e iφ k − e − iφ k (cid:19) , (8)where φ k = 2 π kN . The evolution (6) becomesˆΓ t ( φ k ) = exp h − ˆ H ( φ k ) t i ˆΓ ( φ k ) exp h ˆ H ( φ k ) t i , (9)such that the CM at time t (written in modewise ordering) isΓ t = γ γ · · · γ N − γ − γ . . . ...... . . . . . . ... γ − N · · · · · · γ (10)6ith γ n = 1 N N − X k =0 e inφ k ˆΓ t ( φ k ) . (11)Solving (9), we find that γ n = (cid:18) f n − g n g − n − f n (cid:19) with f n = 1 N N − X k =0 i e inφ k (cid:2) e − iφ k / + e iφ k / (cid:3) sin(2 t sin φ k ) , (12) g n = 1 N N − X k =0 12 e inφ k (cid:2) − e iφ k + (1 + e iφ k ) cos(2 t sin φ k ) (cid:3) d φ . (13)This is the well-known exact solution of the Ising model as found, e.g., in [17,24]. In order to facilitate the evaluation of f n and g n , we consider the limit N →∞ , i.e. we replace 1 N N − X k =0 h ( φ k ) −→ π Z π h ( φ )d φ . We label the functions obtained thereby by f ∞ n and g ∞ n , respectively. Theerror induced by taking the limit can be bounded using the following theorem. Theorem 2 (cf. [25])
Let h : R → R be analytic and π -periodic. Thenthere exists a strip D = R × ( − a, a ) ⊂ C with a > such that h can beextended to a bounded holomorphic and π -periodic function h : D → C . Let M = sup D | h | . Then, ǫ := (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) N N − X j =0 h ( πjN ) − π Z π h ( φ )d φ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ Me aN − . It is straightforward to see that for any a >
0, the addends in both (12) and(13) are holomorphic on D = R × ( − a, a ) and bounded by M = 14 e | n | a (1 + e a )(3 + exp[ t ( e a/ + e − a/ )]) . a = 2 .
9, and taking into account that n ≤ N/ t ≥
4, and N ≥
20, one obtains that | f n − f ∞ n | ≤ ǫ , | g n − g ∞ n | ≤ ǫ ,where ǫ ≤ e − . N +4 . t . (14)Taking the limit in (12) and (13), we find f ∞ n = − π π Z [sin((2 n − α ) + sin((2 n + 1) α )] sin(2 t sin α ) d α = −
12 [ J n − (2 t ) + J n +1 (2 t )]= − nJ n (2 t )2 t (15)and g ∞ n = 12 π Z π [cos(2 nα ) + cos((2 n + 2) α )] cos(2 t sin α ) d α + I n = 12 [ J n (2 t ) + J n +2 (2 t )] + I n = (2 n + 1) J n +1 (2 t )2 t + I n , (16)where I = , I − = − , and I n = 0 otherwise. Here, J k are Bessel func-tions of the first kind, and we have used the recurrence relation [ J n − ( z ) + J n +1 ( z )] = nJ n ( z ) /z . The reduced state of the first L sites is again Gaussian and its correlationmatrix is given by the 2 L × L upper diagonal subblock of Γ t , which we labelby A t . Its entropy can be computed from the normal mode decomposition [26] A t = O L M l =1 (cid:18) λ j − λ j (cid:19) O T , where O ∈ SO(2 L ), as S L ( t ) = L X l =1 h ( λ j ) , with h ( x ) = − x x − − x − x .
8e lower bound h ( x ) ≥ − x as in [23, 27] and thus obtain the bound S L ( t ) ≥ L X l =1 − λ j = L + tr[ A t ] . Writing Γ t = (cid:18) A C − C T B (cid:19) in the 1 , . . . , L and L + 1 , . . . , N partitioning (thus A ≡ A t etc.) and usingthat the purity of the overall state implies Γ t = − A − CC T = − S L ( t ) ≥ tr[ CC T ] = P kl C kl ≡ k C k . We further bound the sum by only considering the lower left and the upperright corner of C , i.e., all entries γ n in (10) with 1 ≤ n ≤ L or N − L ≤ n ≤ N − − L ≤ n ≤ −
1, since γ N − n = γ − n ), and obtain S L ( t ) ≥ L X n = − L | n |k γ n k . (17)We have k γ n k = 2( f n ) + ( g n ) + ( g − n ) ≥ f ∞ n ) + ( g ∞ n ) + ( g ∞− n ) | {z } =: A n − ǫ (2 | f ∞ n | + | g ∞ n | + | g ∞− n | ) | {z } =: B n and with | J n | ≤ / √ | n | ≥ B n ≤ ǫ (cid:18) √ | n | t + 2 | I −| n | | (cid:19) , | n | ≥ . Thus, the total error made in (17) is bounded by12 L X n = − L | n | B n ≤ ǫ √ L ( L + 1)(2 L + 1) t + 1 ! ≤ . L ≤ N/
2, 4 ≤ t ≤ N/
5, and N ≥
20. Thus, S L ( t ) ≥ L X n = − L | n | A n − . ≥ L X k =1 k J k (2 t )(2 t ) − J (2 t )2 t + 14 − . ≥ L X k =1 k J k (2 t )(2 t ) − . . (18)where we have used Eqs. (15) and (16), J − n ( z ) = ( − n J n ( z ), | J (2 t ) | ≤ / √
2, and t ≥
4. We now bound the sum in (18) using Lemma 1 and 2 (seeAppendix ) and obtain S L ( t ) ≥ π (2 t ) − ln(2 t ) − − π π − π − π (2 t ) − L (cid:18) e t L (cid:19) L − . ≥ π t − ln t − ≤ t ≤ eL/ L ≥
10, which proves Theorem 1. (cid:3)
Note that different constraints on t , L , and N can be chosen, resulting indifferent corrections to the leading terms in the expansion of S L ( t ). We gave a rigorous proof to the observation that the entanglement entropyin a simple one-dimensional system increases essentially linearly. Althoughthe proof is tailored to a single exactly solvable model, we think that thisis a widespread property, thus being a requirement which has to be metby every classical simulation method. Of course, there always exist specificcases with tailored representations (e.g. the quasi-free Fermions used here orClifford circuits) close to which one can meet this requirement [28, 29] – evenMPS allow for a description of specific states with only polynomial effortin the entanglement entropy (e.g., if the matrices are tensor products [30]).Yet, these methods seem much more restricted in their applicability thanMPS, so that based on present knowledge it seems unlikely that there isan efficient and generally applicable classical simulation method. Note thatsimilar results can be found by considering any R´enyi entropy with α > global approximation of the state which is a priori notnecessary for obtaining good results for observables with only local non-trivialsupport.
Acknowledgements
We thank B. Horstmann for helpful discussions. This work has been sup-ported by the EU (projects COVAQIAL and SCALA), the cluster of ex-cellence Munich Centre for Advanced Photonics (MAP), the DFG-Forscher-gruppe 635, and the Elite Network of Bavaria project QCCC.
Appendix: Bounds on Bessel sums P n J n ( z ) Lemma 1
For z ≥ , Q ( z ) := ∞ X n =1 n J n ( z ) ≥ π z − z ln z − − π π z − π − π . (19) Proof.—
Differentiating Q using the recurrence relations2 J ′ n ( z ) = J n − ( z ) − J n +1 ( z ) J n ( z ) = z n [ J n − ( z ) + J n +1 ( z )] (20)yields Q ′ ( z ) = z ∞ X n =1 n [ J n − ( z ) − J n +1 ( z )]= z " J ( z ) + ∞ X n =1 nJ n ( z ) | {z } =: f ( z ) . (21)To obtain a lower bound on f ( z ), we differentiate again and find f ′ ( z ) = 2 z ∞ X n =1 [ J n − ( z ) − J n +1 ( z )]= 2 z [ J ( z ) + J ( z )] . (22)As we prove in Lemma 3 later on, J ( z ) + J ( z ) ≥ πz − z , z ≥ . J ( z ) + J ( z ) ≥ ≤ z ≤
1, we obtain by integration f ( z ) ≥ z − π − z for z ≥ f ( z ) ≥ J ( z ) ≥
0, we have Q ′ ( z ) ≥ zf ( z ) / z ≥ (cid:3) Note that the proof can be extended straightforwardly to sums P ∞ n = K . Lemma 2
For z ≤ eK/ , K ≥ , ˜ Q ( z ) = ∞ X n = K +1 n J n ( z ) ≤ Kz (cid:16) ez K (cid:17) K . (23) Proof.—
As in the proof on Lemma 1, we differentiate ˜ Q using the recurrencerelations (20) and find˜ Q ′ ( z ) = z h ( K + 1) J K ( z ) + ( K + 2) J K +1 ( z ) + ∞ X n = K +2 nJ n ( z ) | {z } =: ˜ f ( z ) i and ˜ f ′ ( z ) = 2 z (cid:2) J K +1 ( z ) + J K +2 ( z ) (cid:3) . We now use the lower bound | J n ( z ) | ≤ ( z/ n /n !, which after integrationleads to ˜ Q ( z ) ≤ C ( K, z ) (cid:16) z (cid:17) K +2 with a prefactor C ( K, z ) of order (1 /K !) , which for z ≤ eK/ K ≥
2, andusing k ! ≥ R ∞ k e − t t k d t ≥ ( k/e ) k , straighforwardly reduces to (23). (cid:3) Lemma 3
For z ≥ , J ( z ) + J ( z ) ≥ πz − z . Similar bounds can be proven for any pair J n + J n +1 along the same lines.Note that the asymptotic scaling J n ( z ) + J n +1 ( z ) ∼ /πz is well known [31]. Proof.—
We use the expansion ( [31], Sec. 7.3) J n ( z ) = r πz " cos( z − nπ − π p − X m =0 ( − m ( n, m )(2 z ) m + P ( z ) ! sin( z − nπ − π q − X m =0 ( − m ( n, m + 1)(2 z ) m +1 + Q ( z ) ! with ( n, m ) := Γ( n + m + 1 / m !Γ( n − m + 1 / , where the remainders P and Q lie between zero and the first neglected termof the respective series, given that z >
0, 2 p > n − , and 2 q ≥ n − . With J n ( z ) = S n ( z ) + E n ( z ), where S n ( z ) comes from the series and E n ( z ) fromthe remainders, we have | J n ( z ) | ≥ | S n ( z ) | − | E n ( z ) | and thus J n ( z ) ≥ S n ( z ) − | S n ( z ) || E n ( z ) | . (24)We choose p = q = 1, and first upper bound | S n ( z ) | ≤ √ πz ( n = 0 ,
1) for z ≥ and then, using (24), J ( z ) + J ( z ) ≥ πz − πz − √ πz + 132 πz − √ πz ≥ πz − z for z ≥
1, using | sin( z ) | ≤ | cos( z ) | ≤ (cid:3) References [1] P. W. Shor, SIAM J. Comput. , 1484 (1997), quant-ph/9508027.[2] S. Lloyd, Science , 1073 (1996).[3] A. Sorensen and K. Molmer, Phys. Rev. Lett. , 2274 (1999),quant-ph/9903044.[4] E. Jane, G. Vidal, W. D¨ur, P. Zoller, and J. I. Cirac, Quant. Inf.Comput. , 15 (2003), quant-ph/0207011.[5] D. Porras and J. I. Cirac, Phys. Rev. Lett. , 207901 (2004),quant-ph/0401102.[6] K. G. H. Vollbrecht and J. I. Cirac, Phys. Rev. A , 012324 (2005),quant-ph/0502143.[7] R. Raussendorf, Phys. Rev. A , 052301 (2005), quant-ph/0505122.[8] K. G. Vollbrecht and J. I. Cirac, Phys. Rev. Lett. , 010501 (2007),arXiv:0704.3432. 139] A. Kay, Phys. Rev. A , 030307(R) (2007), arXiv:0704.3142.[10] M. B. Hastings, Phys. Rev. B , 035114 (2007), cond-mat/0701055.[11] M. Hastings, J. Stat. Mech. P08024 (2007), arXiv:0705.2024.[12] A. Kl¨umper, A. Schadschneider, and J. Zittartz, J. Phys. A , L955(1991); Z. Phys. B , 281 (1992).[13] S. R. White, Phys. Rev. Lett. , 2863 (1992).[14] U. Schollw¨ock, Rev. Mod. Phys. , 259 (2005), cond-mat/0409292.[15] G. Vidal, Phys. Rev. Lett. , 040502 (2004), quant-ph/0310089.[16] T. J. Osborne, Phys. Rev. Lett. , 157202 (2006), quant-ph/0508031.[17] P. Calabrese and J. Cardy, J. Stat. Mech. P04010 (2005),cond-mat/0503393.[18] M. Cramer, C.M. Dawson, J. Eisert, and T.J. Osborne, Phys. Rev. Lett. , 030602 (2008).[19] N. Schuch, M. M. Wolf, F. Verstraete, and J. I. Cirac, (2007), Phys.Rev. Lett. , 030504 (2008), arXiv:0705.0292.[20] K. M. R. Audenaert, J. Phys. A , 8127 (2007), quant-ph/0610146.[21] F. Verstraete and J. Cirac, (2004), cond-mat/0407066.[22] S. Bravyi, Quant. Inf. Comput. , 216 (2005), quant-ph/0404180.[23] M. M. Wolf, Phys. Rev. Lett. , 010404 (2006), quant-ph/0503219.[24] S. Sachdev, Quantum Phase Transitions (Cambridge University Press,1999).[25] R. Kress,