Time-evolving block decimation
The time-evolving block decimation (TEBD) algorithm is a numerical scheme used to simulate one-dimensional quantum many-body systems, characterized by at most nearest-neighbour interactions. It is dubbed Time-evolving Block Decimation because it dynamically identifies the relevant low-dimensional Hilbert subspaces of an exponentially larger original Hilbert space. The algorithm, based on the Matrix Product States formalism, is highly efficient when the amount of entanglement in the system is limited, a requirement fulfilled by a large class of quantum many-body systems in one dimension. Introduction. There is nowadays a considerable interest in  the field of quantum theory for computational methods well-suited to the physics of many-body systems. Considering the inherent difficulties of simulating general quantum many-body systems, the exponential increase in parameters with the size of the system, and correspondingly, the high computational costs, one solution would be to look for numerical methods that deal with special cases, where one can profit from the physics of the system. The raw approach, by directly dealing with all the parameters used to fully characterize a quantum many-body system is seriously impeded by the lavishly exponential buildup with the system size  of the amount of variables needed for simulation, which leads, in the best cases, to unreasonably long computational times and extended use of memory. To get around this problem a number of various methods have been developed and put into practice in the course of time, one of the most successful ones being the quantum Monte Carlo method (QMC). Also the density matrix renormalization group (DMRG) method, next to QMC, is a very reliable method, with an expanding community of users and an increasing number of applications to physical systems. When the first quantum computer is plugged in and functioning , the perspectives for the field of computational physics will look rather promising, but until that day one has to restrict oneself to the mundane tools offered by classical computers. While experimental physicists are putting a lot of effort in trying to build the first quantum computer, theoretical physicists are searching, in the field of quantum information theory (QIT), for genuine quantum algorithms, appropriate for problems that would perform badly when trying to be solved on a classical computer, but pretty fast and successful on a quantum one. The search for such algorithms is still going, the best-known (and almost  the only ones found) being the Shor's algorithm, for factoring large numbers, and Grover's search algorithm. In the field of QIT one has to identify the primary resources necessary for genuine quantum computation. Such a resource may be responsible for the speedup gain in quantum versus classical, identifying them means also identifying systems that can be simulated in a reasonably efficient manner on a classical computer. Such a resource is quantum entanglement; hence, it is possible to establish a distinct lower bound for the entanglement needed for quantum computational speedups. Guifré Vidal, then at the Institute for Quantum  Information, Caltech, has recently proposed a scheme useful for simulating a certain category of quantum systems. He asserts that ""any quantum computation with pure states can be efficiently simulated with a classical computer provided the amount of entanglement involved is sufficiently restricted" This happens to be the case with generic Hamiltonians displaying local interactions, as for example, Hubbard-like Hamiltonians. The method exhibits a low-degree polynomial behavior in the increase of computational time with respect to the amount of entanglement present in the system. The algorithm is based on a scheme that exploits the fact that in these one -dimensional systems the eigenvalues of the reduced density matrix on a bipartite split of the system are exponentially decaying, thus allowing us to work in a re-sized space spanned by the eigenvectors corresponding to the eigenvalues we selected. One can also estimate the amount of computational resources required for the simulation of a quantum system on a classical computer, knowing how the entanglement contained in the system scales with the size of the system. The classically (and quantum, as well) feasible simulations are those that involve systems only lightly entangled—the strongly entangled ones being, on the other hand , good candidates only for genuine quantum computations. The numerical method is efficient in simulating real-time dynamics or calculations of ground states using imaginary-time evolution or isentropic interpolations between a target Hamiltonian and a Hamiltonian with an already-known ground state. The computational time scales linearly with the system size, hence many-particles systems in 1D can be investigated. A useful feature of the TEBD algorithm is that it can be reliably employed for time evolution simulations of time-dependent Hamiltonians, describing systems that can be realized with cold atoms in optical lattices, or in systems far from equilibrium  in quantum transport. From this point of view, TEBD had a certain ascendance over DMRG, a very powerful technique, but until recently not very well suited for simulating time-evolutions. With the Matrix Product States formalism being at the mathematical heart of DMRG, the TEBD scheme was adopted by the DMRG community, thus giving birth to the time dependent DMRG , t-DMRG for short. Around the same time, other groups have developed similar approaches in which quantum information plays a predominant role as, for example, in DMRG implementations for periodic boundary conditions , and for studying mixed-state dynamics in one -dimensional quantum lattice systems.<ref name="cond-mat/0406426"> </ref><ref name="cond-mat/0406440"></ref> Those last approaches actually provide a formalism that is more general than the original TEBD approach, as it also allows to deal with evolutions with matrix product operators; this enables the simulation of nontrivial non-infinitesimal evolutions as opposed to the TEBD case, and is a crucial ingredient to deal with higher-dimensional analogues of matrix product states. The decomposition of state. Introducing the decomposition of State. Consider a chain of "N" qubits, described by the function formula_1 This form, known as a Matrix product state, simplifies the  calculations greatly. To understand why, one can look at the Schmidt decomposition of a state, which uses singular value decomposition to express a state with limited entanglement more simply. The Schmidt decomposition. Consider the state of a bipartite system formula_2. Every such state formula_3 can be represented in an appropriately chosen basis as: where formula_5 are formed with vectors formula_6 that make an orthonormal basis in formula_7 and, correspondingly, vectors formula_8, which form an orthonormal basis in formula_9, with the coefficients formula_10 being real and positive, formula_11. This is called the Schmidt decomposition (SD) of a state. In general the  summation goes up to formula_12. The Schmidt rank of a bipartite split is given by the number of non-zero Schmidt coefficients. If the Schmidt rank is one, the split is characterized by a product state. The vectors of the SD are determined up to a phase and the eigenvalues and the Schmidt rank are unique. For example, the two-qubit state: has the following SD: with On the other hand, the state: is a product state: Building the decomposition of state. At this point we probably know enough to try to see how we explicitly build the decomposition (let 's call it "D"). Consider the bipartite splitting formula_18. The SD has the coefficients formula_19. By expanding the formula_20's in the local basis, one can write: The process can be decomposed in three steps, iterated for each bond (and, correspondingly, SD) in the chain: Step 1: express the formula_22's in a local basis for qubit 2: The vectors formula_24 are not necessarily normalized. Step 2: write each vector formula_24 in terms of the at most (Vidal's emphasis) formula_26 Schmidt vectors formula_27 and, correspondingly, coefficients formula_28 The Schmidt eigenvectors are simply: and Rationale. Now, looking at "D", instead  of formula_31 initial terms, there are formula_32. Apparently this is just a fancy way of rewriting the coefficients formula_33, but in fact there is more to it than that. Assuming that "N" is even, the Schmidt rank formula_26 for a bipartite cut in the middle of the chain can have a maximal value of formula_35; in this case we end up with at least formula_36 coefficients, considering only the formula_37 ones, slightly more than the initial formula_31! The truth is that the decomposition "D" is useful when dealing with systems that exhibit a low degree of entanglement, which fortunately is  the case with many 1D systems, where the Schmidt coefficients of the ground state decay in an exponential manner with formula_39: The numerical simulation. The numerical simulation is targeting (possibly time-dependent) Hamiltonians of a system of formula_47 particles arranged in a line, which are composed of arbitrary OQGs and TQGs: It is useful to decompose formula_49 as a sum of two possibly non-commuting terms, formula_50, where Any two-body terms commute: formula_53, formula_54 This is done to make the Suzuki–Trotter expansion (ST) of the exponential operator, named after Masuo Suzuki and Hale Trotter. The Suzuki–Trotter expansion . The Suzuki–Trotter expansion of the first order (ST1) represents a general way of writing exponential operators: or, equivalently Our goal is to make the time evolution of a state formula_58 for a time T, towards the state formula_59 using the n-particle Hamiltonian formula_60. It is rather troublesome, if at all possible, to construct the decomposition formula_57 for an arbitrary n-particle state, since this would mean one has to compute the Schmidt decomposition at each bond, to arrange the Schmidt eigenvalues in decreasing order and to choose the first formula_62 and the appropriate Schmidt eigenvectors. Mind this would  imply diagonalizing somewhat generous reduced density matrices, which, depending on the system one has to simulate, might be a task beyond our reach and patience. Instead, one can try to do the following: i) construct the decomposition formula_57 for a simple initial state, let us say, some product state formula_64, for which the decomposition is straightforward. ii) relate formula_58 to the ground state formula_66 of a Hamiltonian formula_67 by a sufficiently local transformation Q (one that can be expressed as a product of TQGs, for example) formula_68 iii) make an imaginary-time evolution towards the ground state of the Hamiltonian  formula_67, formula_66 according to: or, alternatively, simulate an isentropic evolution using a time-dependent Hamiltonian, which interpolates between the Hamiltonian formula_72, which has the product state formula_64 as its ground state, and the Hamiltonian formula_67; the evolution must be done slow enough, such that the system is always in the ground state or, at least, very close to it. iv)finally, make the time-evolution of the state formula_58 towards formula_59 using the Hamiltonian formula_60: Error sources. The errors in the simulation are resulting from the Suzuki–Trotter approximation and the involved truncation of the Hilbert space. Errors coming from the  Suzuki–Trotter expansion. In the case of a Trotter approximation of formula_79 order, the error is of order formula_80. Taking into account formula_81 steps, the error after the time T is: The unapproximated state formula_83 is: where formula_85 is the state kept after the Trotter expansion and formula_86 accounts for the part that is neglected when doing the expansion. The total error scales with time formula_87 as: The Trotter error is independent of the dimension of the chain. Errors coming from the truncation of the Hilbert space. Considering the errors arising from the truncation of the Hilbert space comprised in  the decomposition "D", they are twofold. First, as we have seen above, the smallest contributions to the Schmidt spectrum are left away, the state being faithfully represented up to: where formula_90 is the sum of all the discarded eigenvalues of the reduced density matrix, at the bond formula_91. The state formula_92 is, at a given bond formula_91, described by the Schmidt decomposition: where is the state formed by the eigenfunctions corresponding to the smallest, irrelevant Schmidt coefficients, which are neglected. Now, formula_96 because they are spanned by vectors corresponding to orthogonal spaces. Using the same argument as for the Trotter  expansion, the error after the truncation is: After moving to the next bond, the state is, similarly: The error, after the second truncation, is: and so on, as we move from bond to bond. The second error source enfolded in the decomposition formula_57 is more subtle and requires a little bit of calculation. As we calculated before, the normalization constant after making the truncation at bond formula_101 formula_102 is: where formula_104. Taking into account the truncated space, the norm is: Taking the difference, formula_106, we get: Hence, when constructing the reduced density matrix, the trace of the matrix is multiplied  by the factor: The total truncation error. The total truncation error, considering both sources, is upper bounded by: When using the Trotter expansion, we do not move from bond to bond, but between bonds of same parity; moreover, for the ST2, we make a sweep of the even ones and two for the odd. But nevertheless, the calculation presented above still holds. The error is evaluated by successively multiplying with the normalization constant, each time we build the reduced density matrix and select its relevant eigenvalues. Section::::"Adaptive" Schmidt dimension. One thing that can save a lot of computational time without  loss of accuracy is to use a different Schmidt dimension for each bond instead of a fixed one for all bonds, keeping only the necessary amount of relevant coefficients, as usual. For example, taking the first bond, in the case of qubits, the Schmidt dimension is just two. Hence, at the first bond, instead of futilely diagonalizing, let us say, 10 by 10 or 20 by 20 matrices, we can just restrict ourselves to ordinary 2 by 2 ones, thus making the algorithm generally faster. What we can do instead is set a threshold for the eigenvalues of the SD , keeping only those that are above the threshold. TEBD also offers the possibility of straightforward parallelization due to the factorization of the exponential time-evolution operator using the Suzuki–Trotter expansion. A parallel-TEBD has the same mathematics as its non-parallelized counterpart, the only difference is in the numerical implementation.