Out of equilibrium dynamics with Matrix Product States

Theoretical understanding of strongly correlated systems in one spatial dimension (1D) has been greatly advanced by the density-matrix renormalization group (DMRG) algorithm, which is a variational approach using a class of entanglement-restricted states called Matrix Product States (MPSs). However, DRMG suffers from inherent accuracy restrictions when multiple states are involved due to multi-state targeting and also the approximate representation of the Hamiltonian and other operators. By formulating the variational approach of DMRG explicitly for MPSs one can avoid errors inherent in the multi-state targeting approach. Furthermore, by using the Matrix Product Operator (MPO) formalism, one can exactly represent the Hamiltonian and other operators. The MPO approach allows 1D Hamiltonians to be templated using a small set of finite state automaton rules without reference to the particular microscopic degrees of freedom. We present two algorithms which take advantage of these properties: eMPS to find excited states of 1D Hamiltonians and tMPS for the time evolution of a generic time-dependent 1D Hamiltonian. We properly account for time-ordering of the propagator such that the error does not depend on the rate of change of the Hamiltonian. Our algorithms use only the MPO form of the Hamiltonian, and so are applicable to microscopic degrees of freedom of any variety, and do not require Hamiltonian-specialized implementation. We benchmark our algorithms with a case study of the Ising model, where the critical point is located using entanglement measures. We then study the dynamics of this model under a time-dependent quench of the transverse field through the critical point. Finally, we present studies of a dipolar, or long-range Ising model, again using entanglement measures to find the critical point and study the dynamics of a time-dependent quench through the critical point.


Introduction
The great success of experimental ultracold atomic physics has made the study of strongly correlated one-dimensional (1D) quantum systems a major avenue of current physics research. Examples of novel condensed matter physics realized with 1D atomic systems include the role of integrability in thermalization [1] and static [2] and dynamic [3] quantum simulators of Hubbard models. Furthermore, as ultracold molecules approach quantum degeneracy [4,5,6], lattice models with complex internal degrees of freedom and long-range interactions become relevant [7,8,9]. As more and more complex models become amenable to study, the need for numerical methods which can adapt to different degrees of freedom, different Hamiltonians, and different dynamical processes thus becomes essential.
In addition to practical interest in understanding and benchmarking atomic and molecular quantum simulators, the ability to simulate the dynamics of 1D systems also provides insight into fundamental questions such as the universality of dynamics approach quantum critical points and the effects of integrability on the thermalization process [10]. The natural setting for studying dynamics near critical points is a quantum quench where one of the parameters of the Hamiltonian is driven through a quantum critical point following a time dependent protocol, for example with θ (t) the step function. Such quenches pose a special difficulty for numerical studies as by definition they involve evolution with a time-dependent Hamiltonian which does not commute with itself at different times. The propagator is then generally a time-ordered exponential whose precise form may be difficult to ascertain. Standard methods such as the Suzuki-Trotter expansion which ignore the time dependence of the Hamiltonian require [11] that the infinitesimal time step used be much less than the fluctuation time-scale of H (t) to be valid, δt ≪ |∂H/∂t| −1 . This can cause simulations with rapid quench rates to become numerically very costly, and invalidates the approach altogether for non-analytic time dependence. Currently the only unbiased method available for the dynamics of quantum systems is exact diagonalization (ED). By unbiased, we refer to the fact that the other methods available for dynamics are generally variational, and so have a bias towards a particular ansatz. ED is limited in an essential way by the exponential growth of the size of the Hilbert space with the physical size of the system. The current state of the art is ∼ 40 spins for spin-1/2 models and 20 sites for a fermionic Hubbard model at half filling. These sizes are often too small for accurate finite-size scaling. An extremely powerful method for the low-energy properties of 1D systems is White's Density Matrix Renormalization Group (DMRG) algorithm, which uses a variational ansatz based on a class of states known as Matrix Product States (MPSs). MPSs will be reviewed in Sec. 2.1. DMRG uses an iterative procedure to develop a set of reduced bases that the full many-body problem is projected into, and then variationally minimizes an energy functional in this reduced space, enlarging it if necessary. DMRG uses an implicit MPS representation, which is to say that the state is not stored explicitly. This also means that the Hamiltonian and other operators in the calculation are stored in an approximate way, as they are represented within the reduced basis describing the variational state. This does not cause problems in practice when a single state is sought using the DMRG process. In fact, one can show that the algorithm to variationally find the ground state is identical when phrased in the implicit formulation of DMRG and when using an explicit MPS representation for the variational state, other than the representation of the Hamiltonian [12]. However, because of the exact representation of operators independent of the state, MPSs can put rigorous bounds on distances from exact quantum states such as eigenstates by considering quantities like the variance ψ| Ĥ − E 2 |ψ withĤ the Hamiltonian operator and E the energy expectation of the MPS |ψ . In contrast, DMRG can only return the distance of the variational state from the approximate operatorĤ in the given variational basis, and is unable to determine how well the given variational basis represents the true operator. A particularly clear indication of the failures this can cause is given in Ref. [13] where time evolution of a particular initial state in DMRG fails because the Hamiltonian has no nonzero matrix elements in the initial DMRG basis.
The situation becomes much different when multiple states are sought using the DMRG procedure. In this case the reduced density matrix used to determine the optimal reduced bases for the algorithm is a convex sum of the reduced density matrices for the desired states. This is called multi-state targeting. In contrast, an explicit MPS representation stores each of the desired states separately as an MPS. In multi-state targeting, none of the states can be represented with the same accuracy available if DMRG targeted that state alone. The MPS representation also deals automatically with the fact that each state has its own optimal bases for representation, whereas in DMRG these bases are all tied together by the multi-state targeting. In this work we present two algorithms which take advantage of MPSs' ability to deal with multiple states, eMPS to find excited states of 1D Hamiltonians and tMPS to simulate the dynamics of a generic time-dependent Hamiltonian. In the first algorithm a projector orthogonal to a set of lower-lying eigenstates is constructed from their MPS representations and used to orthogonalize a variational state against this set. In the second algorithm Krylov vectors in a Lanczos approximation to the matrix exponential are stored separately as MPSs and combined in an optimal way only at the end of the calculation. While Krylov-based MPS approaches have been used [13,66] to study time-dependent processes, the errors in these approaches were set by time derivative of the Hamiltonian. In contrast, by taking explicit account of the time ordering of the propagator, the errors in our approach are set only by commutators of the Hamiltonian at different times, and hence allow for larger time steps. Because of the explicit MPS representation, we are able to put bounds on the errors in each step of the calculations.
Finally, MPSs have a natural operator-valued extension known as Matrix Product Operators (MPOs) which allow for the exact representation of all operators used in the calculation. We present a general framework for constructing MPOs from a set of rules which is independent of the nature of the microscopic degrees of freedom. This allows for the templating of 1D Hamiltonians for general purpose software. In addition, the ability to perform arithmetic operations on MPOs exactly enables us to perform time-evolution using our tMPS algorithm with knowledge only of the MPO form of the Hamiltonian and the time-dependent functional form of its parameters. To emphasize the general nature of our algorithms, we include a generic simulation protocol for the out-of-equilibrium dynamics of strongly correlated 1D systems using the algorithms presented in this paper.
The remainder of this paper is organized as follows. In Sec. 2 we review the theory of MPSs, MPOs, and their canonical forms. In addition to providing a canonical form for operators within the matrix product formalism, we define finite state automaton rules for MPOs and demonstrate how 1D Hamiltonians can be constructed from a small set of such rules. In Sec. 3 we review the algorithm for finding ground states using MPSs as variational ansätze, and in Sec. 4 we present the eMPS algorithm which extends the ground state search to general excited states. Sec. 5 discusses how to extract observable quantities from MPSs. In Sec. 6 we discuss methods for time evolution with MPS. In particular, we provide the tMPS algorithm to time evolve an MPS using only the MPO representation of a Hamiltonian and the functional form of its time-dependent parameters. We contrast our approach with other Krylov subspace approaches and identify the possible sources of error. In Sec. 7 we present two case studies. The first is of the Ising model in a transverse field, where we study both the statics and the dynamics of a linear quench of the transverse field through the quantum critical point. The second is of a dipolar, or long-range Ising model in a transverse field, where we also determine the critical point from the statics and study a linear quench of the transverse field. Finally, in Sec. 8, we conclude. Details concerning numerically exact solutions for the Ising model which are used to benchmark our algorithms are given in Appendix A.

