Observations Outside the Light-Cone: Algorithms for Non-Equilibrium and Thermal States
aa r X i v : . [ qu a n t - ph ] J a n Observations Outside the Light-Cone: Algorithms for Non-Equilibrium and ThermalStates
M. B. Hastings , Center for Nonlinear Studies and Theoretical Division,Los Alamos National Laboratory, Los Alamos, NM, 87545 Kavli Institute of Theoretical Physics, University of California, Santa Barbara, CA, 93106
We apply algorithms based on Lieb-Robinson bounds to simulate time-dependent and thermalquantities in quantum systems. For time-dependent systems, we modify a previous mapping toquantum circuits to significantly reduce the computer resources required. This modification is basedon a principle of “observing” the system outside the light-cone. We apply this method to studyspin relaxation in systems started out of equilibrium with initial conditions that give rise to veryrapid entanglement growth. We also show that it is possible to approximate time evolution undera local Hamiltonian by a quantum circuit whose light-cone naturally matches the Lieb-Robinsonvelocity. Asymptotically, these modified methods allow a doubling of the system size that one canobtain compared to direct simulation. We then consider a different problem of thermal properties ofdisordered spin chains and use quantum belief propagation to average over different configurations.We test this algorithm on one dimensional systems with mixed ferromagnetic and anti-ferromagneticbonds, where we can compare to quantum Monte Carlo, and then we apply it to the study ofdisordered, frustrated spin systems.
I. INTRODUCTION
Matrix product and density-matrix renormalization group methods provide one of the most powerful ways ofsimulating one-dimensional quantum systems. In addition to ground state properties[1], they have been extendedto thermal and open systems[2] and dynamical problems[3]. The reason for the success of these algorithms is thatin many cases the appropriate quantum state can be well-approximated by a matrix product state, giving a verycompact representation of the state of the system.In some cases, we even have theorems that quantify the accuracy of matrix product states. For ground states oflocal quantum systems with a gap, the ability to represent ground states as matrix product states follows from boundson the entanglement entropies[4]. Related results are available for thermal systems[5, 6] and for non-equilibrium statesobtained by starting with a factorized state and evolving under a local Hamiltonian for a time t ; in the first casethe bond dimension needed to get a good matrix product approximation to the desired state scales exponentiallyin β , while in the second case it scales exponentially in t [7]. The two results [5, 7] are constructive proofs, whichgive an algorithm to find the matrix product state. All of these constructive proofs rely heavily on Lieb-Robinsonbounds[8, 9, 10, 11].In contrast to these constructive proofs, the matrix product algorithms used in practice are variational: they involveoptimizing over different matrix product states to find the best one. This works well because in many practical casesthe entanglement grows much more slowly than the upper bounds set by theory. For example, systems described byconformal field theory have an entanglement entropy growing only logarithmically with system size[12], while the arealaw bound[4] gives no useful result for these systems due to the absence of a gap. Similarly, for many initial conditions,evolution under a local Hamiltonian gives an entanglement entropy growing only logarithmically in time[13, 14, 15]while the theoretical upper bound gives an entanglement entropy growing linearly in time[16].In this paper, we argue that there are many situations in which algorithms based on Lieb-Robinson bounds arethe best technique. We first look at the case of time evolution in systems out of equilibrium. Here, there are initialconditions for which the entanglement entropy is known to grow linearly in time[17], in accordance with conformalfield theory predictions[18]. Roughly speaking, logarithmic entropy growth tends to occur in cases where we candivide the chain into a small number of subchains such that the initial state is an eigenstate of the Hamiltonian oneach subchain; for example, starting an XXZ spin chain in a state in which all the spins on the left half of the chainare up and all the spins on the right half of the chain are down leads to a logarithmic entropy growth[15]. On theother hand, the linear entropy growth tends to occur in cases where the initial state differs from an eigenstate ofthe Hamiltonian on every subsystem of the full spin chain. For example, starting an XXZ spin chain in an initialcondition in which the spins alternate between up and down (a Neel state, the ground state when the Ising term isthe only term in the Hamiltonian) leads to a linear entropy growth[17].For system with linear entropy growth, matrix product methods will require a bond dimension growing exponentiallyin time to obtain accurate results. At this point, both variational matrix product and constructive, Lieb-Robinson-based methods require resources growing exponentially in time. The question, then, is how to obtain the smallestexponential. To some extent, the Lieb-Robinson based methods (such as [7] and the methods below) are “worst-case”:the bond dimension depends on theoretical upper bounds for arbitrary local Hamiltonians, while the matrix productmethods can adaptively find better representations. On the other hand, there are some disadvantages to matrixproduct methods. To get a rough idea of the resources required, let us consider a system of N spins, each of spin-1 / t involves writing the initialcondition down in some basis and then directly simulating it (we discuss below different methods for doing this),requiring resources scaling as 2 N t using sparse matrix methods. A matrix product algorithm can avoid truncationerror for this system by using matrices with bond dimension 2 N/ [19, 20]. However, the algorithm must then performsingular value decompositions and eigenvalue calculations, taking a time which scales as the cube of these matrices,and hence of order 2 N/ , which is slower. Of course, the matrix product methods are really only useful if they areable to represent the system with a smaller bond dimension. In this case, though, even if such a representation exists,the algorithm must find it, and this can pose a problem. The algorithm for simulation of non-equilibrium systemsdepends on breaking the time evolution into a series of small Trotter steps; if the Trotter steps are too long this canlead to Trotter error, while if they are too short, the truncation error can grow rapidly: even if there is a good statethe algorithm may not find it[15].We can use this scaling of the difficulty with N to get an idea of the scaling of computation effort with time, usinga Lieb-Robinson bound on the group velocity, v LR . Suppose we wish to compute the expectation value of a localobservable, such as a spin on a site, at a time t f , starting from a factorized state at time t = 0. The number of spinsin the past light-cone of this spin is 2 v LR t , and the Lieb-Robinson bounds imply that the effect of spins outside thelight-cone is exponetially small. Thus, it suffices to simulate only the dynamics of the N = 2 v LR t f + O (log( ǫ )) spinsclosest to the given spin in order to compute the expectation value to an accuracy ǫ . This requires an effort scalingexponentially in time as t v LR t . Similarly, if the entanglement entropy grows linearly in time, the matrix productmethods also require an effort scaling exponentially in time.Our main result in this paper is the light-cone quantum circuit algorithm, an application of the Lieb-Robinsonbound that allows us to simulate the evolution of local observables with resources growing asymptotically as only N t v LR , by some statistical sampling. This allows twice as large systems as the direct method. We analyze the entropygrowth in these systems and argue that matrix product methods are also less efficient for long time simulation. Weapply the light-cone quantum circuit algorithm then to the problem of spin relaxation in spin chains started in theNeel state. The physical idea behind the light-cone quantum circuit algorithm is as follows: to find the state of agiven spin at a time t f , we only have to track the dynamics within the past light-cone of the spin. For times t closeto zero, the past light-cone includes roughly 2 v LR t f spins, but at these early times the entanglement is small andhence the computational effort should be less. For times t close to t f , the past light-cone includes few spins and henceshould be easier to simulate.The paper is organized as follows. We first derive the light-cone quantum circuit algorithm. We then apply it tospin relaxation, and study oscillations of the central spin, decay of the envelope of the oscillations, and also seeminglyrandom oscillations of the central spin in chains where boundary effects become important. We then derive a relatedquantum circuit method, the corner transfer quantum circuit which may be useful for studying the evolution of globalobservables in highly entangled non-equilibrium states.We then turn to a different problem, presenting one other application of Lieb-Robinson methods, using the quantumbelief propagation algorithm[22] to study thermal states in disordered systems. The quantum belief propagationalgorithm explicitly constructs a matrix product state for a thermal quantum system. While it manipulates operators,rather than states, and thus can be computationally expensive, it has other advantages. It has no Trotter error, makingit fast and accurate at high temperatures: it can obtain quantities such as the susceptibility peak to higher accuracyusing fewer resources than methods such as transfer matrix DMRG[23], although at low temperatures it breaks down,with the resources required scaling exponentially with the temperature. In this case, the exponential scaling withthe temperature is again related to a linear relationship between a time scale, in this case β = 1 /T , and a lengthscale. It can be applied to random systems, where transfer matrix DMRG cannot be used because of a lack oftranslation invariance. A good test of variational matrix product[2] methods on this kind of system is lacking, sowe cannot compare here. We apply the quantum belief propagation algorithm to two random systems, one withoutfrustration where we can compare to quantum Monte Carlo and one with frustration where Monte Carlo methods arenot applicable. II. QUANTUM CIRCUIT METHODS
In this section we present the various quantum circuit methods. We begin by reviewing previous work, and thenderive the light-cone quantum circuit algorithm, and apply it to a problem of spin relaxation. We then use previousresults on the entanglement entropy growth to estimate the computational resources required for different approachesto this problem, and finally we present the corner transfer quantum circuit method, an extension which allows accessto global quantities.
A. Background
To understand our algorithm, we first review the ideas in [7] which gives a construction of a matrix productoperator approximation to the time evolution operator, exp( − iHt ), using resources exponential in t . We consider alocal Hamiltonian H = X i h i , (1)where each i acts on sites i, i + 1. This Hamiltonian obeys a Lieb-Robinson bound: given any operator O which hassupport on set of sites X , the operator exp( iHt ) O exp( − iHt ) can be written, with exponentially small error, as anoperator acting on the set of sites i within distance v LR t of X , where v LR is the Lieb-Robinson group velocity.To simulate the system for a time t , we divide the system into blocks of length l , where l is slightly larger than2 v LR t (the error in the approximation will be exponentially small in l − v LR t ). We then let H = H + H ′ where H is the sum of the Hamiltonians on each block and H ′ is the Hamiltonian connecting the blocks: H = X k i ≤ ( k +1) l − X i = kl +1 h i , (2) H ′ = X k h kl , (3)where the sum ranges over integers k .We then write exp( − iHt ) = (cid:16) T exp[ − i Z t exp( − iH t ′ ) H ′ exp( iH t ′ )] (cid:17) exp( − iH t ) , (4)where T denotes that the exponential is time-ordered. The operator exp( − iH t ) is equal to the product Q k U k , where U k is a unitary operator acting on sites kl + 1 , kl + 2 , ..., kl : U k = exp( − i i ≤ ( k +1) l − X i = kl +1 h i t ) . (5)Using the Lieb-Robinson bounds, we can approximate exp( − iH t ′ ) h kl exp( iH t ′ ) by an operator h kl ( t ′ ) loc whichhas support on sites kl − l/ , ..., kl + l/ − | t ′ | ≤ t and for l > v LR t . Then the operator T exp[ − i R t exp( − iH t ′ ) H ′ exp( iH t ′ )] can then be approximated by a product Q k V k where V k = T exp[ − i Z t h kl ( t ′ ) loc ] . (6)The key in this construction is that the intervals kl − l/ , ..., kl + l/ − k . Thisconstruction expresses the time evolution as a quantum circuit:exp( − iHt ) ≈ Y k V k Y k U k . (7)The support of the operators U k , V k and the circuit is shown in Fig. 1.Precise error bounds can be given using Lieb-Robinson bounds for the error in Eq. (7). To get an error of order ǫ in the propagator (7), we only need to take l = 2 v LR t + O (log( N/ǫ )). In what follows, we will not make detailed errorestimates, since Lieb-Robinson error estimates are fairly simple and are by now standard in the literature; when we saythat it suffices to take a length scale “of order” v LR t to obtain an approximation to a given local quantity, we mean thatby taking the length scale v LR t + O (log( N/ǫ )) we can obtain an error of order ǫ in the state; when we are computing U U U1 2 3
V1 V2
44 164 16 4 4 16 4 4 16
FIG. 1: Support of the operators U k , V k and the quantum circuit for N = 12, l = 4. The small numbers at the bottomof operators U k , V k represent the bond dimensions required to represent these operators as matrix product operators; themaximum product of these across any bond is 16 = 4 l/ . expectation values of local quantities, to obtain an error of order ǫ we need a length scale v LR t + O (log( v LR t/ǫ )), sothat the error bound does not depend on N in this case.Suppose we want to apply the quantum circuit procedure to compute the time evolution of some state Ψ . Forsimplicity, let Ψ be a factorized state (later we discuss the case where Ψ is a matrix product state both in thisprocedure and using our algorithm, and we find that using the idea of “observation” discussed below the case ofmatrix product state initial conditions presents no additional difficulty). The operator U k is an operator acting on 2 l sites. Any such operator can be written as a matrix product operator with a bond dimension equal to 4 l/ = 2 l [19].This maximum bond dimension is achieved halfway across the interval of length 2 l , while at any point a distance d from one of the ends of the interval the bond dimension is only 2 d , as shown in Fig. 1.The same holds for V k , and so the maximum product of bond dimensions across any bond is 2 l (a slightly worseestimate of 2 l was found for this construction in [7] since the fact that the bond dimension may vary with positionwas not taken into account).There are several problems, however, with implementing the above method in practice, and it is these problemswhich we overcome with our methods, the light-cone quantum circuit method and the corner transfer quantum circuitmethod. The first method is most appropriate for computing local quantities (such as a spin or energy expectationvalue) while the second method is most appropriate for finding a good global approximation to the ground state.The first problem is that operator equations of motion are computationally expensive in practice. To this end, wewill modify the procedure to deal only with state vectors, rather than operators. The second problem is that thevelocity of the quantum circuit does not obviously match the Lieb-Robinson velocity, in the following sense: given arbitrary operators U k , V k , supported as described above, the product Q k V k Q k U k can propagate information bya distance of 2 l in each time step. Since l is roughly v LR t , this means that such a quantum circuit could have aneffective velocity roughly 2 v LR . Of course, the operators U k , V k are not arbitrary operators, but still we would like tofix this problem; we will show how to do this with the corner transfer quantum circuit method below which also leadsto improved estimates on the maximum matrix product state dimension needed.However, the real problem with this method is that it doesn’t lead to any improvement over naive simulation whenit comes to computing local observables. The main problem we consider in this section is the following: we starta spin chain at time t = 0 in a factorized state Ψ , and then evolve under a local Hamiltonian for to a final time t f , at which point we wish to compute some local observable, such as S zi , the z -expectation value of spin i . By theLieb-Robinson bounds, we can approximate S zi ( t ) by considering only a subchain of the full chain: we consider onlysites i − l, ..., i + l , where l is slightly larger than v LR t f . We then define Ψ ′ to be the appropriate factorized stateon this subchain and evolve Ψ ′ for a time t f using the Hamiltonian H ′ acting on the subchain. We then compute h Ψ ′ | exp( iH ′ t f ) S zi exp( − iH ′ t f ) | Ψ ′ i . Using sparse matrix methods to compute the time evolution of Ψ ′ , this requiresa computational effort of order l l , which scales as 2 v LR t f . The quantum circuit method discussed above would alsorequire simulations on intervals of length 2 l , and hence leads to no improvement when computing this local quantity. M−l −l/2 −1 0 1 l/2 lH L H RH’L H’RH
FIG. 2: Support of the operators H L , H R , H M , H ′ L , H ′ R . B. Light-Cone Quantum Circuit Algorithm
We now show how to reduce the computational effort to an amount of order 2 v LR t f , allowing the time scale tobe twice as large, by combining Lieb-Robinson bounds with statistical sampling. While we focus on this section onstarting in a factorized state, we later discuss the case of starting in a matrix product state.We now derive our algorithm, which we call the light-cone quantum circuit method as it avoids keeping track ofcertain degrees of freedom outside the light-cone by making certain “observations” to reduce the computational effort.We first define a subchain of length 2 l + 1 and an initial state Ψ ′ on that subchain as above. We label the sites in thesubchain by − l, ..., , ..., + l . We let H ′ be the Hamiltonian on the subchain and we write H ′ = H L + H R + H B , (8)where H L acts on the left half of the chain (sites − l, ..., − H R acts on the right half of the chain, and the boundaryHamiltonian H B acts on sites − , ,
1. Thus, [ H L , H R ] = 0. We define H M to act on the middle half of the chain: itis supported on sites − l/ , ..., + l/ l even for simplicity). See Fig. 2.Using the Lieb-Robinson bounds in the same way as in the quantum circuit method above we can approximate thetime evolution for a time t i = t f / ′ ( t f /
2) = exp( − iH ′ t f / ′ ≈ exp( − iH M t f /
2) exp[ − i ( H M − H B ) t f /
2] exp[ iH L t f /
2] exp[ iH R t f / ′ . (9)Note that the decomposition (9) is only good for times of order t f /
2. We wish to compute the expectation value of S z , or any other observable on site 0, at time t f . Again using the Lieb-Robinson bounds, the expectation value ofthis is approximately equal to h Ψ ′ ( t f / | exp( iH M t f / O exp( − iH M t f / | Ψ ′ ( t f / i . (10)Combining Eqs. (9,10) and our use of the Lieb-Robinson bounds to approximate the expectation value of O on thefull chain by its expectation value on the subchain, we have h Ψ | O ( t ) | Ψ i ≈ h ˜Ψ | O | ˜Ψ i , (11)where ˜Ψ = exp( − iH M t f / ′ ≈ exp( − iH M t f ) exp[ − i ( H M − H B ) t f /
2] exp[ iH L t f /
2] exp[ iH R t f / ′ . (12)The operator H M − H B is a sum of two operators, H ′ L and H ′ R , where H ′ L acts on sites − l/ , ..., − H ′ R acts onsites +1 , ..., + l/
2. Therefore,˜Ψ = exp( − iH M t f / ′ ≈ exp( − iH M t f ) (cid:16) exp[ iH ′ L t f /
2] exp[ − iH L t f / (cid:17)(cid:16) exp[ iH ′ R t f /
2] exp[ − iH R t f / (cid:17) Ψ ′ . (13)Eq. (13) is not yet useful computationally, since it will require an effort of order lt f l to compute the evolution ofthe state Ψ ′ . We now describe the light-cone quantum circuit method to compute the expectation value: first, writeΨ ′ = Ψ L ⊗ Ψ C ⊗ Ψ R , (14)where Ψ L , Ψ R are states on the left and right half of the chain, and Ψ C is a state on the center site. We compute thestates Ψ ′ L = (exp[ iH ′ L t f /
2] exp[ − iH L t f / (cid:17) Ψ L , (15)and Ψ ′ R = (exp[ iH ′ R t f /
2] exp[ − iH R t f / (cid:17) Ψ R , (16)which requires an effort of order lt f l . Next, we introduce a complete orthonormal basis of states on the sites − l, ..., − l/ −
1, which we label φ L ( α ), and another a complete basis of states on sites l/ , ..., + l labelled φ R ( α ).Then we decompose Ψ ′ L and Ψ ′ R as: Ψ ′ L = X α A ( α ) φ L ( α ) ⊗ ξ L ( α ) , (17)Ψ ′ R = X α A ( α ) φ R ( α ) ⊗ ξ R ( α ) , where ξ L ( α ) is some normalized state on sites − l/ , ..., − ξ R ( α ) is some normalized state on sites +1 , ..., + l/ φ L ( α ) need not be eigenvectors of any reduced density matrix and the states ξ L ( α ) need not be orthogonalto each other, as Eq. (17) is not a Schmidt decomposition. Thus, from Eqs. (13,15,16,17), h ˜Ψ | O | ˜Ψ i = X α L X α R | A ( α L ) | | A ( α R ) | E ( O, α L , α R ) , (18)where E ( O, α L , α R ) = h ξ L ( α L ) ⊗ Ψ C ⊗ ξ R ( α R ) | exp( iH M t f ) O exp( − iH M t f ) | ξ L ( α L ) ⊗ Ψ C ⊗ ξ R ( α R ) i . (19)Eq. (18) is at the heart of the light-cone quantum circuit approach. Numerically we proceed as follows: first, wecompute | A ( α L ) | and | A ( α R ) | for all α L , α R . This requires an effort of order 2 l . We then do a statistical sampling:we randomly pick an α L and an α R according to the probability distributions | A ( α L ) | and | A ( α R ) | , and computethe average E ( O, α L , α R ). We repeat this procedure many times, to average over different choices of α L , α R .The computational effort required then scales as only l l , or roughly t v LR t as claimed. Asymptotically this allowsdouble the time. The time scales linearly in the number of iterations of statistical sampling, which we denote N it .However, if O has bounded operator norm, then E ( O, α L , α R ) has bounded moments, and so by the central limittheorem, the number of iterations required still scales only polynomially in the error. C. Results on Non-Equilibrium Dynamics
We now discuss results from this method, as well as some implementation details. We consider evolution under theXXZ Hamiltonian H = X i ( S xi S xi +1 + S yi S yi +1 + ∆ S zi S zi +1 ) . (20)For ∆ = 0, this problem can be mapped to free fermions by a Jordan-Wigner transformation and solved exactly. Weuse this as a check on our results later.We used as a starting point the Neel state, with spins alternating up and down, and we computed the timedependence of S z ( t ) for the central spin. The main numerical effort is to compute the evolution of a state undera Hamiltonianm, which we did using a combination of short steps with a series method. For example, to computeexp( − iH L t )Ψ L we divide the time t into shorter intervals of time t , and compute exp( − iH L t )Ψ L = Ψ L − it H L Ψ L − ( t H L / L + ... , keeping a fixed number of terms in this series. We then repeat this procedure ( t/t ) times. Toobtain negligible error for a chain of 20 sites with t = 1 required going to roughly 40-th order for | ∆ | ≤
1, while for∆ = 2 slightly longer series were required. A more sophisticated way of doing the time evolution would be to builda tridiagonal Hamiltonian in the Krylov space spanned by Ψ L , H L Ψ L , ..., H kL Ψ L for some k , and then evolve exactlywith this Hamiltonian[25].Another important point of numerical simulation is the use of symmetries. We can choose the states φ L ( α ) and φ R ( α ) to be eigenstates of total S z . Then, since the state Ψ ′ is an eigenstate of total S z , the state ξ L ( α L ) ⊗ Ψ C ⊗ ξ R ( α R )is also an eigenstate of total S z , which allows us to use symmetries when computing the evolution of the state. Sincemost of the numerical time is consumed statistically sampling E ( O, α L , α R ), we build the sparse matrix for theHamiltonian H M in each spin sector once, before doing the sampling, and then run the sampling.It is also possible, although we did not implement it, to take into account reflection symmetry. Since α L and α R are chosen independently, the state ξ L ( α L ) ⊗ Ψ C ⊗ ξ R ( α R ) does not have reflection symmetry. However, if both theHamiltonian H M and the operator O = S zi have reflection symmetry, then it is useful to write ξ L ( α L ) ⊗ Ψ C ⊗ ξ R ( α R ) = Ψ S ( α L , α R ) + Ψ A ( α L , α R ) , (21)where Ψ S , Ψ A are symmetric and anti-symmetric states. Then, E ( O, α L , α R ) = h Ψ S ( α L , α R | exp( iH M t ) O exp( − iH M t ) |i Ψ S ( α L , α R ) i (22)+ h Ψ A ( α L , α R | exp( iH M t ) O exp( − iH M t ) |i Ψ A ( α L , α R ) i , and so we can statistically sample one of the two terms on the right-hand side of Eq. (22). Note that on each iterationwe randomly choose an α L and an α R and then randomly choose a term in Eq. (22), rather than repeatedly samplingEq. (22).As the algorithm is described above, the initial computation of the states Ψ L ( α ) and Ψ R ( α ) depends on the finaltime. For each final time t , we have to compute a new set of states Ψ L ( α ) and Ψ R ( α ) and then do the statisticalsampling. However, in fact, we can speed the algorithm up at a slight cost in accuracy: we fix a given t f and on eachstatistical sample we compute the stateexp( − iH M t f ) | ξ L ( α L ) ⊗ Ψ C ⊗ ξ R ( α R ) i , (23)to evaluate the expectation value in Eq. (19). We then act on this state with the operator exp( iH M δt ) for some small δt . We then use this new state to compute an approximation to the expectation value at time t f − δt . We then acton that state with exp( − iH M δt ) to compute an approximation to the expectation value at time t f − δt , and so on.Since the computational cost of performing the time evolution of a state under a Hamiltonian is proportional to thetime evolved, these additional steps are relatively cheap, for δt << t f . There is a small cost in accuracy: in general,to compute expectation values at a time t , we can do the initial evolution for a time t i , and then evolve further fortime t − t i . To make the effect of boundary conditions as small as possible, we would like to have both t i and t − t i as small as possible, which is why above we chose to evolve for a time t i = t f /
2. However, if δt << t f , then we arenot far away from the ideal choice of t i by initially evolving for time t f / t f / − δt . Wefollowed this procedure in the numerical work below, with δt = 0 .
25 and taking t f to be spaced with integer stepsusing t = 1. This accounts for some of the slight kinks in the curves after every integer value of t .The spin-wave velocity of Hamiltonian (20) for ∆ ≤ v sw = ( π/
2) sin( θ ) /θ, (24)where cos( θ ) = ∆[24]. On the other hand, in the development above, we used a Lieb-Robinson velocity v LR , where v LR ≥ v sw . (25)Using the Lieb-Robinson bounds, we showed that we could accurately simulate for a time t using length scales l = v LR t .As mentioned above, we do not give precise error estimates, but it is not hard to give rigorous estimates of the error.However, the Lieb-Robinson bound is actually fairly conservative because of Eq. (25), and thus in practice lengthscales v sw t suffice to get good results.We begin by illustrating results for the XY chain, with ∆ = 0. In Fig. 3, we illustrate results from exact simulationsof the XY chain for various sizes. We consider chains with open boundary conditions with N = 35 , ,
101 and achain with periodic boundary conditions with N = 36. For the chains with open boundary conditions, we plot theaverage of the central spin as a function of time, while for the periodic boundary conditions we plot the average ofan arbitrarily chosen spin as a function of time. The large exact simulations are possible because this chain can bemapped to free fermions by a Jordan-Wigner transformation. Later we will present comparison of these results tolight-cone quantum circuit results. For now, we discuss aspects of Fig. 3 which show the influence of boundary effects.In the range of times, the chain with N = 101 shows no effect of the boundary. There are oscillations of the spin withfrequency ω roughly 2, with an envelope decaying as 1 / √ t so as h S zi i ∼ √ t cos( ωt + θ ) , (26) Time -0.4-0.200.20.4 < S z > f o r ce n t r a l s p i n N=101N=51N=35N=36
FIG. 3: Time dependence of h S z i for the central spin as a function of time. Curves are N = 101 (black), N = 51 (red), N = 35(green), and N = 36 (blue). for θ ≈ (3 / π . Simulations on longer chains show that as long as N is less than t/
2, the decaying oscillations ofEq. (26) continue to hold. In regard to the 1 /t / decay of oscillations found here numerically, it is interesting to notea similar power-law decay found for a different system of free bosons, where fluctuations about a maximal entropystate were proven to decay at least as fast as 1 /t / [21].For N larger than t/
2, the boundary conditions become important, as one can see in the curves for N = 35 and N = 51, which start to show deviations from the N = 101 curve for times roughly 16 and 22 respectively. This is nosurprise, since the distance between the central spin and the boundary is N/
2, and v sw = 1 for ∆ = 0. We will seelater for ∆ = 0 that in general boundary effects become important for N/ tv sw . Interestingly, periodic boundaryconditions offer no improvement, since the N = 36 periodic curve deviates at the same time as the N = 35 opencurve.Another interesting effect is that once the boundary conditions become important, the expectation value showswild oscillations, which no longer decrease in magnitude. This may be a consequence of the fact that the chain isintegrable. It would be interesting to see for a non-integrable system whether such oscillations occur or not; thesimplest statistical assumption for a non-integrable system is that the state at long times would be a random purestate satisfying the conservation laws of total S z and total energy, and thus the expectation value of h S zi i for any i would show only exponentially small fluctuations about the average spin.We now consider the application of the light-cone quantum circuit method to this chain. We considered l = 18 and20, and did N it = 1000 iterations of statistical sampling, as shown in Fig. 4. The results for exact simulations with N = 35 and N = 101 are also shown for comparison. We see that the light-cone quantum circuit method with l = 18is accurate over the same range of times as the exact simulation with N = 35, while the light-cone quantum circuitmethod with l = 20 improves on this result. By increasing l from 18 to 20, we increase the range of times by roughly2, while increasing N by 2 would increase the range of times by only 1.There are statistical fluctuations in the light-cone quantum circuit results in Fig. 4, due to random fluctuations in E ( O, α L , α R ) for different choices of α L , α R . In Fig. 5 we plot the rms fluctuation in E ( O, α L , α R ) as a function oftime, sampling this expectation value with the probability distribution | A ( α L ) | | A ( α R ) | . The results in Fig. 4 are anaverage over N it = 1000 samples, so the spread on each data point in Fig. 4 is equal to 1 / √ .
005 in the worst case. A very interesting point is that for t ≤
5, the rms fluctuation isnegligible, while for times t ≥
9, the rms fluctuations are up to their full value, with a sharp bend in the curve (plottedon a log scale) around t = 8 or t = 9. This is a consequence of the maximum spin-wave velocity: the influence of theregions on sites − l, ..., − l/ l/ , ..., + l takes a time of order l/ v sw to reach the central region. There are sillsome fluctuations for t ≤ l/ v sw = 9, but they decay rapidly until they are negligible at very short time. Time -0.4-0.200.20.4 < S z > f o r ce n t r a l s p i n N=101l=18N=35l=20
FIG. 4: Time dependence of h S z i for the central spin as a function of time for ∆ = 0. Exact curves are N = 101 (black) and N = 35 (green). Light-cone quantum circuit curves are l = 18 (red) and l = 20 (blue). Time R M S f l u c t u a ti on s i n E FIG. 5: RMS fluctuations in E ( O, α L , α R ) as described in text as a function of time for l = 18. We then applied the light-cone quantum circuit method to chains with ∆ = 0 . , , N = 20. For larger ∆ = 0 . , t − . for ∆ = 0 .
5. For ∆ = 2, no oscillations are seen. The effects of statistical noise are much more noticeable, sincethe magnitude of the spin is much less. All simulations were done with N it = 1000 statistical samples, except thesimulation with ∆ = 0 . , l = 22 has only N it = 250 samples, and the simulations with ∆ = 1 , l = 16 and ∆ = 1 , l = 180 Time -0.4-0.200.2 < S z > f o r ce n t r a l s p i n exact, N=20l=18l=20l=2210 15 20-0.0500.05 FIG. 6: Time dependence of h S z i for the central spin as a function of time for ∆ = 0 .
5. Exact curve is N = 20 (black).Light-cone quantum circuit curves are l = 18 (red), l = 20 (blue), and l = 22 (green). had N it = 10000 and N it = 3000 samples respectively. The spin-wave velocity is larger for these chains, so thesimulation breaks down at an earlier time than for ∆ = 0; we again see that the simulations work for t < ∼ N/ (2 v sw ) , (27)or t < ∼ l/v sw , (28)in the exact and light-cone methods respectively. D. Entanglement Entropy
From the predictions in [18] and the numerical work in [17], we can compare the difficulty of doing a similarcalculation doing time-dependent DMRG or time evolving bond decimation. The entanglement entropy is predictedto grow linearly in time. In [17], the prefactor of the numerical growth was determined in a quench of an XXZ chainfrom an Ising coupling with ∆ = ∆ to an Ising coupling with strength ∆ . The prefactor depended on ∆ and waslargest in the case of ∆ = ∞ , the case we have considered here where the initial condition is a Neel state. There,in a quench to ∆ = 0, the entanglement entropy (with logs taken to base 2) was observed to grow a little fasterthan 0 . v sw t . This implies that the minimum size of the bond dimension needed in a matrix product state is at least2 . v sw t , which requires an effort going as 2 . v sw t . In contrast, the numerical effort required to directly simulate achain of length N scales as 2 N . To get accurate results for a time of order t , requires N = 2 v sw t , or an effort 2 v sw t .Using the techniques in the present paper, the effort can be reduced to of order 2 v sw t .When comparing to [17], note that there is a difference in normalization of the Hamiltonian by a a factor of 4,but since the results in [17] are expressed in terms of the spin-wave velocity multiplied by time, the time axis is thesame for ∆ = 0. For ∆ >
0, the time axis in [17] differs from the time axis here, since [17] multiplies the time by therelevant spin-wave velocity which is greater than unity. For ∆ = 1, the spin-wave velocity using our normalization isequal to π/ >
1, and for ∆ = 0 . √ /
4, so the times in the present paper shouldbe multiplied by such a factor greater than unity when comparing to [17].In [17], other initial conditions were considered, other than just the Neel state. These initial conditions were chosento be the ground state of various other XXZ Hamiltonians. As the initial ∆ was reduced, the entropy growth was1
Time -0.100.10.20.30.40.5 < S z > f o r ce n t r a l s p i n exact, N=20l=16l=18l=205 10-0.02-0.0100.010.02 FIG. 7: Time dependence of h S z i for the central spin as a function of time for ∆ = 1. Exact curve is N = 20 (black). Light-conequantum circuit curves are l = 16 (red), l = 18 (green), and l = 20 (blue). found to still be linear in time, but with a smaller prefactor. Such simulations could still be carried out with ourmethod as follows: first, use DMRG as done in [17] to find the ground state on a long chain. Then, suppose we areinterested in local observable on a region of length l . To find how these evolve a time t after a quench, locate someregion of length l + 2 vt in the chain.We then “observe” the state of the system outside this region to statistically sample different pure states withinthe region, and then evolve the pure states within the region. This is done as follows. The DMRG ground state canbe written in the form Ψ = X αβ A αβ | Ψ αL i ⊗ | Ψ αβM i ⊗ | Ψ βR i , (29)where | Ψ αL i are a set of orthonormal states on the chain to the left of the region, | Ψ βR i are a set of orthonormal statesto the right of the region, Ψ αβM are a set of states (normalized to unity but not necessarily orthogonal) on the givenregion of length l + 2 vt , and A αβ are a set of amplitudes. We choose an α and a β according to the probability P = | A αβ | , (30)and then evolve the state Ψ αβM using our present algorithm. This corresponds to observing the matrix product statein the basis | Ψ αL i ⊗ | Ψ βR i . Repeating this statistical sampling many times we obtain the desired quantities on the localregion. Note that on each iteration we statistically sample α, β as well as doing the sampling above. The statisticalsampling of the state outside the region may be justified using Lieb-Robinson bounds as before. Further, when amatrix product state is written in the canonical form, the bond variables naturally have the orthonormality propertyused above to do the statistical sampling, giving a state in the form (29).In some cases, if the entanglement entropy grows linearly in time but with a sufficiently small prefactor, it may beworth using the light-cone ideas above, but doing the initial evolution for a time t i using matrix product methodsinstead of quantum circuit methods, as follows. Suppose the matrix product methods require an effort t ft , for somenumber f , to simulate for a time t . Then, to compute an observable at time t f , we simulate a subchain of length2 v sw t f for a time t i using matrix product methods. We then statistically sample states on a smaller subchain oflength 2 v sw ( t f − t i ) and perform the simulation of that subchain for time t f − t i exactly. The total effort is then t O (2 ft i + 2 v sw ( t f − t i )) . (31)2 Time < S z > f o r ce n t r a l s p i n exact, N=20l=18l=200 5 1000.1 FIG. 8: Time dependence of h S z i for the central spin as a function of time for ∆ = 2. Exact curve is N = 20 (black). Light-conequantum circuit curves are l = 18 (red) and l = 20 (green). Choosing t i = t f / (1 + f ) to minimize the computational cost, we find that the cost scales as t O (2 f ′ t f ) , (32)with f ′ = 1 f − + (2 v sw ) − ) . (33) E. The Corner Transfer Quantum Circuit Method
In this subsection we introduce the corner transfer quantum circuit method. It is primarily of theoretical, ratherthan practical interest. For calculation of local quantities (such as the expectation value of the spin on a single site)the light-cone method above is less work. However, the corner transfer method does give an approximation to thefull wavefunction, and may be less work than variational matrix product methods in cases where the entanglemententropy grows rapidly.We now define the quantum circuit that approximates exp( − iHt ). We define a length scale l ′ << v LR t : the errorin our quantum circuit approximation to exp( − iHt/
2) will be exponentially small in l ′ , while the maximal velocityof information propagation for our quantum circuit will be v max = v LR + O ( l ′ /v LR t ) . (34)We construct the quantum circuit in two steps, first presenting a quantum circuit that approximates exp( − iHt/ U k that describe the time evolution under a time dependent Hamiltonian: each U k contains ata time t only the interaction terms which are contained within one of the triangles with a flattened top and jagged sidessurrounded by a dashed line shown in in Fig. 9(a). That is, we break the time t/ n = ⌈ v LR ( t/ ⌉ subintervalsof time at most 1 /v LR . We place the center of the triangles on sites kl , for k integer, where l = 2 l ′ + v LR t. (35)3 V1 2 3 4 5 6 7 8UUU1,01,11,2
A)B)C)
U U UV V
FIG. 9: (a) The dashed lines show the region of space-time used in defining operators U k , with space on the horizontal axisand time on the vertical axis. Only three such regions are shown, but the pattern repeats over the entire system. We also showthis as a quantum circuit writing U = U , , U , , U , where each operator U ,n computes the exponential of the Hamiltonianin a given time slice. Construction is shown for l ′ = 2 , n = 3. (b) Action of operators V k as discussed in text. (c) Iteratingmany rounds of the corner transfer quantum circuit. On the bottom row we label U and V ; on the row above, a ˜ U sits aboveeach V and a ˜ V above each U . Then we define U k,n for 0 ≤ n ≤ n − U k,n = exp[ − i i ≤ kl + l ′ / n X i = kl − l ′ / − n h i ( t/ n )] , (36)and define U k = U k, U k, U k, ...U k,n − . (37)4We then define the operators V k as follows. We define U Lk,n = exp[ − i kl X i = kl − l ′ / − n h i ( t/ n )] , (38) U Rk,n = exp[ − i kl + l ′ / n X i = kl h i ( t/ n )] . That is, U L,Rk,n represent evolution under the left or right half of the triangles in in Fig. (9(a). Then, we define V k = exp( − i ( k +1) l + l/ − X i =( k +1 / l − l/ h i t/ (cid:16) U Rk, U Rk, ...U Rk,n − (cid:17) † (cid:16) U Lk +1 , U Lk +1 , ...U Lk +1 ,n − (cid:17) † . (39)That is, each V k “undoes” the evolution in the shaded triangles as shown in Fig. 9(b), and then performs the fullevolution in the rectangle bordered by the solid line.We now approximate exp( − iHt/ ≈ Y k V k Y k U k . (40)Using Lieb-Robinson bounds, one can show that the error in this approximation is k exp( − iHt/ − Y k V k Y k U k k ≤ O ( t X i k h i k exp( −O ( l ′ )) . (41)In the second step of construction the corner transfer quantum circuit, we define another approximation toexp( − iHt/ U k,n = exp[ − i i ≤ ( k +1 / l + l ′ / n X i =( k +1 / l − l ′ / − n h i ( t/ n )] , (42)and define ˜ U k = ˜ U k, ˜ U k, ˜ U k, ... ˜ U k,n − . (43)Here, the centers of the upward facing triangles in Fig. 9(a) are shifted by l/
2, compared to (37). We then define ˜ V k in analogy to Eq. (39) with the centers again shifted by l/ − iHt/ ≈ Y k ˜ V k Y k ˜ U k . (44)In Fig. 9(c) we show multiple rounds of the quantum circuit. Except for the row of triangles on the top and bottomboundaries, the quantum circuit looks like diamond-shaped patches in space time. They are rotated 45-degrees fromthe patches in [7], justifying the name “corner transfer”; the 45-degree rotation is the key improvement over [7] leadingto the bound on information propagation (34). III. QUANTUM BELIEF PROPAGATION
In this section we apply quantum belief propagation to disordered systems. We begin with a brief review of quan-tum belief propagation, focusing on the computational effort required, and then discuss modifications for disorderedsystems. We applied this procedure to two different disordered systems: one a chain with no frustration where MonteCarlo is available for comparison[29], and the other a frustrated spin system with disorder[30].Quantum belief propagation[22] is a method for constructing a matrix product density operator for a thermalstate of a quantum system on a line or other loopless lattice. The algorithm depends on a parameter l and theapproximation to the thermal state has the form in one dimension of: ρ = (cid:16) O † N − l +1 ...O † O † (cid:17) ρ (cid:16) O O ...O N − l +1 (cid:17) , (45)5where the operators O i act on sites i, ..., i + l − ρ has the form ρ = ρ ...l ⊗ ⊗ ... ⊗ ρ ...l being a density matrix on sites 1 ...l . The implementation of the algorithm depends on tracing out sites,a process analogous to the observations discussed above. Suppose we wish to compute the partition function. If allof the O i are the same, as they would be in a system without disorder, then we can consider the completely positivemap ρ → tr ( O † i ρO i ) ⊗ , (46)where tr ( ... ) denotes a trace over the first site, and ρ is a density operator on an interval of length l . We start withthis map at ρ = ρ , and iterate it until we reach the end of the chain. We then compute the trace of the final ρ andthis is the of the partition function. A similar procedure can be applied to compute expectation values.The cost of the procedure scales exponentially in l , as it requires diagonalizing matrices of dimension 2 l [26]. Theprocedure is effective down to inverse temperatures β ∼ l /v LR . (47)The physical intuition is that the algorithm keeps quantum effects only up to a length scaling l . However, as we willsee, the algorithm is capable of keeping track of classical correlations on much longer length scales; this should be nosurprise since for classical Hamiltonians which are the sum of commuting operators, such as the Ising Hamiltonian H = P i S zi S zi +1 , quantum belief propagation exactly reproduces classical transfer matrix methods which can exhibitcorrelation lengths exponentially large in β .For a disordered system, the operators O i may vary as a function of i , but they depend only on the bonds on sites i...i + l −
1. Below, we consider a system in which the bonds can assume only discrete values; in the first systemthere are 2 · ( l − / different possible choices for the set of bonds on sites i...i + l − l − possible choices. We then use the following algorithm: first, we pre-compute the operator O for each possiblechoice of bonds. We then randomly generate a configuration of bonds, and iterate the map (46) above, choosing theappropriate O at each step. This requires an effort which scales linearly in system size.Interestingly, the algorithm seems to work better at low temperatures on disordered systems than on orderedsystems. The reason is probably the following: when deriving the algorithm in [22], we used Lieb-Robinson boundswith velocity v LR . However, just as in the non-equilibrium case above where the actual velocity v sw is less than v LR ,in these thermal systems the actual velocity may again be less than v LR . For systems with disorder, localizationeffects may further reduce the velocity, and even change the ballistic spreading of the wavepacket to a slower growth.Interestingly, this phenomenon is known by different names in condensed matter, where it is called many-bodylocalization[27], and quantum information theory, where it is called a Lieb-Robinson bound for a disordered system[28]. A. Results on Disordered Systems
The first disordered system we considered has the Hamiltonian H = N/ X i =1 J ~S i − · ~S i + N/ X i =1 J i S i · S i +1 , (48)where each spin is spin-1 / J > J i = J F < p and J i = J A > − p .This model has both ferromagnetic and anti-ferromagnetic couplings, and was proposed to model the compound( CH ) CHN H Cu ( Cl x Br − x ) , where the probability p = x [33, 34, 35, 36, 37]. The case considered is J = 1 and | J i | = 2. We performed simulations on this chain with l = 5 , , p = 0 . , . , . , . l = 5 , , ,
11 for p = 0 ,
1. For the random cases ( p = 0 . , . , . , .
8) we considered chains of 100000 sites and computed the uniformsusceptibility by the change in partition function in response to a weak applied field. In this way, the susceptibilitywas self-averaging. For the pure cases, we considered shorter chains, and computed the applied susceptibility bymeasuring partition functions as in [22]. The results are shown in Figs. 10,11. For small l , qualitatively wrong resultsare seen at low temperature, with the susceptibility diverging at low enough temperature. However, as l is increased,the accuracy extends to lower temperature. The difference between the curves for l = 7 and l =9 is small for T greater than roughly 1 /
7. In this region also, we agree well with Monte Carlo data.Data was taken over several different β , with a step of 0 .
25 in β . As in [22], we set the perturbation A (followingnotation of [22]) to equal (1 / h l − ,l − + h l − ,l ) rather than A = h l − ,l .The second model we considered is [30] a model with frustration and disorder, where quantum Monte Carlo resultsare not available. The study in [30] was motivated by experimental studies on CuGeO [31] and Cu Ge O − xH χ u , p = χ u , p = β χ u , p = . FIG. 10: Uniform susceptibility as a function of temperature for different p, l . Top: l = 5 (black), l = 7 (red), l = 9 (green).Middle: l = 5 (black), l = 7 (red), l = 9 (green), l = 11 (blue). Bottom: l = 5 (black), l = 7 (red), l = 9 (green), l = 11(blue). where second neighbor interactions may be important. We were interested in a case where second neighbor interactionswould be very important, so we considered the Hamiltonian H = N − X i =1 J i ~S i · ~S i +1 + N − X i =1 K i ~S i · S i +2 , (49)and in the pure case we considered J i = 1 , K = 1 /
2. This is a Majumdar-Ghosh chain with a dimerized groundstate. In the disordered case we chose J i = 0 . J i = 1 .
1, with probability 1 / J = 0 .
75 or J = 1 . K i = (1 / J i J i +1 , (50)which correlates the second neighbor interaction with the nearest neighbor interaction as described in [29].We studied l = 5 , , β for the pure case, computed from the second derivative of the partition function (probably a very slightly moreaccurate method is to take the first derivative of the energy as in [22]). A strong difference is seen between l = 5 and l = 7 above β ∼ .
25, but the l = 7 and l = 9 curves are almost identical. This indicates that by going to l = 9we have succeeded in converging the specific heat in l for β ≤ l = 5 and l = 7 ,
9, but only slight difference can be seen between l = 7 , β = 8. Again, the results seem to beconverged in l by going to l = 9 in the range of temperatures we consider. The disordered curves show again that l = 5 is too small, but for l = 7 , β = 10.A very slightly higher susceptibility is seen in the random case compared to the pure case. Finally, we consider thedimer susceptibility, defined by χ dimer = β h (cid:16) N/ − X i =1 ~S i · ~S i +1 − ~S i +1 · ~S i +2 (cid:17) i . (51)7 χ u , p = . χ u , p = . β χ u , p = . FIG. 11: Uniform susceptibility as a function of temperature for p = 0 . , . , . l = 5 (black), l = 7 (red), l = 9 (green). This shows a large difference between the pure and disordered cases. The pure cases again show agreement between l = 7 , χ dimer /β growing exponentially in β . The disordered cases show χ dimer /β saturating as a function of β .The saturation of the dimer-dimer correlation function in the disordered case is no surprise. The disorder locallybreaks the Z symmetry between different ordering patterns of the dimers, and is relevant for this one dimensionalsystem. The fact that the uniform susceptibility shows only very slight difference between pure and disordered systemsis more surprising; asymptotically, the uniform susceptibility should decay exponentially in β in the pure system andshould increase as β/ log ( β ) in the random singlet phase[30]. IV. DISCUSSION
The main result in this paper is the light-cone quantum circuit method. We have tested this method numericallyon a free system, with ∆ = 0, and on interacting systems with ∆ = 0. We have found decaying oscillations in theexpectation value of the spin. In future, this technique will be useful for studying non-integrable systems, to see ifthey relax to a thermal state[38].We can study two-dimensional systems by considering them as wide one-dimensional systems; this allows us todouble the number of spins, but only leads to a factor of √ l = 18, we find resultscomparable to a system of N = 35 or N = 37 sites. The computational cost to exactly evolve a given system is roughlycomparable to that required to do exact diagonalization using Lanczos methods on that system: Lanczos methods andexact evolution both require sparse matrix-vector multiplication, but the number of multiplications needed to reachconvergence for the time evolution may be larger than that needed to reach convergence for ground state properties.Thus, we expect that sizes around 35 sites, especially given the low symmetry of the present system, are around theupper limit for exact methods now, while we carried the light-cone quantum circuit method up to l = 22. Further, theasymptotic analysis of time requirements applies also to memory requirements: the memory requirements of an exactsolution on N sites scale as N N , while while the memory requirement of the light-cone quantum circuit method scaleas l l , so, regardless of what N can be obtained using an exact solution, it should be possible to obtain the same l , upto a difference of a couple sites, in the light-cone quantum circuit method. The main additional cost in the light-cone8 C l =5l =7l =90 2 4 6 8 100.10.2 χ u l =5, purel =5, random0 2 4 6 8 10 β χ d i m e r l =5, pure FIG. 12: Top: specific heat for the pure system. Middle: uniform susceptibility for the pure system ( l = 5 , , l = 5 , , l = 5 , , l = 5 , , quantum circuit method is the need to run many times to obtain statistical samples, but this is a problem which canbe parallelized. Acknowledgments—
I thank R. Melko for useful comments on implementing sparse matrix-vector multiplication.I thank the KITP for hospitality while this research was conducted. This research was supported in part by theNational Science Foundation under Grant No. PHY05-51164. This work supported by U. S. DOE Contract No.DE-AC52-06NA25396. [1] U. Schollwoeck, Rev. Mod. Phys. , 259 (2005).[2] F. Verstraete, J. J. Garc´ıa-Ripoll, J. I. Cirac, Phys. Rev. Lett. , 207204 (2004); M. Zwolak and G. Vidal, Phys. Rev.Lett. , 207205 (2004).[3] G. Vidal, Phys. Rev. Lett. , 147902 (2003); G. Vidal, Phys. Rev. Lett. , 040502 (2004).[4] M. B. Hastings, J. Stat. Mech., P08024 (2007).[5] M. B. Hastings, Phys. Rev. B 73, 085115 (2006).[6] M.M. Wolf, F. Verstraete, M.B. Hastings, and J.I. Cirac, arXiv:0704.3906.[7] T. J. Osborne, Phys. Rev. Lett. 97, 157202 (2006).[8] E. H. Lieb and D. W. Robinson, Commun. Math. Phys. , 251 (1972).[9] M. B. Hastings and T. Koma, Commun. Math. Phys. , 781 (2006).[10] B. Nachtergaele and R. Sims, Commun. Math. Phys. , 119 (2006).[11] M. B. Hastings, Phys. Rev. B , 104431 (2004).[12] G. Vidal, J. I. Latorre, E. Rico, and A. Kitaev, Phys. Rev. Lett. , 227902 (2003).[13] V. Eisler and I. Peschel, J. Stat. Mech. P06005 (2007).[14] P. Calabrese and J. Cardy, J. Stat. Mech, P10004 (2007).[15] D. Gobert, C. Kollath, U. Schollwoeck, G. Schuetz, Phys. Rev. E , 036102 (2005).[16] S. Bravyi, M. B. Hastings, and F. Verstraete, Phys. Rev. Lett. , 050401 (2006); J. Eisert and T. J. Osborne, Phys. Rev.Lett. , 150404 (2006).[17] G. De Chiara, S. Montangero, P. Calabrese, R. Fazio, J. Stat. Mech.P03001 (2006).[18] P. Calabrese and J. Cardy, J. Stat. Mech. P04010 (2005).[19] S. Ostlund and S. Rommer, Phys. Rev. Lett , 3537 (1995); M. Fannes, B. Nachtergaele, and R. F. Werner, Comm.Math. Phys. , 443 (1992); F. Verstraete, D. Porras, and J. I. Cirac, Phys. Rev. Lett. , 227205 (2004). [20] To understand why 2 N/ states suffice, consider how much entanglement is possible across each bond. If we divide thechain of N sites into a subchain of n sites and another subchain of N − n sites, with n ≤ N − n , then there are only 2 n states available on the subchain and so it is possible to perform a Schmidt decomposition and find a matrix product staterepresentation with 2 n states across that bond. The worst case is midway across the chain, when n = N/
2. If we wish towrite matrix product operators instead of states, then these estimates are increased to 4 N/ .[21] M. Cramer, C.M. Dawson, J. Eisert, T.J. Osborne, arXiv:0703314.[22] M. B. Hastings, Phys. Rev. B Rapids 76, 201102 (2007).[23] X. Wang and T. Xiang, Phys. Rev B , 5061 (1997).[24] V. E. Korepin, N. M. Boloiubov, and A. G. Izergin, Quantum Inverse Scattering Method and Correlation Functions ,Cambridge University Press (1993).[25] The Krylov space methods may lead to some improvement in the time required compared to the series methods. However,since we are considering a situation in which the entropy increases rapidly, and hence the state exp( iH L t )Ψ L has very littleoverlap with Ψ L , it seems likely that it would be necessary to take a fairly large k to obtain accurate results for large t .Therefore, we are not sure how much advantage would actually be obtained by using Krylov methods here.[26] The quantum belief propagation method, as currently implemented, requires manipulating operators, rather than purestates, which increases the computational effort.[27] D. M. Basko, I. L. Aleiner, and B. L. Altshuler, Ann. Phys. , 1126 (2006); D. M. Basko, I. L. Aleiner, and B. L.Altshuler, cond-mat/0602510.[28] C. K. Burrell, T. J. Osborne, Phys. Rev. Lett. 99, 167201 (2007).[29] P. Zhang, Zh. Xu, H. Ying, J. Dai, cond-mat/0502543, conference on ”Non-Perturbative Quantum Field Theory: Latticeand Beyond”, Guangzhou, P.R. China, Nov. 2004.[30] E. Yusuf and Kun Yang, Phys. Rev. B 68, 024425 (2003).[31] J. E. Lorenzo et. al. Phys. Rev. B , 1278 (1994).[32] M. Hase, K. Ozawa, and N. Shinya, Phys. Rev. B , 214421 (2003).[33] H. Manaka, I. Yamada, and H. A. Katori, Phys. Rev. B , 104408 (2001).[34] H. Manaka, I. Yamada, and K. Yamaguchi, J. Phys. Soc. Jpn. , 564 (1997).[35] H. Manaka and I. Yamada, J. Phys. Soc. Jpn. , 1908 (1997).[36] H. Hida, J. Phys. Soc. Jpn. , 688.[37] T. Nakamura, J. Phys. Soc. Jpn.72