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 state 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 Time-evolving block decimation is a numerical algorithm that can efficiently simulate the time evolution of one dimensional quantum systems having limited entanglement entropy. Naively, to time evolve a system characterized by a Hamiltonian H ^ {\displaystyle {\hat {H}}} one would directly exponentiate the Hamiltonian to get the time evolution operator U ^ ( t ) = e − i H ^ t / ℏ {\displaystyle {\hat {U}}(t)=e^{-i{\hat {H}}t/\hbar }} and apply this to an initial state: ψ ( t ) = U ^ ( t ) ψ ( t = 0 ) {\displaystyle \psi (t)={\hat {U}}(t)\psi (t=0)} However, as the number of degrees of freedom of the system grows it quickly becomes computationally infeasible to perform the associated matrix exponentiation and matrix-vector multiplication. For example, if ψ {\displaystyle \psi } represents a system of n {\displaystyle n} qubits then the Hilbert space in which ψ {\displaystyle \psi } resides has dimension 2 n {\displaystyle 2^{n}} , meaning matrix operations are effectively intractable for all but the smallest values of n {\displaystyle n} . TEBD presents an efficient scheme for performing time evolution by limiting itself to a much smaller subspace of the configuration space. There are several other noteworthy examples of ways to get around this exponential scaling, including quantum Monte Carlo and the density matrix renormalization group. Guifré Vidal proposed the scheme while at the Institute for Quantum Information, Caltech. 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 for a wide suite of Hamiltonians characterized by local interactions, 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 one to work in a re-sized space spanned by the eigenvectors corresponding to the selected eigenvalues. 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. Other groups have developed similar approaches in which quantum information plays a predominant role: for example, in DMRG implementations for periodic boundary conditions, and for studying mixed-state dynamics in one-dimensional quantum lattice systems. 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 | Ψ ⟩ ∈ H ⊗ N {\displaystyle |\Psi \rangle \in H^{{\otimes }N}} . The most natural way of describing | Ψ ⟩ {\displaystyle |\Psi \rangle } would be using the local M N {\displaystyle M^{N}} -dimensional basis | i 1 , i 2 , . . , i N − 1 , i N ⟩ {\displaystyle |i_{1},i_{2},..,i_{N-1},i_{N}\rangle } :
| Ψ ⟩ = ∑ i = 1 M c i 1 i 2 . . i N | i 1 , i 2 , . . , i N − 1 , i N ⟩ {\displaystyle |\Psi \rangle =\sum \limits _{i=1}^{M}c_{i_{1}i_{2}..i_{N}}|{i_{1},i_{2},..,i_{N-1},i_{N}}\rangle }
where M is the on-site dimension. The trick of TEBD is to re-write the coefficients c i 1 i 2 . . i N {\displaystyle c_{i_{1}i_{2}..i_{N}}} :
c i 1 i 2 . . i N = ∑ α 1 , . . , α N − 1 = 0 χ Γ α 1 [ 1 ] i 1 λ α 1 [ 1 ] Γ α 1 α 2 [ 2 ] i 2 λ α 2 [ 2 ] Γ α 2 α 3 [ 3 ] i 3 λ α 3 [ 3 ] ⋅ . . ⋅ Γ α N − 2 α N − 1 [ N − 1 ] i N − 1 λ α N − 1 [ N − 1 ] Γ α N − 1 [ N ] i N {\displaystyle c_{i_{1}i_{2}..i_{N}}=\sum \limits _{\alpha _{1},..,\alpha _{N-1}=0}^{\chi }\Gamma _{\alpha _{1}}^{[1]i_{1}}\lambda _{\alpha _{1}}^{[1]}\Gamma _{\alpha _{1}\alpha _{2}}^{[2]i_{2}}\lambda _{\alpha _{2}}^{[2]}\Gamma _{\alpha _{2}\alpha _{3}}^{[3]i_{3}}\lambda _{\alpha _{3}}^{[3]}\cdot ..\cdot \Gamma _{\alpha _{N-2}\alpha _{N-1}}^{[{N-1}]i_{N-1}}\lambda _{\alpha _{N-1}}^{[N-1]}\Gamma _{\alpha _{N-1}}^{[N]i_{N}}}
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 | Ψ ⟩ ∈ H A ⊗ H B {\displaystyle \vert \Psi \rangle \in {H_{A}\otimes H_{B}}} . Every such state | Ψ ⟩ {\displaystyle |{\Psi }\rangle } can be represented in an appropriately chosen basis as:
| Ψ ⟩ = ∑ i = 1 M A | B a i | Φ i A Φ i B ⟩ {\displaystyle \left\vert \Psi \right\rangle =\sum \limits _{i=1}^{M_{A|B}}a_{i}\left\vert {\Phi _{i}^{A}\Phi _{i}^{B}}\right\rangle }
where | Φ i A Φ i B ⟩ = | Φ i A ⟩ ⊗ | Φ i B ⟩ {\displaystyle |{\Phi _{i}^{A}\Phi _{i}^{B}}\rangle =|{\Phi _{i}^{A}}\rangle \otimes |{\Phi _{i}^{B}}\rangle } are formed with vectors | Φ i A ⟩ {\displaystyle |{\Phi _{i}^{A}}\rangle } that make an orthonormal basis in H A {\displaystyle H_{A}} and, correspondingly, vectors | Φ i B ⟩ {\displaystyle |{\Phi _{i}^{B}}\rangle } , which form an orthonormal basis in H B {\displaystyle {H_{B}}} , with the coefficients a i {\displaystyle a_{i}} being real and positive, ∑ i = 1 M A | B a i 2 = 1 {\textstyle \sum \limits _{i=1}^{M_{A|B}}a_{i}^{2}=1} . This is called the Schmidt decomposition (SD) of a state. In general the summation goes up to M A | B = min ( dim ( H A ) , dim ( H B ) ) {\displaystyle M_{A|B}=\min(\dim({H_{A}}),\dim({H_{B}}))} . 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:
| Ψ ⟩ = 1 2 2 ( | 00 ⟩ + 3 | 01 ⟩ + 3 | 10 ⟩ + | 11 ⟩ ) {\displaystyle |{\Psi }\rangle ={\frac {1}{2{\sqrt {2}}}}\left(|{00}\rangle +{\sqrt {3}}|{01}\rangle +{\sqrt {3}}|{10}\rangle +|{11}\rangle \right)}
has the following SD:
| Ψ ⟩ = 3 + 1 2 2 | ϕ 1 A ϕ 1 B ⟩ + 3 − 1 2 2 | ϕ 2 A ϕ 2 B ⟩ {\displaystyle \left|{\Psi }\right\rangle ={\frac {{\sqrt {3}}+1}{2{\sqrt {2}}}}\left|{\phi _{1}^{A}\phi _{1}^{B}}\right\rangle +{\frac {{\sqrt {3}}-1}{2{\sqrt {2}}}}\left|{\phi _{2}^{A}\phi _{2}^{B}}\right\rangle }
with
| ϕ 1 A ⟩ = 1 2 ( | 0 A ⟩ + | 1 A ⟩ ) , | ϕ 1 B ⟩ = 1 2 ( | 0 B ⟩ + | 1 B ⟩ ) , | ϕ 2 A ⟩ = 1 2 ( | 0 A ⟩ − | 1 A ⟩ ) , | ϕ 2 B ⟩ = 1 2 ( | 1 B ⟩ − | 0 B ⟩ ) {\displaystyle |{\phi _{1}^{A}}\rangle ={\frac {1}{\sqrt {2}}}(|{0_{A}}\rangle +|{1_{A}}\rangle ),\ \ |{\phi _{1}^{B}}\rangle ={\frac {1}{\sqrt {2}}}(|{0_{B}}\rangle +|{1_{B}}\rangle ),\ \ |{\phi _{2}^{A}}\rangle ={\frac {1}{\sqrt {2}}}(|{0_{A}}\rangle -|{1_{A}}\rangle ),\ \ |{\phi _{2}^{B}}\rangle ={\frac {1}{\sqrt {2}}}(|{1_{B}}\rangle -|{0_{B}}\rangle )}
On the other hand, the state:
| Φ ⟩ = 1 3 | 00 ⟩ + 1 6 | 01 ⟩ − i 3 | 10 ⟩ − i 6 | 11 ⟩ {\displaystyle |{\Phi }\rangle ={\frac {1}{\sqrt {3}}}|{00}\rangle +{\frac {1}{\sqrt {6}}}|{01}\rangle -{\frac {i}{\sqrt {3}}}|{10}\rangle -{\frac {i}{\sqrt {6}}}|{11}\rangle }
is a product state:
| Φ ⟩ = ( 1 3 | 0 A ⟩ − i 3 | 1 A ⟩ ) ⊗ ( | 0 B ⟩ + 1 2 | 1 B ⟩ ) {\displaystyle \left|\Phi \right\rangle =\left({\frac {1}{\sqrt {3}}}\left|0_{A}\right\rangle -{\frac {i}{\sqrt {3}}}\left|1_{A}\right\rangle \right)\otimes \left(\left|0_{B}\right\rangle +{\frac {1}{\sqrt {2}}}\left|1_{B}\right\rangle \right)}
Building the decomposition of state At this point we know enough to try to see how we explicitly build the decomposition (let's call it D). Consider the bipartite splitting [ 1 ] : [ 2.. N ] {\displaystyle [1]:[2..N]} . The SD has the coefficients λ α 1 [ 1 ] {\displaystyle \lambda _{{\alpha }_{1}}^{[1]}} and eigenvectors | Φ α 1 [ 1 ] ⟩ | Φ α 1 [ 2.. N ] ⟩ {\displaystyle \left|{\Phi _{\alpha _{1}}^{[1]}}\right\rangle \left|{\Phi _{\alpha _{1}}^{[2..N]}}\right\rangle } . By expanding the | Φ α 1 [ 1 ] ⟩ {\displaystyle \left|{\Phi _{\alpha _{1}}^{[1]}}\right\rangle } 's in the local basis, one can write:
| Ψ ⟩ = ∑ i 1 , α 1 = 1 M , χ Γ α 1 [ 1 ] i 1 λ α 1 [ 1 ] | i 1 ⟩ | Φ α 1 [ 2.. N ] ⟩ {\displaystyle |{\Psi }\rangle =\sum \limits _{i_{1},{\alpha _{1}=1}}^{M,\chi }\Gamma _{\alpha _{1}}^{[1]i_{1}}\lambda _{\alpha _{1}}^{[1]}|{i_{1}}\rangle |{\Phi _{\alpha _{1}}^{[2..N]}}\rangle }
The process can be decomposed in three steps, iterated for each bond (and, correspondingly, SD) in the chain: Step 1: express the | Φ α 1 [ 2.. N ] ⟩ {\displaystyle |{\Phi _{\alpha _{1}}^{[2..N]}}\rangle } 's in a local basis for qubit 2:
| Φ α 1 [ 2.. N ] ⟩ = ∑ i 2 | i 2 ⟩ | τ α 1 i 2 [ 3.. N ] ⟩ {\displaystyle |{\Phi _{\alpha _{1}}^{[2..N]}}\rangle =\sum _{i_{2}}|{i_{2}}\rangle |{\tau _{\alpha _{1}i_{2}}^{[3..N]}}\rangle }
The vectors | τ α 1 i 2 [ 3.. N ] ⟩ {\displaystyle |{\tau _{\alpha _{1}i_{2}}^{[3..N]}}\rangle } are not necessarily normalized. Step 2: write each vector | τ α 1 i 2 [ 3.. N ] ⟩ {\displaystyle |{\tau _{\alpha _{1}i_{2}}^{[3..N]}}\rangle } in terms of the at most (Vidal's emphasis) χ {\displaystyle \chi } Schmidt vectors | Φ α 2 [ 3.. N ] ⟩ {\displaystyle |{\Phi _{\alpha _{2}}^{[3..N]}}\rangle } and, correspondingly, coefficients λ α 2 [ 2 ] {\displaystyle \lambda _{{\alpha }_{2}}^{[2]}} :
| τ α 1 i 2 [ 3.. N ] ⟩ = ∑ α 2 Γ α 1 α 2 [ 2 ] i 2 λ α 2 [ 2 ] | Φ α 2 [ 3.. N ] ⟩ {\displaystyle |\tau _{\alpha _{1}i_{2}}^{[3..N]}\rangle =\sum _{\alpha _{2}}\Gamma _{\alpha _{1}\alpha _{2}}^{[2]i_{2}}\lambda _{{\alpha }_{2}}^{[2]}|{\Phi _{\alpha _{2}}^{[3..N]}}\rangle }
Step 3: make the substitutions and obtain:
| Ψ ⟩ = ∑ i 1 , i 2 , α 1 , α 2 Γ α 1 [ 1 ] i 1 λ α 1 [ 1 ] Γ α 1 α 2 [ 2 ] i 2 λ α 2 [ 2 ] | i 1 i 2 ⟩ | Φ α 2 [ 3.. N ] ⟩ {\displaystyle |{\Psi }\rangle =\sum _{i_{1},i_{2},\alpha _{1},\alpha _{2}}\Gamma _{\alpha _{1}}^{[1]i_{1}}\lambda _{\alpha _{1}}^{[1]}\Gamma _{\alpha _{1}\alpha _{2}}^{[2]i_{2}}\lambda _{{\alpha }_{2}}^{[2]}|{i_{1}i_{2}}\rangle |{\Phi _{\alpha _{2}}^{[3..N]}}\rangle }
Repeating the steps 1 to 3, one can construct the whole decomposition of state D. The last Γ {\displaystyle \Gamma } 's are a special case, like the first ones, expressing the right-hand Schmidt vectors at the ( N − 1 ) t h {\displaystyle (N-1)^{th}} bond in terms of the local basis at the N t h {\displaystyle N^{th}} lattice place. As shown in, it is straightforward to obtain the Schmidt decomposition at k t h {\displaystyle k^{th}} bond, i.e. [ 1.. k ] : [ k + 1.. N ] {\displaystyle [1..k]:[k+1..N]} , from D. The Schmidt eigenvalues, are given explicitly in D:
| Ψ ⟩ = ∑ α k λ α k [ k ] | Φ α k [ 1.. k ] ⟩ | Φ α k [ k + 1.. N ] ⟩ {\displaystyle |{\Psi }\rangle =\sum _{\alpha _{k}}\lambda _{{\alpha }_{k}}^{[k]}|{\Phi _{\alpha _{k}}^{[1..k]}}\rangle |{\Phi _{\alpha _{k}}^{[k+1..N]}}\rangle }
The Schmidt eigenvectors are simply:
| Φ α