Matrix Product States
The Hilbert space of a quantum mechanical many-body system is exponentially large in the physical size of the system, for example the number of unit cells in a lattice or the number of particles. Stated another way, a state picked at random from the Hilbert space of a quantum many-body system will have entanglement (as quantified by the Schmidt measure [14]) which grows exponentially with the system size. In contrast to this random state, it has been shown that the class of states which are physically relevant in the sense that they can be prepared from some reference state by generic time evolution in polynomial time [11] or are useful for quantum computation [15] is much smaller than the full Hilbert space. In 1D, the physically relevant class of states appears to be those which have entanglement which is either constant or polynomially growing as a function of system size. A convenient representation of states with entanglement restricted in this manner is known as matrix product states (MPSs) [16,17,12].
We consider our physical system to comprised of a 1D lattice of sites, where each site i is a ddimensional Hilbert space H i spanned by the vectors {|i }. We will refer to d as the local dimension, and take all sites to be identical for simplicity. We define an MPS on a lattice with L sites as where the object A [k]i k is a χ k × χ k+1 matrix (with χ 1 = χ L+1 ) and Tr denotes the matrix trace. We will refer to the maximum linear dimension of any of the matrices A [k]i k , max k χ k , as the bond dimension of the MPS, and denote this quantity by χ. χ may be used as a measure of the entanglement of the state [12]. In this work we will focus on systems with open boundary conditions (OBC). MPS algorithms can also be devised for systems with periodic boundary conditions, as discussed in Refs. [18,19,7,20], but these algorithms have worse scaling and are generally less numerically stable than their OBC counterparts. For OBC, χ 1 = χ L+1 = 1, and arguments using the Schmidt decomposition demonstrate that χ k ≤ min(d k−1 , d L−k ) [21]. MPSs have been used for many years to represent exact ground states of parent Hamiltonians [22] which are formed from projectors onto local high-symmetry subspaces [23,24]. However, it was not until the great success of the density matrix renormalization group algorithm (DMRG) pioneered by White [25] that MPSs became valuable as variational ansätze in their own right [18]. Why are MPSs useful as variational ansätze? It has been shown [26] that the ground states of gapped 1D systems have bipartite entanglement which does not depend on the system size. Such states can be represented exactly as MPSs with a fixed bond dimension [27]. This is an example of an area law [28]; the entanglement between two disjoint subsystems depends only on the boundary of the two regions and not on their volume. For systems near a quantum critical point which is described by a conformal field theory (CFT), this area law is subject to weak logarithmic violations such that the entropy of entanglement between two subsystems of size L is given by the Calabrese-Cardy formula [29,30] where a is a constant and c is the central charge of the underlying CFT. Here ∼ denotes scaling equivalence in the bulk of an infinite system. In finite systems there are often oscillating boundary and finite size contributions [31,32]. Hence, the bond dimension of an MPS describing a conformally invariant critical system is given as χ L ∼ exp S L ∼ e a L c/6 , which grows only polynomially in the system size. Typical values of c range from 1/2 for the Ising model [33] and 1 for the Bose-Hubbard model [34] to 2 for more exotic phases like the gapless Mott insulator of the JK model [35]. We note that, strictly speaking, finding an MPS which approximates the ground state of an arbitrary 1D Hamiltonian to an accuracy which is an inverse polynomial in the system size is still NP-complete [36], but practical experience demonstrates that this method is extremely useful and robust for physical systems of interest. We adopt the following conventions for the representation of tensors: we use roman indices for indices which correspond to physical states and greek indices for indices which are summed over in the is represented as a contraction over rank three tensors with two rank-two boundary tensors.
A particularly useful means to visualize MPSs and manipulations with them is provided by tensor network diagrams like those shown in Fig. 1 [37]. Here a rank-k tensor is represented by a point with k lines extending from it. Each line represents one of the indices of the tensor. Whenever a line connects two points, that index is summed over, and disconnected lines represent free indices. Hence, an MPS can be represented as a chain of rank-3 tensors as in Fig. 1(c). Note that the first and last MPS tensors are rank two because we have assumed OBC and so χ 1 = χ L+1 = 1.
We note that the MPS definition Eq. (2) does not uniquely specify the tensors A. That is, we can insert an invertible matrix X and its inverse X −1 between any two MPS matrices without altering the physical content of the state: This is referred to as gauge freedom in the literature [21]. For OBC, we can specify the state uniquely ‡ by choosing a site k, which we call the orthogonality center of the MPS, and requiring that all sites i ‡ The state is unique up to possible degeneracies in the Schmidt decomposition. to the left and right of k satisfy the left and right gauge conditions, respectively. In these expressions, I is the appropriately dimensioned identity matrix. These conditions are shown in graphical notation in Fig. 2(a) and (b), respectively. § Graphically it is clear that the norm squared of the state is as shown in Fig. 2(c), and so this site carries all information about the norm of the state. This particular canonical form for an MPS is called mixed canonical form [12]. The mixed canonical form is crucial for improving the speed and numerical stability of variational algorithms with MPSs. We can impose the left gauge conditions via the following recursion: where Eq. (12) represents the singular value decomposition (SVD) ofÃ with Σ a diagonal matrix of singular values and U and V unitary matrices. Because U returned from the SVD is unitary, Eq. (8) is satisfied by construction. Similarly, the recursion for the right gauge conditions is § Here we also establish the graphical convention that downwards pointing lines correspond to Hermitian conjugates of tensors.
Note that any matrix decomposition ofÃ which returns a unitary matrix as part of the decomposition will suffice in place of the SVD. In particular, the QR decomposition [38] is particularly efficient when the rank ofÃ is not required.  To put a general state into mixed canonical form with orthogonality center k we begin at site 1 and iterate Eqs. (11)-(14) until we reach site k, then start at site L and iterate Eqs. (15)-(18) until we again reach site k.
Finally, we note that the set of all MPSs with a fixed bond dimension χ is not a vector space, as the sum of two MPSs with bond dimensions χ A and χ B has a bond dimension χ which is bounded by the sum of the two bond dimensions χ ≤ χ A + χ B . This can be seen from considering the sum of the two states |0 . . . 0 and |1 . . . 1 , with MPS representations The matrices A i and B i have a bond dimension of 1, as these are product states. Their sum is which has a bond dimension of 2.

