On entropy growth and the hardness of simulating time evolution

The simulation of quantum systems is a task for which quantum computers are believed to give an exponential speedup as compared to classical ones. While ground states of one-dimensional systems can be efficiently approximated using Matrix Product States (MPS), their time evolution can encode quantum computations, so that simulating the latter should be hard classically. However, one might believe that for systems with high enough symmetry, and thus insufficient parameters to encode a quantum computation, efficient classical simulation is possible. We discuss supporting evidence to the contrary: We provide a rigorous proof of the observation that a time independent local Hamiltonian can yield a linear increase of the entropy when acting on a product state in a translational invariant framework. This criterion has to be met by any classical simulation method, which in particular implies that every global approximation of the evolution requires exponential resources for any MPS based method.


Introduction
One of the main challenges in the field of Quantum Computation is to determine the boundaries between classical and quantum computers and pinpoint the computational problems for which quantum devices will be more powerful than classical ones. The research on this topic goes in two directions: On the one hand, several efficient quantum algorithms for mathematical problems which are believed to be hard classically have been found, the most prominent being Shor's celebrated algorithm for factoring and discrete logarithms [1]. On the other hand, the idea to use quantum simulators to simulate quantum systems has gained growing attention, as simulating quantum systems is arguably the most natural thing to do with a quantum computer [2]. Moreover, in the near future quantum devices are far more likely to be used for simulating quantum systems than for factoring large numbers, both due to the much smaller number of qubits required and the weaker demands on the control of the system. In particular, there is no need to have a universal quantum computer available [3][4][5].
In fact, while universal quantum computation requires a dynamics which is either inhomogeneous in time [6,7] or in space [8,9] (or both, as in the standard network model), physically interesting dynamics to simulate can be homogeneous in time, i.e., generated by a time-independent Hamiltonian, and translational invariant.
Regarding static properties translational invariance appears to considerably simplify the problem and enables us in many cases to find efficient classical approximation algorithms. For instance ground states of local onedimensional Hamiltonians, in particular in the presence of an energy gap, have an efficient classical description [10,11] in terms of so-called Matrix Product States (MPS) [12], which are used as a variational ansatz in the the Density Matrix Renormalization Group (DMRG) algorithm [13,14].
If static properties of a quantum system are efficiently simulatable on a classical computer in the presence of translational symmetry, why not also time evolution governed by a time-independent local Hamiltonian?
The present paper is devoted to providing evidence for the hardness of classically simulating quantum mechanical time evolution for systems which are homogeneous in time and space. Hence, quantum computers (or simulators) may indeed yield an exponential speed-up compared to their classical counterparts.
Clearly, we cannot hope for a complexity theoretic separation between classical an quantum due to the lack of parameters needed to encode a computational task. We will rather show first that a sufficiently fast (e.g. linear) increase of the entanglement entropy makes it hard -and for MPS based approaches [15,16] impossible -to even approximate the state of the system efficiently. That such a linear increase can occur even in the simplest one-dimensional systems was observed in [17] based on arguments of conformal field theory together with numerical diagonalization of chains up to 100 sites. The main part of the present paper is devoted to a rigorous proof of this observation.

The hardness of simulating time evolution
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 invariant product state of spin-1 2 particles, and let it evolve under a translational invariant nearest neighbor Hamiltonian. For this system, we prove a linear lower 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 been studied in [18], where similar effects were observed. We then combine this linear bound with the continuity inequality for the von Neumann entropy which leads to a linear lower bound on the block entropy of any approximation of the state [19]. As the computational effort needed to deal with a general MPS grows exponentially in the entropy, it would thus increase exponentially 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 2 logarithm. The natural logarithm will be denoted by ln.