Matrix Product Operators
The natural operator generalization of MPSs is a Matrix Product Operator (MPO), defined aŝ Here dimensional matrix, and we will again refer to the maximum value of χ O as the bond dimension of the operator. Note that this bond dimension χ O need not be the same as the bond dimension χ appearing in the MPS representation, Eq. (2). That an MPO takes MPSs to MPSs can be seen clearly from the graphical representation of Fig. 3. We also see from this representation that the bond dimension of the MPS representing the product of an MPO and an MPS is the product of the bond dimensions of the MPO and the MPS. Because the states |i 1 . . . i L are tensor products, we can also use the notation where each one of the objects That is to say, we can consider the matrices which appear in the matrix-product ansatz of an MPO to be operator-valued.
It is tempting to look for canonical forms for MPOs just as we did for MPSs, but the relevant norm for MPOs is the Frobenius norm Ô 1 ,Ô 2 = Tr(Ô † 1Ô 2 ) which scales exponentially in the local dimension with the number of lattice sites. Thus, for a typical many-body system with an exponentially large Hilbert space, the elements of the orthogonality center can differ in magnitude greatly, causing a catastrophic loss of precision during orthogonalization. However, most physically relevant MPOs such as one-dimensional Hamiltonians can be written down exactly in terms of an MPO canonical form which is analogous to an LU decomposition [39]. To discuss this canonical form, it is useful to recast an MPO as a finite state automaton (FSA) [40,41].
To recast an MPO as a FSA, we first enumerate all of the physical operators we use to define our local Hilbert space, O = {Ô α }, where α labels distinct operators. We will call this set our operator alphabet (OA). As examples, the OA for the Ising model would be Î ,σ x ,σ z withÎ the identity operator andσ ν the Pauli operator along the ν th Cartesian direction. The OA for the Bose-Hubbard model [42,43] would be Î ,b † ,b,n , whereb is a bosonic destruction operator andn =b †b the number operator. The particular matrix representation of the OA fixes the local basis states {|i }.
Using the OA, we now introduce a set of FSA rules {R p {Ô p 1 , . . . ,Ô pn }, {h p } , w p } which generate strings of the n operators {Ô p 1 , . . . ,Ô pn } ∈ O weighted by scalar w p and variables {h p }. We will call the variables {h p } Hamiltonian parameters and the scalars w p weights. Each rule consists of three operator-valued matrices in the case where the {h p } do not depend on position. The first matrix is the rightmost matrix in the MPO representation, W [L] , and represents the initial configuration. Next, we have the MPO matrix of the bulk of the chain W [j] , 2 ≤ j ≤ L − 1, which takes an input vector of operators on k sites and produces a vector of operators on k + 1 sites according to some deterministic pattern. Finally, we have the leftmost matrix in the MPO representation, W [1] , which extends the operators according to the patten of W [j] and returns a 1 × 1 operator-valued matrix. The trace of this matrix as in Eq (23) is the desired Hamiltonian term. The generalization to position-dependent Hamiltonian parameters requires L − 2 matrices in place of the bulk matrix W [j] , 2 ≤ j ≤ L − 1, but the only modification is that h p becomes h p (j). As a concrete example, consider the Ising model The Hamiltonian consists of two rules. The first is a site rule R site (σ x , h, −1) which generates the string −h iσ x i . The three matrices which provide this rule are As can be verified, and so this rule produces the desired operator. Similarly, there is a bond rule R bond ({σ z ,σ z }, J, −1) given by which produces −J i,j σ z iσ z j , with i, j denoting a sum over nearest neighbors i and j. The full Hamiltonian is given by the direct sum of the matrices. Collecting rows of the direct sum which are exactly parallel, we have the MPO representation of the full operator This construction can be readily extended to general sums of nearest-neighbor and on-site interactions with Hamiltonian by summing up the individual rules: The bond dimension of the MPO representation of this Hamiltonian is 2 + n B . MPO representations are not restricted to nearest-neighbor terms. Exponentially decaying terms of the form ¶ can also be accommodated with the rule R exp ({Ô 1 ,Ô 2 }, {J, λ}, −1) We can also produce a general monotonically decaying term 35) ¶ Note in this form that the nearest neighbor coupling is J and only longer ranged couplings contain λ.
by approximating the term as a sum of exponentials to a desired distance r cutoff . That is, we minimize the functional with respect to the exponential weights a i and decay parameters b i , where a and b are the elements a i and b i , respectively, arranged as vectors and n exp is a convergence parameter controlling the number of exponentials used in the expansion [44,45,46]. Note that the fit is only guaranteed to be accurate to r cutoff while the term has infinite range. This does not cause difficulties in practice for decaying functions, as the resulting fit is also guaranteed to be monotonically decaying. The decaying function rule is then just a sum of these n exp exponential rules, . In contrast to the other rules presented above this rule is not exact, but the number of exponentials can be increased to any desired accuracy + -a few tens of exponentials suffice for a 1/r 3 interaction on a lattice of 50 sites [46]. While the bond dimension increases linearly with the number of exponentials, the number of nonzero terms in the MPO also grows only linearly with the number of exponentials, and so operations with the MPO scale well as n exp increases. However, the eigenstate of a Hamiltonian with larger n exp may be more highly entangled than with smaller n exp , leading to an increase in χ and longer runtimes. That is to say, the relationship between the bond dimension of an MPO and the bond dimension of an MPS representing an eigenstate of the MPO is difficult to predict.
While the pure functional interaction Eq. (35) is appealing from a theoretical point of view, in practical applications infinite range interactions do not represent a consistent level of approximation. That is to say, at some distance interactions are screened, where the screening length is set by energetics or lifetime constraints in the case of ultracold atoms or molecules [47]. Hence, we provide the finiteranged rule R FiniteFunction ({Ô 1 ,Ô 2 }, {h, f, r cutoff }, w) which generates the Hamiltonian and is given by For a vector f which contains zero elements a small modification must be made, but the bond dimension remains the same. A small set of rules like those presented above allow us to template operators for MPS simulations. That is, given a method to produce a Hamiltonian from a set of basic rules, we can generate a wide + By accuracy we mean that the functional Eq. (36) is smaller than a given tolerance. variety of different Hamiltonians by specifying different OAs, weights, and Hamiltonian parameters as inputs to a program. The resulting MPOs are lower triangular and usually very sparse, and so sparse matrix structures can be used for efficiency. Beyond the flexibility that MPOs provide to MPS algorithms, they also allow for arithmetic operations such as addition and multiplication to be performed exactly, albeit at the expense of a growing bond dimension. The sum of two operators expressed as MPOs has MPO matrices which are the direct sums of the constituent matrices, and the product of two MPOs has MPO matrices which are the direct product of the constituent MPO matrices. Both of these operations preserve the sparse lower triangular structure of the MPOs.

Variational Ground State Search
We now turn to using MPSs as variational ansätze for the eigenstates of a Hamiltonian expressed as an MPO. The ground state is found by minimizing the functional with respect to the parameters of |ψ , where E is a Lagrange multiplier enforcing normalization. The general minimization of this functional is an NP-hard problem, so we instead adopt a local search heuristic that has proven to work well in practice. Let us consider fixing all parameters in the MPS except for a contiguous block of s MPS tensors A [j] . . . A [j+s −1] . The single-site (s = 1) and two-site (s = 2) algorithms are the most commonly used variants. We then find the extremum as which corresponds to the diagrammatic equation shown in Fig. 4. Here, the partial derivative with respect to a tensor is defined to be a tensor whose elements are the partial derivatives with respect to the elements of the tensor. If we assume that the block of tensors to be optimized contains the orthogonality center then the rightmost diagram reduces to the block of tensors being optimized and the leftmost diagram is the action of the effective Hamiltonian on this same block. Thus, minimization consists of finding the eigenvector corresponding to the smallest eigenvalue of the effective Hamiltonian eigenvalue problem where the effective Hamiltonian iŝ L and R are the partial overlaps of the Hamiltonian MPO with the state as in Fig. 4, and W [j] is the MPO tensor at site j of the Hamiltonian.Ĥ [j] eff represents the Hamiltonian for the variational degrees of freedom in the block to be optimized with the rest of the state held fixed. This justifies our use of E as the eigenvalue, as E obtained from the solution of this equation is the current best estimate for the energy. We can view Eq. (42) as a linear eigenvalue problem by combining indices as The linear dimension of this matrix representation ofĤ [j] eff is χ j d s χ j+s−1 , and so a solution of this problem with dense methods would require O(χ 3 j d 3s χ 3 j+s−1 ) basic operations, leading to a very slow algorithm Figure 4. Variational ground state search in diagrammatic notation. Here we display the single site (s = 1) effective Hamiltonian eigenvalue problem for simplicity. The contractions between the block of tensors to be optimized, the MPO, and the rest of the diagram have been omitted to accentuate the structure of the effective Hamiltonian. The rightmost equality follows from assuming that the variational site is the orthogonality center.
of order O(χ 6 ). In contrast, by taking advantage of the separable form of the effective Hamiltonian Eq. (43) multiplication of our block of tensors by the effective Hamiltonian can be done in O(χ 3 ) time [12]. Thus, sparse eigensolvers such as the Lanczos [38] or Davidson [48] algorithms, which require only matrix-vector multiplies, should be employed to solve this eigenvalue problem. The general algorithm for ground state search is thus as follows. We begin with an initial state with orthogonality center at site k. We choose a block of tensors containing k and optimize them by solving the effective Hamiltonian eigenvalue problem. We then shift the orthogonality center and the block of tensors one site to the right and again optimize. We continue shifting to the right until we reach the right boundary. We then reverse direction, shifting the orthogonality center and the block of tensors to be optimized to the left and solving the effective Hamiltonian eigenvalue problem until we reach the left boundary, at which point we reverse direction again. A single iteration of this optimization cycle which affects each tensor twice is called a sweep, and sweeping is continued until convergence. In addition, using the MPO form of the Hamiltonian, it is possible to develop a caching algorithm for the overlaps L and R such that the solution of this problem requires O (L) scaling in the number of lattice sites [41]. * To do so, we begin the iteration with a guess for the ground state |ψ assumed to have orthogonality center k. We then use the left recursion to generate the L overlaps up to k and the right recursion to generate the R overlaps down to k + s. Here the square braces indicate the order in which the contraction should be performed to achieve ideal scaling. Once the eigenvalue problem has been solved and the orthogonality center shifted, we use the recurrence Eq. (46) to update the overlaps when we * This scaling does not account for possible L dependence of the bond dimension χ such as exists for conformal critical points.
are sweeping to the right and the recurrence Eq. (48) to update the overlaps when we are sweeping to the left. Convergence is achieved when the variance with E the energy eigenvalue, drops below a user-specified tolerance ǫ. Given the MPO form of the Hamiltonian, the variance operator∆ ≡Ĥ 2 − E 2 can be constructed by constructing an MPO whose matricesW [i] consist of the direct product of the corresponding matrices fromĤ, , and then subtracting −ÎE 2 /L from the lower leftmost element of eachW [i] , whereÎ is the identity operator. This representation is exact, in contrast to DMRG-based approaches where the basis of the Hamiltonian is tied together with the basis of the state. The variance is a much better measure of convergence of the state than the so-called discarded weight which is used to measure convergence of the two-site DMRG algorithm. This is because it is a property of the actual MPS state and not of the eigenvalue. As a note of caution, the variance only guarantees that the state found is an eigenstate to the given tolerance, it does not specify that it is the ground state. This has not proven to cause problems in practice for non-disordered systems.
In summary, the complete algorithm for variational ground state search is: (i) Input: Input a HamiltonianĤ in MPO form, an initial guess |ψ for the ground state in MPS form, and a tolerance ǫ for the variance.
(ii) Initialization: Construct the LR overlaps using the recursions Eq. (46) and (48). (iv) Check convergence: Using the most recent estimate of the energy eigenvalueẼ from the last effective Hamiltonian solution, construct the variance operator∆ and find the variance. If ∆ < ǫ, exit, otherwise return to a.

Variational Excited State Search
We now turn to finding excited states. We find the n th excited state by minimizing the functional where {|φ k } are the n − 1 lower-lying eigenstates ofĤ and the {µ k } are Lagrange multipliers enforcing the orthogonality constraints ψ|φ k = 0. Again fixing a block of s tensors, the minimization of this Figure 5. Linear forms enforcing orthogonality. (a) The overlap ψ|φ k in diagrammatic notation. The thin red lines correspond to ψ| and the thick black lines to |φ k . (b) The linear form F [j]{k} in diagrammatic notation for the single-site case. As before, we leave the bottom tensor uncontracted to accentuate the definitions of the LR overlaps.
functional with respect to this block is given by the projected effective Hamiltonian eigenvalue problem whereP [j] is a projector into the space orthogonal to the {|φ k }. Given the states {|φ k } as MPSs, we construct these projectors as follows. The diagram corresponding to the overlap ψ|φ k is shown in Fig. 5(a), with the bold lines corresponding to |φ k and the thin lines to ψ|. This is a linear form in all of the MPS tensors of ψ|, and so the condition that |ψ be orthogonal to this state for the given block of tensors with all others held fixed may thus be stated as Here the linear form enforcing orthogonality F [j]{k} is shown diagrammatically in Fig. 5(b). We can construct the projectorP [j] aŝ where (N −1 ) kk ′ is the kk ′th element of the inverse of the Gram matrix This Gram matrix inverse is important to ensure idempotency of the projector. As before, direct construction of the projected effective Hamiltonian leads to an algorithm which scales poorly as O(χ 6 ). Hence, it is important to use sparse methods which require only the application ofP andĤ on some block of tensors A [j] . . . A [j+s −1] . Direct application ofP [j] as written requires O(χ 4 ) operations and also scales quadratically in the number of eigenstates desired N E due to the double sum in Eq. (54). To find a total of N E eigenstates by this method thus requires O(N 3 E χ 4 ) operations, which is unacceptably slow. A simple idea to reduce this scaling would be to find the eigenvectors of the inverse Gram matrix and reexpress the projectors F [j]{k} in terms of them, rendering the double sum a single sum. However, while the Gram matrix N is Hermitian and positive semidefinite it may also be badly conditioned and singular. A numerically stable alternative to this idea is to construct the Moore-Penrose pseudoinverse [38] of the Gram matrix where V is the matrix with the eigenvectors of N as columns and λ are the n p eigenvalues of N which are greater than n √ λ max ǫ, where n is the linear dimension of N , λ max its largest eigenvalue, and ǫ the machine precision. We use this pseudoinverse to transform to a new set of linear forms such thatP Often, the dimension of the set G [j]{µ} is much smaller than N E . The diagonalization of the Gram matrix requires O(N 3 E ) operations, independent of χ, and its construction and the construction of G in Eq. (57) both require O(χ 2 ) operations. The operation ofP [j] on the variational block of tensors is noŵ which is linear in N E and scales only as O(χ 2 ). Thus, the dominant scaling for typical parameters χ ≫ N E is still the O(χ 3 ) scaling of the effective Hamiltonian multiply, and the algorithm to find N E excited states scales as O(N E χ 3 ). A sweeping approach is used as in the ground state search algorithm, and the iteration is stopped when the variance drops below a user-specified tolerance. As before, the variance does not guarantee that the state found is the next lowest-lying eigenstate, but this does not usually cause problems in practice.
As with the LR overlaps used in the variational ground state search, one can also cache the overlaps LR used to construct the linear forms F using the recursions and where B [ℓ]{k} is the MPS tensor of |φ k at site ℓ. The linear forms are constructed using these overlaps as see Fig. 5(b). The variational ground state algorithm presented above is essentially equivalent to standard DMRG, aside from the calculation of the variance [12]. When finding excited states, however, DMRG-based approaches target multiple excited states in a single MPS, which causes the bond dimensions to grow and the quality of each individual eigenstate to degrade. Furthermore, as the ground state and all excited states are solved together in that approach, the sparse eigensolver must be able to converge interior eigenvalues, which is known to be troublesome [38,39]. We call the present algorithm, which is a sparse and numerically stable variant of that proposed in Ref. [49] for PBC, eMPS to accentuate the difference.
In our experience, there are two main limitations of eMPS. The first is that it is difficult to construct good variational guesses for the excited states in contrast to the ground state where the infinite size MPS algorithm [25,50] is applicable. Here, the usefulness of the variance becomes readily apparent, as the discarded weight can be 10 −12 or less while the variance is of order 10 −2 in early sweeps. The second is that that the area law considerations which make MPS algorithms so practical for ground states do not in general apply to bulk eigenstates, and so the bond dimension required to accurately represent a general eigenstate may be exponential in the system size, rendering eMPS inapplicable.
The ability to find excited states is useful in many contexts. It provides access to the dynamical gap for determining the location of second order quantum phase transitions [51] and Kibble-Zurek scalings, even when the gap is not between different symmetry sectors. It can help in understanding the structure of conformal field theories by providing access to the primary scaling fields [52]. Excited states yield the structure function and other dynamic response functions of low-lying excitations. Such response functions are of great use for comparing to experimental measurements. Finally, by considering more complex functionals such as ψ|(Ĥ − ǫ) 2 |ψ − λ ψ|ψ for minimization, one can determine level spacing statistics in a desired energy range for systems much larger than are amenable to exact diagonalization. Such studies are immensely useful in discussions of integrability and quantum chaos, as well as investigations of the thermalization hypothesis [53,54,55].