Setup and entropy scaling
Theorem 1 Consider a chain with an odd number N of spins with periodic boundary conditions, corresponding to a Hilbert space (C 2 ) ⊗N , which at time t = 0 is in the state |ψ(0) = | ↑ . . . ↑ ≡ |1 . . . 1 and evolves under the Ising with ρ L (t) = tr L+1,...,N |ψ(t) ψ(t)|. Then for N ≥ 20, L ≥ 10, 4 ≤ t ≤ eL/4, and t ≤ N/5, We postpone the proof of the Theorem to Sec. 3. Note that although the lower bound is independent on L -the intuition behind being that entanglement arises at the boundaries -the maximally attainable entanglement is proportional to the length of the block due to the constraint t ≤ eL/4.

Exponential scaling of bond dimension
Let us now see what the entropy scaling implies for an ansatz state |φ ≡ |φ(t) which we use to approximate |ψ ≡ |ψ(t) . Let ρ L ≡ ρ L (t) = tr B |ψ ψ|, σ L ≡ σ L (t) = tr B |φ φ| be the reduced states on a block A of length L (B is the complement of A). The error we make in this approximation is where we have used that the trace norm η 1 = tr |η| is contractive under the partial trace. We now exploit Fannes' continuity inequality for the von Neumann entropy in its improved version by Audenaert [20], is the binary entropy, and thus Using the lower bound on S(ρ L (t)) ≡ S L (t) of Theorem 1, we have (under the corresponding conditions on N, L, and t) that To get an optimal bound, we choose L as small as possible for the given t, L = 4t/e, which implies the new constraint t ≥ 5e/2, and obain In order to turn this into a lower bound on the size of an MPS |φ(t) used to approximate the state at time t, note that the entropy of any block is naturally upper bounded by the rank of the reduced state ρ L , which in turn is at most D 2 . Here, D is the so-called bond dimension of the MPS, which is polynomially related to the complexity of the MPS [14]. Then, (2) gives i.e. the bond dimension -and thus the resources for a straight forward MPS simulation of the time evolution -will scale exponentially with time as soon as the accuracy ǫ < ǫ 0 = 2e/3π ≈ 0.577. Note that the same argument holds similarly for lattices in higher spatial dimensions. As long as the entropy S(ρ A (t)) of a region A has a linear leading term, there is an accuracy ǫ 0 beyond which S(σ A (t)) has to grow linearly as well. Hence, approaches based on projected entangled pair states (PEPS) [21] or variants thereof are burdened with the same problem.

Proof of Theorem 1
We now prove Theorem 1. Let us briefly sketch the main steps: First, we map the spin chain to a one-dimensional system of free fermions, which can be solved exactly, and bound the error made by going to the thermodymamic limit N → ∞. The entropy of a contiguous subblock of length L is equal to the entropy of the corresponding fermionic modes; using a parabola as a bound on the binary entropy we can lower bound the block entropy by a quadratic function in the correlation matrix. Combining this with the exact solution for the thermodynamic limit, we obtain a lower bound on the block entropy which is essentially of the form with J n Bessel functions of the first kind. In the appendix we bound this sum by deriving a lower bound for L = ∞ and an upper bound on the error made 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 ≥ 20, L ≥ 10, 4 ≤ t ≤ eL/4, and t ≤ N/5. In addition we can assume L ≤ N/2 as the overall state is pure.

Diagonalization of the Hamiltonian
We start by mapping the spin operators to fermionic Majorana operators via the Jordan-Wigner transformation We will refer to c 0 , . . . , c N −1 as the position and to c N , . . . , c 2N −1 as the momentum operators; each pair c k , c N +k defines a fermionic mode. For the initial state |ψ(0) = |1 · · · 1 we have

and the Hamiltonian
is equal to the spin Hamiltonian (1) of Theorem 1 on the −1 eigenspace of Π = σ (1) z · · · σ (N ) z , while for the +1 eigenspace, the coupling between spins N and 1 has opposite sign. Since Π|ψ(0) = −|ψ(0) for odd N and [H, Π] = 0, (5) indeed describes the evolution of the Hamiltonian of Theorem 1. Note that the following results equally hold for even N where they describe an Ising system with a flipped coupling across the boundary.
Since the initial state is a Gaussian and the Hamiltonian a quadratic function in the Majorana operators, the system remains in a Gaussian state throughout its evolution and is thus completely characterized by its second moments, the antisymmetric correlation matrix [22,23] Γ kl = − i 2 tr (ρ[c k , c l ]) which for the initial state is where we defined the Hamiltonian matrix H by We now apply Fourier transformations F kl = e 2πi kl/N / √ N to both the position and the momentum subspace. Thereby, Γ 0 and H become block diagonal with blocksΓ where φ k = 2π k N . The evolution (6) becomeŝ such that the CM at time t (written in modewise ordering) is Solving (9), we find that This is the well-known exact solution of the Ising model as found, e.g., in [17,24].