Calculation of observables
We now turn to how we can extract information from a state expressed as an MPS. We do so by the expectation values of Hermitian operators, or observables. We will demonstrate how to compute observables of three different types: local observables, two-point correlation functions, and general MPOs.
We define a local observable as an operator which acts only on the Hilbert space of a single site: ii ′ |i i ′ |. If this site corresponds to the orthogonality center of |ψ then the expectation value reduces to The overall scaling for fixed site index k is O (χ 2 d 2 ).

A two-point correlation function is an expectation value of the form Ô [q] †Ô[r]
where we take q < r without loss of generality. If the orthogonality center of the MPS, k, lies within the range q ≤ k ≤ r, then we can evaluate the expectation value using only the tensors in this range. The most efficient way to proceed is first to form the matrix recursively generate for ℓ = 1, . . . , r − q + 1, and then evaluate For fixed q and r, the algorithm scales as O (χ 3 d + χ 2 d 2 ).
To compute the expectation of a general many-body observableÔ expressed as an MPO we start from the right (left) boundary and follow the recursion Eq. (46) (Eq. (48)) all the way to the opposite boundary, at which point the remaining 1 × 1 × 1 tensor is the expectation value. The overall scaling is O ( Entanglement measures such as the bond entropy can be calculated from the singular values Σ of the MPS tensor A [j]i αβ as when this tensor is the orthogonality center. These singular values are computed automatically as part of the algorithm to shift the orthogonality center, see Sec. 2.1.

Time evolution with MPSs
We now turn our attention to a variational solution of the time-dependent Schrödinger equation using MPSs. The general strategy is to find some representation of the propagator over some time interval [t, t + δt],Û (t, t + δt), and variationally optimize the functional with respect to the MPS tensors of |ψ(t + δt) . Several complications arise in this case which were not present in the earlier algorithms. The first practical consideration is that the MPO form of the propagator may be difficult and very expensive to calculate. The second difficulty is more physical; the time-dependent state following a global quench of a Hamiltonian parameter has entanglement which generally grows linearly in time [56]. This causes the bond dimension χ to grow exponentially in time, and so there is some finite time where an MPS simulation will exhaust the available computational resources. However, many important questions regarding non-equilibrium dynamics can still be answered by considering moderately sized systems and short times. In addition, consideration of a situation in which the Hamiltonian changes only locally can greatly increase the accessible system sizes and simulation times [57,58,59]. The most common approach to time evolution for MPSs is to use the Suzuki-Trotter expansion or its higher order variants to construct a series of two-site propagators which can be constructed and applied easily. This is the basis of the equivalent [12] TEBD [60] and tDMRG [61,62] algorithms. Herê H n is the nearest-neighbor bond Hamiltonian acting on sites n and n + 1. This approach is no longer viable when the Hamiltonian has longer-ranged terms, and attempts to accommodate such longerranged terms often exhibit poor scaling [63,64] and require Hamiltonian-specialized implementation, resulting in inefficient, sometimes prohibitively inefficient code. Krylov-based time evolution, which will form the basis for our approach, has been considered in both DMRG [65] and MPS [13] variants for the time-independent case. We note that the latter approach has been used [13,66] to study time dependent systems, but this necessitated very small time steps set by the rate of change of the Hamiltonian in order to provide accurate results. Our approach generalizes the latter method to the time-dependent case where the error is independent of the rate of change of the Hamiltonian and demonstrates how the algorithm can be formulated entirely in terms of FSA rules for MPOs.

Commutator-free Magnus Expansions
The propagator of a general time-dependent Hamiltonian which does not commute with itself at different times is given as a time-ordered exponential whose most well-known form is the Dyson serieŝ This formulation of the propagator is not convenient numerically, as the Dyson series is an asymptotic series and so it can be difficult to determine an appropriate criteria for termination of the series. Furthermore, keeping only a finite number of terms in the Dyson series does not preserve the Lie group structure of the propagator; that is, the finite approximation is not unitary. An alternative approach which produces unitary approximations to the propagator was given by Magnus [67] who used the ansatzÛ (0, t) = exp −itΩ (t) (77) to define the Magnus serieŝ where the n th term is of order t n in the sense that its power series in t starts with t n . The termΩ n (t) involves n nested integrations over n − 1 nested commutators ofĤ (t) at different times. Explicitly, the first few terms are: While approximations obtained from truncating the series yield unitary propagators, these expressions are still formidable numerically, involving nested commutators and multidimensional integrals. The commutators pose a special difficulty for MPOs, as exact multiplication of MPOs involves multiplication of the bond dimensions of the MPOs and hence the algorithm scales exponentially in the number of terms kept in the series. Optimization algorithms which attempt to variationally shrink the bond dimension of an MPO sum or product such as those proposed in Ref. [46] may also be used, but these become numerically unstable for large systems, and when MPOs are subtracted as in commutators large cancellations can cause these algorithms to become stuck far from the variational optimum. Hence, rather than work directly with the Magnus series, Eq. (78), we start from ansätze of the formÛ where each one of theΩ i is a linear combination ofĤ at different times in the interval [t, t + δt], and require that our ansatz matches the Magnus expansion (equivalently, the full propagator) up to order δt N +1 . We will call such an ansatz a commutator-free Magnus expansion (CFME) [68,69]. This ansatz has a number of features which make it desirable for our purposes. It is exactly unitary and so the norm is conserved. Also, provided that we consider the case where only the Hamiltonian parameters change in time and the operators are time-independent, the sums of the Hamiltonian at different times can be represented exactly as an MPO using the rules of Sec. 2.2. Thus, the need for complex operations with MPOs vanishes. Finally, because the ansatz takes into account the time dependence of the Hamiltonian explicitly, the time step is not necessarily fixed by the rate of variation of the Hamiltonian, allowing for more coarse stepping in time with fixed error. Following Ref. [70], the procedure for generating an N th -order CFME is to expand the function H (t) in terms of (shifted) Legendre polynomials P n , The orthogonality properties of the Legendre polynomials allow the nested integration to be done exactly, leaving a series of nested commutators of theĤ n . Working in a Hall basis [71], this series of commutators is matched with the original Magnus expansion to yield the order conditions f i,n such that We note that these order conditions are independent ofĤ (t) by construction, and so are set by the choice of CFME alone. As we only require the result to be valid to order δt N +1 , the integration required for the coefficientsĤ n may be performed using Gauss-Legendre quadrature of order N/2 + 1. The end result of the analysis is that an N th order expansion with s exponentials may be written aŝ where x m and w m are the points and weights for Gauss-Legendre quadrature [72]. In this work we use a fourth order expansion with three exponentials (N = 4, s = 3). The optimal order conditions for this case, obtained in Ref. [70], are with f s−i+1,n = (−1) n+1 f i,n . Order conditions for higher order expansions may also be found in Ref. [70].
We now consider that our time-dependent Hamiltonian MPO is constructed from a set of FSA rules R p Ô p 1 , . . . ,Ô pn , {h p (t)} , w p in which the OA and the weights are chosen without loss of generality to be time-independent. In this case, the expansion Eq. (83) is applied individually to each Hamiltonian parameter h p (t), resulting in the parameters {h p n }. Now, because of the canonical decomposition of Sec. 2.2, the MPO forms ofΩ i from Eq. (86) at time t can be constructed exactly using the FSA rule set R p Ô p 1 , . . . ,Ô pn , N/2+1 m=1 g i,m h p (t + x m δt) , w p . We note that each one of these operators has the same bond dimensions as the original Hamiltonian, and the updates of operatorsΩ i at each time step can be done in O (Lχ O ) time, which is essentially negligible.
The fact that we can construct theΩ i from the same FSA rules as the Hamiltonian implies that our CFME ansatz is equivalent to evolving our system according to piecewise constant Hamiltonians of the same form but with differing Hamiltonian parameters. Additionally, as also occurs in high-order Suzuki-Trotter expansions, evolution backwards in time may occur. Finally, we note that even terms in the Hamiltonian whose parameters do not vary in time have their magnitude altered by Eq. (86), as m g i,m = 1 in general.

Krylov subspace propagation
Using the CFME Eq. (85) we never need to explicitly form an MPO representation of the propagator provided we can find an MPS representation of the exponential of an MPO applied to an MPS. We find such a representation from minimizing functionals of the form where, importantly,Ω has a known MPO representation. We do so by forming a Krylov subspace approximation to the exponential [73] in which the Krylov vectors are represented as MPSs. Specifically, we do so via the Lanczos algorithm for the matrix exponential, which can be stated as (iv) Finalize: i |v i Here T (j) is the symmetric tridiagonal matrix with the α i , 1 ≤ i ≤ j on the diagonal and β i , 1 ≤ i ≤ j − 1 on the superdiagonal. It is important to use a matrix exponentiation method which produces a unitary matrix to machine precision in order to not lose the Lie group structure. Because of the small linear dimensions of the matrix T (j) , exponentiation by direct diagonalization is practical. An a posteriori estimate for convergence of the Lanczos recursion is that 2β j−1 c (j) j < ǫ, where ǫ is the tolerance [73]. This can be compared with residual estimates in the ordinary Lanczos algorithm for finding eigenvalues. A rigorous bound on the approximation ||ψ krylov − |ψ | ≤ 12 exp − (ρδt) 2 16n eρδt 4n n can be established [74] when n ≥ ρδt/2 with n the number of Lanczos vectors and ρ = |E max − E min | the spectral width ofΩ. This estimate shows that for typical tolerances ǫ = 10 −6 to 10 −10 , 6 to 20 Lanczos vectors suffice.
As stated before, MPSs do not form a vector space and so the multiplication byΩ, the orthogonalization, and the final summation cannot be done exactly while keeping the bond dimension of our MPS fixed. However, just as with the eigenstate search, we can devise variational algorithms for these three operations which are iteratively performed until a desired tolerance is reached and use this tolerance to bound the bond dimension of our time-evolved MPS. We begin by briefly reviewing the standard algorithm [12] for finding the optimal MPS |φ representing a sum of MPSs N k=1 c k |ψ k to a given tolerance, as the other algorithms are similar but more complex. In this case we have a set of LR overlaps defined between our variational state φ| and the states |ψ k as in Eqs. (61) and (63). We now sweep through the lattice and make the replacement where the F [j]{k} are formed as in Eq. (64), see also Fig. 5(b). The orthogonality center of |φ is then shifted, the LR overlaps updated, and sweeping continued. Convergence can be monitored via with Re denoting the real part. Because we do not have to solve an eigenequation at each iteration, this algorithm is often much less costly than the iterative eigenstate search. Also, when we have that the coefficient vector c and all of the {|ψ k } have length 1, we can normalize the state |φ at the end of the calculation if required. The algorithm to variationally fit an MPS toΩ|ψ is similar. In this case we have a set of LR overlaps defined via the recursions where the MPS tensors of |φ are denoted by A and those of |ψ denoted by B. We now sweep through the lattice and make the replacement where the effective Hamiltonian is formed from the LR overlaps as in Eq. (43). Convergence can be monitored via Here, ψ|Ĥ 2 |ψ can be computed in a manner similar to the variance, and need only be computed once at the beginning of the calculation. We have also assumed that the block of tensors in Eq. (97) contains the orthogonality center. Again, this algorithm is often much less costly than the iterative eigenstate search.
We now turn to steps (iii)(d) and (iii)(e) of our algorithm.
Step (iii)(d) is usually stated for ordinary vector spaces as as α j = v j |r and β j−1 = v j−1 |r and so this is equivalent to classical Gram-Schmidt orthogonalization. Hence, we could implement step (iii)(d) by using the fitting algorithm Eq. (90) to find the MPS closest to |r − α j |v j − β j−1 |v j−1 . However, we have found that the following algorithm, which is closely related to eMPS, often converges more quickly and also is applicable to step (iii)(e). In our method we look for the optimal MPS |φ representing |ψ but also subject to the constraints that φ|ψ k = 0 for some set {|ψ k }. We start by copying the state |ψ to a variational ansatz |φ . We then construct overlaps between the state φ| and |ψ , which we call LR♯ and a set of overlaps between φ| and |ψ k , which we call LR. We then sweep through the lattice and make the replacement with B the MPS tensors of |ψ and A the MPS tensors of |φ . We then apply the projector into the space orthogonal to the |ψ k by constructing the set G [j]{µ} , µ = 1, . . . , p, as in Eq. (57) and performing for µ = 1, . . . , p. Using the fitting algorithm of Eq. (90) corresponds to replacing the LR overlaps, which are between the variational state φ| and the set {|φ k }, with a set of LR overlaps between the state ψ| and the set {|φ k }. Our algorithm, which amounts to fitting followed by modified Gram-Schmidt, uses information about the distance between the variational state and those to be orthogonalized against to determine operations, and hence often converges more quickly and is more stable. Convergence can be monitored by ensuring that φ|ψ k are orthogonal to a precision ǫ via If one requires additional truncation of the bond dimension, one can switch to the ordinary fitting algorithm Eq. (90) at this point, using a new variational state |ζ to fit to |φ . We now pause to consider the sources of error in the time-propagation routine. First, because the CFME expansion Eq. (85) is only of order δt N +1 , the error incurred in using this form of the propagator is ǫ CFME = ct final δt N , where t final is the final time reached. Thus, as the final time desired becomes longer, smaller time steps should be taken in order to keep the error fixed. The coefficient c can be determined by using this known scaling and decreasing the time step by a constant factor. Factors in the range 2 1/N to 3 1/N are practical. This strategy can also be used to devise adaptive time-step applications such as those used widely in ordinary differential equation solvers. Next, there is the error ǫ Krylov incurred in the Krylov subspace approximation to the exponential. As discussed above, this error can be minimized by adding more and more Lanczos vectors. This error also involves the time step δt, and so when determining the CFME expansion error constant c one should be careful that ǫ Krylov < ǫ CFME . Finally, there are errors resulting from the variational fitting of MPSs in steps (iii)(b), ♯ In this initialization all of the LR are Kronecker deltas provided that |ψ has an orthogonality center.
(iii)(d), and (iii)(e) of the Lanczos algorithm for the matrix exponential. These can be reduced by lowering the variational tolerances, but this is done typically at the expense of a larger bond dimension χ and hence a slower algorithm and more memory usage.

Simulation Protocol
We are now in a position to devise a complete, generic protocol for the time evolution of a 1D quantum system.
Here J is the coupling energy, h is a transverse magnetic field, and the {σ i } are the Pauli spin operators on site i. We choose this model because its dynamics are amenable to numerically exact study using the time-dependent Bogoliubov-de Gennes formalism (see Appendix A for a review) and so we were able to carefully check convergence of our results. The statics have all been converged to eight digits, and the dynamics at all times to at least four digits. Here we refer to convergence of local quantities such as the energy or density of defects. Nonlocal quantities such as the bond entropy will not have this same accuracy, but numerical tests show that the qualitative behavior is unaffected. We begin with a discussion of the statics. In Fig. 6(a) we demonstrate the gaps from the ground state to the two lowest eigenstates, computed via eMPS. The variances are smaller than the point size in this case. The upper (lower) curve corresponds to even (odd) parity, while the ground state has even parity, where parity is defined as simultaneous inversion of all spins P = iσ x i . Hence, the relevant gap for discussing the quantum phase transition is in fact the gap to the second excited state, shown in green. † † Both gaps close at criticality, and this can cause the first excited state returned by eMPS to be a mixture of these two states. This will not affect the energies so long as the variance remains small, but it can affect other observable properties of the states. There are two ways to remove this nuisance. The first, most complex, and most preferable is to use an MPS representation in which the state is explicitly Z 2 invariant [39]. The second is to add a field coupling to the parity −h p iσ x i to cause the relevant even-parity subspace to become lower in energy than the odd-parity subspace. This operator is an MPO with bond dimension 1, the MPO equivalent of a product state. The closing of the gap at the known critical point h = 1 is linear in 1/L, indicating a conformally invariant critical theory with dynamical critical exponent z = 1.
We venture to determine the central charge of the critical theory by fitting to the Calabrese-Cardy formula in two ways. In the first, we fit to the finite-size formula 103) † † The even parity gap for h > 1 is in fact twice the demonstrated gap in green, but the essential piece is the closing of the gap at criticality. for fixed L and variable i, and in the second we fix i at L/2 and fit S L/2 to this formula for various L. Near criticality, the presence of nonzero c indicates a curvature of S i , while in the gapped phases S i obeys an area law and is hence a constant apart from finite-size effects. The bond entropy in the bulk approaches the correct limiting values of log 2 as h → 0 and 0 as h → ∞. The first fit, shown in Fig. 6(a), provides us with a clear indicator of the critical region by the spike in the central charge c. However, the precise determination of c for a finite size system in this case is noisy, likely due to strong finite-size effects. Once we have narrowed down where the critical region is, the second fit, shown in Fig. 6(b), allows us to extract the anticipated value c = 1/2 much more precisely. If the same scaling analysis is attempted at a point which is not the critical point, the bond entropy saturates and c → 0 as L → ∞. We can understand this as a large but finite correlation length ξ. For L/2 < ξ, the system appears to be conformal and we see scaling of the bond entropy with L. For L > ξ the bond entropy saturates and this scaling breaks down, indicating that the given region is not critical. We note that in this analysis we have used no properties which are specific to this system e.g. correlation functions of an order parameter to extract the critical behavior.
We now turn to the dynamics. The Ising model has also been a subject of interest for dynamics as a testbed for the Kibble-Zurek hypothesis (KZH) that equilibrium properties determine nonequilibrium properties following a quench across a quantum critical point. This was studied in Refs. [75,76,77] using the quench protocol A useful quantity for determining how non-adiabatic a particular quench is in this case is the density of defects which is the density of Bogoliubov quasiparticles at zero magnetic field. In addition to the density of defects, universal scaling has also been predicted in the heat, or non-adiabatic part of the energy, where |ψ g.s. (t) represents the instantaneous ground state ofĤ(t). In addition to these quantities, which are amenable to Bogoliubov-de Gennes analysis, we also compute the time-dependent bond entropy.
Our results are shown in Fig. 7. We first discuss the heat, as shown in panels (a) and (d). Panel (a) displays the heat as a function of time, and demonstrates a sharp change in the behavior of the system as we pass through the critical point. This is especially true of the longest quenches. In panel (d) we investigate the heat as a function of the quench rate both at the time when h takes on its critical value, t c = 4τ /5, and at the final time when h = 0. The large difference indicates that non-adiabatic processes continue after we have passed from the critical region back into the gapped ferromagnetic region. Thus, the universal scaling of the heat may be difficult to determine if the critical point itself is not known sufficiently well. We now move on to the density of defects as shown in panels (b) and (e). In panel (b) we see that the density of defects at the final time decreases slowly to zero as τ → ∞; that is, when the quench becomes perfectly adiabatic. This is in accordance with the KZH prediction. Panel (e) demonstrates the large disparity between the density of defects at the critical time and the final time for all but the most rapid of quenches. Finally, in panels (c) and (f) we investigate the bond entropy. In panel (c) we see the bond entropy of the central bond as a function of time. As the quench becomes more adiabatic, the entropy increases more towards the ferromagnetic limiting value of log 2. However, for very slow quenches, the bond entropy reaches this value before the end of the quench and then begins to oscillate. In panel (f) we show the bond entropy as a function of the bond index at the critical time. Bulk curvature such as that seen at criticality in Fig. 6(d) is absent, indicating that the time-evolved state is still quite far from the conformal ground state.