The thermodynamic limit
In order to facilitate the evaluation of f n and g n , we consider the limit N → ∞, i.e. we replace We label the functions obtained thereby by f ∞ n and g ∞ n , respectively. The error induced by taking the limit can be bounded using the following theorem.
Theorem 2 (cf. [25]) Let h : R → R be analytic and 2π-periodic. Then there exists a strip D = R × (−a, a) ⊂ C with a > 0 such that h can be extended to a bounded holomorphic and 2π-periodic function h : D → C. Let M = sup D |h|. Then, It is straightforward to see that for any a > 0, the addends in both (12)  Applying the above theorem for a = 2.9, and taking into account that n ≤ N/2, t ≥ 4, and N ≥ 20, one obtains that Taking the limit in (12) and (13), we find and [cos(2nα) + cos((2n + 2)α)] cos(2t sin α) dα + I n where I 0 = 1 2 , I −1 = − 1 2 , and I n = 0 otherwise. Here, J k are Bessel functions of the first kind, and we have used the recurrence relation 1 2 [J n−1 (z) + J n+1 (z)] = nJ n (z)/z.

Lower bound for the block entropy
The reduced state of the first L sites is again Gaussian and its correlation matrix is given by the 2L × 2L upper diagonal subblock of Γ t , which we label by A t . Its entropy can be computed from the normal mode decomposition [26] where O ∈ SO(2L), as We lower bound h(x) ≥ 1 − x 2 as in [23,27] and thus obtain the bound Writing in the 1, . . . , L and L + 1, . . . , N partitioning (thus A ≡ A t etc.) and using that the purity of the overall state implies Γ 2 t = −1 1 [22,23], we find A 2 − CC T = −1 1 and thus We further bound the sum by only considering the lower left and the upper right corner of C, i.e., all entries γ n in (10) We have =:Bn and with |J n | ≤ 1/ √ 2 for |n| ≥ 1, Thus, the total error made in (17) is bounded by where the latter follows using (14) together with L ≤ N/2, 4 ≤ t ≤ N/5, and N ≥ 20. Thus, where we have used Eqs. (15) and (16), J −n (z) = (−1) n J n (z), |J 1 (2t)| ≤ 1/ √ 2, and t ≥ 4. We now bound the sum in (18) using Lemma 1 and 2 (see Appendix ) and obtain using 4 ≤ t ≤ eL/4 and L ≥ 10, which proves Theorem 1.
Note that different constraints on t, L, and N can be chosen, resulting in different corrections to the leading terms in the expansion of S L (t).

Discussion
We gave a rigorous proof to the observation that the entanglement entropy in a simple one-dimensional system increases essentially linearly. Although the proof is tailored to a single exactly solvable model, we think that this is a widespread property, thus being a requirement which has to be met by every classical simulation method. Of course, there always exist specific cases with tailored representations (e.g. the quasi-free Fermions used here or Clifford circuits) close to which one can meet this requirement [28,29] -even MPS allow for a description of specific states with only polynomial effort in the entanglement entropy (e.g., if the matrices are tensor products [30]). Yet, these methods seem much more restricted in their applicability than MPS, so that based on present knowledge it seems unlikely that there is an efficient and generally applicable classical simulation method. Note that similar results can be found by considering any Rényi entropy with α > 1 instead of the von Neumann entropy [19]. Also observe that in our argument we required a good global approximation of the state which is a priori not necessary for obtaining good results for observables with only local non-trivial support.
Note that the proof can be extended straightforwardly to sums ∞ n=K .