Case Study: Dipolar Ising chain
In this section we go beyond exactly solvable models and study a dipolar Ising chain Such models are relevant to the study of ultracold molecules in optical lattices, where the dipole-dipole interaction falls away as 1/r 3 with r the distance between dipoles [78,7,8,9,47]. Here the cutoff |j − i| ≤ 6 represents a consistent order of approximation in going from a Hubbard-type model with dipolar interactions to a spin model. We stress that all results obtained in this section were obtained using the same implementation as the last section.
We first turn to the statics of this model, shown in Fig. 8. Many of the features are similar to those of the nearest-neighbor Ising model. The most important differences are that the critical region is shifted towards larger h with respect to the nearest-neighbor Ising model as seen in panels b) and c). This indicates increased stability of the ferromagnetic phase, in accordance with expectations. Using these points as a guide, we determine the critical point to be h c = 1.362 ± 0.01, as shown by the scaling in panel d). Additionally, the energy of the first even parity excited state deviates from pure z = 1 linear behavior near h = 0, indicating interactions between quasiparticles which were noninteracting in the nearest-neighbor Ising model.
We now turn to dynamics, following the same quenching protocol Eq. (104) as above. The results are shown in Fig. 9. We reiterate that the dynamics of this model cannot be handled by Bogoliubov-de Gennes methods, nor straightforwardly with standard tDMRG/TEBD. The density of defects no longer represents the density of quasiparticles at the critical point, but we compute it for comparison with the results of the nearest-neighbor Ising model. Because of the larger MPO bond dimensions and the more rapid growth of bond dimension for this model, we restrict our analysis to short times Jτ / ≤ 5, though an optimized implementation could reach longer times. The basic features are similar to the dynamics of the nearest-neighbor Ising model. One quantitative difference is that, because the critical point is reached at an earlier time than in the nearest-neighbor Ising model, oscillations in the central bond entropy occur for shorter quench times.

Conclusions
The power of matrix product state algorithms over DMRG-based algorithms is most readily apparent when multiple states are involved, as each state may be represented as a separate matrix product state in the former approach. Because matrix product states with a fixed bond dimension do not form a vector space, a set of matrix product states carries more information at smaller numerical cost than the same set represented as a multi-state targeted basis in DMRG. We have presented two algorithms, eMPS and tMPS, which use this property to find eigenstates and perform time evolution of strongly correlated 1D quantum systems. eMPS uses a set of eigenstates stored as separate matrix product states to define a projector into the space orthogonal to this set. We use this projector to explicitly orthogonalize a variational state against previously determined eigenvectors in order to find excited states. The explicit matrix product state representation allows us to store the excited states much more accurately than with standard DMRG, and allows allows us to ensure global orthogonality between the eigenstates to a desired tolerance. The variance, which is computed exactly using the matrix product operator representation of the Hamiltonian, gives strict error bars on the energies obtained with this procedure. tMPS avoids the need for an explicit representation of the propagator by using a commutatorfree Magnus expansion and then building successive Krylov subspace approximations to the matrix exponentials which appear in the expansion. Each vector in the Krylov subspace is stored as a separate matrix product state to maximize efficiency. Furthermore, the operatorsΩ i have exact representations as matrix product operators with the same bond dimension as the Hamiltonian. Our algorithm eliminates the need for Hamiltonian-specialized implementation of dynamics. Additionally, by carefully accounting for the time dependence of the Hamiltonian with a commutator-free Magnus expansion, the error in our algorithm depends only on commutators of the Hamiltonian with itself at different times and not on its derivatives. As with eMPS, the errors are rigorously accounted for by considering distance functionals with the variational state.
The matrix product operator forms of 1D Hamiltonians can be obtained using a small set of finite state automaton rules such as the site, bond, and finite function rules. Using matrix product operator arithmetic, we can add together the various terms in a Hamiltonian from these rules to form a complete canonical MPO representation. This representation allows for templating of Hamiltonians which depend only on the type of interactions and not on the microscopic constituents of the lattice model. Furthermore, given the time-dependent form of the Hamiltonian parameters, one can use the same template to form the operatorsΩ i which appear in tMPS at negligible numerical cost. We used our algorithms to study both the nearest-neighbor Ising model in a transverse field and a long-range Ising model in a transverse field. By the closing of the gap obtained with eMPS we determined that the critical points of these models were conformal, and so we used the Calabrese-Cardy formula for the bond entropy of conformal systems to locate the critical point and its associated central charge. The known result h c = 1 was verified for the nearest-neighbor case, and the critical point was shifted deeper into the paramagnetic region h c = 1.362 ± 0.01 for the long-range case, indicating a stabilization of the ferromagnetic phase. We used tMPS to study the dynamics of these models following a linear quench of the transverse field from the paramagnetic phase through the critical point into the ferromagnetic phase. We saw strong non-adiabatic effects in the heat as quenching continued into the ferromagnetic region, scaling of the density of defects consistent with the Kibble-Zurek hypothesis, and the oscillation of the bond entropy near its limiting ferromagnetic phase value for nearly adiabatic quenches. This work was supported by the Heidelberg Center for Quantum Dynamics, the Alexander von Humboldt Foundation and the National Science Foundation under Grant PHY-0903457. We also acknowledge the Golden Energy Computing Organization at the Colorado School of Mines for the use of resources acquired with financial assistance from the National Science Foundation and the National Renewable Energy Laboratories. We thank Erman Bekaroglu for discussions and a thorough reading of the manuscript.
Appendix A. Time evolution of exact solution of 1D transverse-field quantum Ising model for comparison with tMPS The solution of the statics of the transverse-field quantum Ising model is covered in standard texts [51]. However, for comparison with tMPS we require a description of the dynamics, and so we present the dynamical case here. To find the exact solution of the Ising model, we affect the Jordan-Wigner transformation where the fermionic operatorsĉ i satisfy the anticommutation relations {ĉ i ,ĉ † j } = δ ij , {ĉ i ,ĉ j } = {ĉ † i ,ĉ † j } = 0. This transforms the Ising model into the fermion Hamiltonian As this is a quadratic form in fermion operators, it may be diagonalized by a canonical (Bogoliubov) transformation [79,80] which provides the set of Bogoliubov-de Gennes equations ǫ k u k = Au k + Bv k (A.6) where u k are the elements of {u ik , i = 1, . . . , L} arranged as a vector and likewise for v k . The matrices A and B are real and tridiagonal, with the nonzero matrix elements A i,i = 2h, A i,j = −J, |i − j| = 1 and B i,i+1 = −B i+1,i = −J. The pairs (u ik , v ik ) with positive energy ǫ k , ǫ 1 ≤ ǫ 2 ≤ . . . ≤ ǫ L , and the normalization |u k | 2 + |v k | 2 = 1 define quasiparticle operatorŝ which bring the Hamiltonian into the diagonal form Corresponding to every such pair is another pair (ũ ik ,ṽ ik ) = (v ik , u ik ) withǫ k = −ǫ k which defines the conjugate operatorγ † k . Writing u k and v k together as a composite vector, the Bogoliubov-de Gennes equations take the form of a real symmetric eigenproblem of dimension 2L: Evolution under the fermion Hamiltonian Eq. (A.4) does not preserve the number of fermions N F but it does preserve the fermionic parity (−1) N F . Because the ground state is the Bogoliubov vacuum it contains no fermions, and so the first accessible excited state consists of two Bogoliubov excitations, one in each of the lowest two modes. The gap between the ground and first excited states is thus ǫ 1 + ǫ 2 .
We now consider the Heisenberg equations of motion for the fermi operators where u i (t) and v i (t) subject to the time-dependent Bogoliubov-de Gennes equations i d dt andγ k andγ † k diagonalize the Hamiltonian at the initial time. Equivalently, u and v define timedependent Bogoliubov operatorŝ such that the time-evolved state |ψ(t) is the Bogoliubov vacuum of this set, i.e.γ k (t)|ψ(t) = 0. To compare with the MPS simulations, we note that the energy at time t is Similarly, the density of defects is (v i,k (t) − u i,k (t)) u ⋆ i+1,k (t) + v ⋆ i+1,k (t) . (A.17) Time evolution thus reduces to the solution of a 2L×2L time-dependent matrix differential equation which we solve using a CFME as in Sec. 6.1. Because the dimensions of the system are much smaller than those of the associated MPS problem we are able to take very small time steps, and so the results obtained in the method may be considered to be numerically exact, see Fig.A1. Here we demonstrate that the error incurred in the energy as a function of time scales with the infinitesimal time step δt as δt 4 using our fourth-order CMFE. Hence, by decreasing δt, any desired degree of accuracy may be met.