Introducing the decomposition of State

Consider a chain of N qubits, described by the function ${\displaystyle |\Psi \rangle \in H^{{\otimes }N}}$. The most natural way of describing ${\displaystyle |\Psi \rangle }$ would be using the local ${\displaystyle M^{N}}$-dimensional basis ${\displaystyle |i_{1},i_{2},..,i_{N-1},i_{N}\rangle }$:

${\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 ${\displaystyle c_{i_{1}i_{2}..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 ${\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:

${\displaystyle |{\Psi }\rangle =\sum \limits _{i=1}^{M_{A|B}}a_{i}|{\Phi _{i}^{A}\Phi _{i}^{B}}\rangle }$

where ${\displaystyle |{\Phi _{i}^{A}\Phi _{i}^{B}}\rangle =|{\Phi _{i}^{A}}\rangle \otimes |{\Phi _{i}^{B}}\rangle }$ are formed with vectors ${\displaystyle |{\Phi _{i}^{A}}\rangle }$ that make an orthonormal basis in ${\displaystyle H_{A}}$ and, correspondingly, vectors ${\displaystyle |{\Phi _{i}^{B}}\rangle }$, which form an orthonormal basis in ${\displaystyle {H_{B}}}$, with the coefficients ${\displaystyle a_{i}}$ being real and positive, ${\displaystyle \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 ${\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:

${\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:

${\displaystyle |{\Psi }\rangle ={\frac {{\sqrt {3}}+1}{2{\sqrt {2}}}}|{\phi _{1}^{A}\phi _{1}^{B}}\rangle +{\frac {{\sqrt {3}}-1}{2{\sqrt {2}}}}|{\phi _{2}^{A}\phi _{2}^{B}}\rangle }$

with

${\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:

${\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:

${\displaystyle |{\Phi }\rangle ={(}{\frac {1}{\sqrt {3}}}|{0_{A}}\rangle -{\frac {i}{\sqrt {3}}}|{1_{A}}\rangle {)}\otimes {(}|{0_{B}}\rangle +{\frac {1}{\sqrt {2}}}|{1_{B}}\rangle {)}}$

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 ${\displaystyle [1]:[2..N]}$. The SD has the coefficients ${\displaystyle \lambda _{{\alpha }_{1}}^{[1]}}$ and eigenvectors ${\displaystyle |{\Phi _{\alpha _{1}}^{[1]}}\rangle |{\Phi _{\alpha _{1}}^{[2..N]}}\rangle }$. By expanding the ${\displaystyle |{\Phi _{\alpha _{1}}^{[1]}}\rangle }$'s in the local basis, one can write:

${\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 ${\displaystyle |{\Phi _{\alpha _{1}}^{[2..N]}}\rangle }$'s in a local basis for qubit 2:

${\displaystyle |{\Phi _{\alpha _{1}}^{[2..N]}}\rangle =\sum _{i_{2}}|{i_{2}}\rangle |{\tau _{\alpha _{1}i_{2}}^{[3..N]}}\rangle }$

The vectors ${\displaystyle |{\tau _{\alpha _{1}i_{2}}^{[3..N]}}\rangle }$ are not necessarily normalized.

Step 2: write each vector ${\displaystyle |{\tau _{\alpha _{1}i_{2}}^{[3..N]}}\rangle }$ in terms of the at most (Vidal's emphasis) ${\displaystyle \chi }$ Schmidt vectors ${\displaystyle |{\Phi _{\alpha _{2}}^{[3..N]}}\rangle }$ and, correspondingly, coefficients ${\displaystyle \lambda _{{\alpha }_{2}}^{[2]}}$:

${\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:

${\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 ${\displaystyle (N-1)^{th}}$ bond in terms of the local basis at the ${\displaystyle N^{th}}$ lattice place. As shown in,^[1] it is straightforward to obtain the Schmidt decomposition at ${\displaystyle k^{th}}$ bond, i.e. ${\displaystyle [1..k]:[k+1..N]}$, from D.

The Schmidt eigenvalues, are given explicitly in D:

${\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:

${\displaystyle |{\Phi _{\alpha _{k}}^{[1..k]}}\rangle =\sum _{\alpha _{1},\alpha _{2}..\alpha _{k-1}}\Gamma _{\alpha _{1}}^{[1]i_{1}}\lambda _{\alpha _{1}}^{[1]}\cdot \cdot \Gamma _{\alpha _{k-1}\alpha _{k}}^{[k]i_{k}}|{i_{1}i_{2}..i_{k}}\rangle }$

and

${\displaystyle |{\Phi _{\alpha _{k}}^{[k+1..N]}}\rangle =\sum _{\alpha _{k+1},\alpha _{k+2}..\alpha _{N}}\Gamma _{\alpha _{k}\alpha _{k+1}}^{[k+1]i_{k+1}}\lambda _{\alpha _{k+1}}^{[k+1]}\cdot \cdot \lambda _{\alpha _{N-1}}^{N-1}\Gamma _{\alpha _{N-1}}^{[N]i_{N}}|{i_{k+1}i_{k+2}..i_{N}}\rangle }$

Rationale

Now, looking at D, instead of ${\displaystyle M^{N}}$ initial terms, there are ${\displaystyle {\chi }^{2}{\cdot }M(N-2)+2{\chi }M+(N-1)\chi }$. Apparently this is just a fancy way of rewriting the coefficients ${\displaystyle c_{i_{1}i_{2}..i_{N}}}$, but in fact there is more to it than that. Assuming that N is even, the Schmidt rank ${\displaystyle \chi }$ for a bipartite cut in the middle of the chain can have a maximal value of ${\displaystyle M^{N/2}}$; in this case we end up with at least ${\displaystyle M^{N+1}{\cdot }(N-2)}$ coefficients, considering only the ${\displaystyle {\chi }^{2}}$ ones, slightly more than the initial ${\displaystyle M^{N}}$! 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 ${\displaystyle \alpha }$:

${\displaystyle \lambda _{{\alpha }_{l}}^{[l]}{\sim }e^{-K\alpha _{l}},\ K>0.}$

Therefore, it is possible to take into account only some of the Schmidt coefficients (namely the largest ones), dropping the others and consequently normalizing again the state:

${\displaystyle |{\Psi }\rangle ={\frac {1}{\sqrt {\sum \limits _{{\alpha _{l}}=1}^{{\chi }_{c}}{|\lambda _{{\alpha }_{l}}^{[l]}|}^{2}}}}\cdot \sum \limits _{{{\alpha }_{l}}=1}^{{\chi }_{c}}\lambda _{{\alpha }_{l}}^{[l]}|{\Phi _{\alpha _{l}}^{[1..l]}}\rangle |{\Phi _{\alpha _{l}}^{[l+1..N]}}\rangle ,}$

where ${\displaystyle \chi _{c}}$ is the number of kept Schmidt coefficients.

Let's get away from this abstract picture and refresh ourselves with a concrete example, to emphasize the advantage of making this decomposition. Consider for instance the case of 50 fermions in a ferromagnetic chain, for the sake of simplicity. A dimension of 12, let's say, for the ${\displaystyle \chi _{c}}$ would be a reasonable choice, keeping the discarded eigenvalues at ${\displaystyle 0.0001}$% of the total, as shown by numerical studies,^[4] meaning roughly ${\displaystyle 2^{14}}$ coefficients, as compared to the originally ${\displaystyle 2^{50}}$ ones.

Even if the Schmidt eigenvalues don't have this exponential decay, but they show an algebraic decrease, we can still use D to describe our state ${\displaystyle \psi }$. The number of coefficients to account for a faithful description of ${\displaystyle \psi }$ may be sensibly larger, but still within reach of eventual numerical simulations.

One-qubit gates acting on qubit k

The OQGs are affecting only the qubit they are acting upon, the update of the state ${\displaystyle |{\psi }\rangle }$ after a unitary operator at qubit k does not modify the Schmidt eigenvalues or vectors on the left, consequently the ${\displaystyle \Gamma ^{[k-1]}}$'s, or on the right, hence the ${\displaystyle \Gamma ^{[k+1]}}$'s. The only ${\displaystyle \Gamma }$'s that will be updated are the ${\displaystyle \Gamma ^{[k]}}$'s (requiring only at most ${\displaystyle {O}(M^{2}\cdot \chi ^{2})}$ operations), as

${\displaystyle \Gamma _{\alpha _{k-1}\alpha _{k}}^{'[k]i_{k}}=\sum _{j}U_{j_{k}}^{i_{k}}\Gamma _{\alpha _{k-1}\alpha _{k}}^{[k]j_{k}}.}$

Two-qubit gates acting on qubits k, k+1

The changes required to update the ${\displaystyle \Gamma }$'s and the ${\displaystyle \lambda }$'s, following a unitary operation V on qubits k, k+1, concern only ${\displaystyle \Gamma ^{[k]}}$, and ${\displaystyle \Gamma ^{[k+1]}}$. They consist of a number of ${\displaystyle {O}({M\cdot \chi }^{3})}$ basic operations.

Following Vidal's original approach, ${\displaystyle |{\psi }\rangle }$ can be regarded as belonging to only four subsystems:

${\displaystyle {{H}=J{\otimes }H_{C}{\otimes }H_{D}{\otimes }K}.\,}$

The subspace J is spanned by the eigenvectors of the reduced density matrix ${\displaystyle \rho ^{J}=Tr_{CDK}|\psi \rangle \langle \psi |}$:

${\displaystyle \rho ^{[1..{k-1}]}=\sum _{\alpha }{(\lambda _{\alpha }^{[k-1]})}^{2}|{\Phi _{\alpha }^{[1..{k-1}]}}\rangle \langle {\Phi _{\alpha }^{[1..{k-1}]}}|=\sum _{\alpha }{(\lambda _{\alpha }^{[k-1]})^{2}}|{\alpha }\rangle \langle {\alpha }|.}$

In a similar way, the subspace K is spanned by the eigenvectors of the reduced density matrix:

${\displaystyle \rho ^{[{k+2}..{N}]}=\sum _{\gamma }{(\lambda _{\gamma }^{[k+1]})^{2}}|{\Phi _{\gamma }^{[{k+2}..N]}}\rangle \langle {\Phi _{\gamma }^{[{k+2}..N]}}|=\sum _{\gamma }{(\lambda _{\gamma }^{[k+1]})^{2}}|{\gamma }\rangle \langle {\gamma }|.}$

The subspaces ${\displaystyle H_{C}}$ and ${\displaystyle H_{D}}$ belong to the qubits k and k + 1. Using this basis and the decomposition D, ${\displaystyle |{\psi }\rangle }$ can be written as:

${\displaystyle |{\psi }\rangle =\sum \limits _{\alpha ,\beta ,\gamma =1}^{\chi }\sum \limits _{i,j=1}^{M}\lambda _{\alpha }^{[C-1]}\Gamma _{\alpha \beta }^{[C]i}\lambda _{\beta }^{[C]}\Gamma _{\beta \gamma }^{[D]j}\lambda _{\gamma }^{[D]}|{{\alpha }ij{\gamma }}\rangle }$

Using the same reasoning as for the OQG, the applying the TQG V to qubits k, k + 1 one needs only to update

${\displaystyle \Gamma ^{[C]}}$, ${\displaystyle \lambda }$ and ${\displaystyle \Gamma ^{[D]}.}$

We can write ${\displaystyle |{\psi '}\rangle =V|{\psi }\rangle }$ as:

${\displaystyle |{\psi '}\rangle =\sum \limits _{\alpha ,\gamma =1}^{\chi }\sum \limits _{i,j=1}^{M}\lambda _{\alpha }\Theta _{\alpha \gamma }^{ij}\lambda _{\gamma }|{{\alpha }ij\gamma }\rangle }$

where

${\displaystyle \Theta _{\alpha \gamma }^{ij}=\sum \limits _{\beta =1}^{\chi }\sum \limits _{m,n=1}^{M}V_{mn}^{ij}\Gamma _{\alpha \beta }^{[C]m}\lambda _{\beta }\Gamma _{\beta \gamma }^{[D]n}.}$

To find out the new decomposition, the new ${\displaystyle \lambda }$'s at the bond k and their corresponding Schmidt eigenvectors must be computed and expressed in terms of the ${\displaystyle {\Gamma }}$'s of the decomposition D. The reduced density matrix ${\displaystyle \rho ^{'[DK]}}$ is therefore diagonalized:

${\displaystyle \rho ^{'[DK]}=Tr_{JC}|{\psi '}\rangle \langle {\psi '}|=\sum _{j,j',\gamma ,\gamma '}\rho _{\gamma \gamma '}^{jj'}|{j\gamma }\rangle \langle {j'\gamma '}|.}$

The square roots of its eigenvalues are the new ${\displaystyle \lambda }$'s. Expressing the eigenvectors of the diagonalized matrix in the basis :${\displaystyle \{|{j\gamma }\rangle \}}$ the ${\displaystyle \Gamma ^{[{D]}}}$'s are obtained as well:

${\displaystyle |{\Phi ^{'[{DK}]}}\rangle =\sum _{j,\gamma }\Gamma _{\beta \gamma }^{'[{D}]j}\lambda _{\gamma }|{j\gamma }\rangle .}$

From the left-hand eigenvectors,

${\displaystyle \lambda _{\beta }^{'}|{\Phi _{\beta }^{'[{JC}]}}\rangle =\langle {\Phi _{\beta }^{'[{DK}]}}|{\psi '}\rangle =\sum _{i,j,\alpha ,\gamma }(\Gamma _{\beta \gamma }^{'[{D}]j})^{*}\Theta _{\alpha \gamma }^{ij}(\lambda _{\gamma })^{2}\lambda _{\alpha }|{{\alpha }i}\rangle }$

after expressing them in the basis ${\displaystyle \{|{i\alpha }\rangle \}}$, the ${\displaystyle \Gamma ^{[{C}]}}$'s are:

${\displaystyle |{\Phi ^{'[{JC}]}}\rangle =\sum _{i,\alpha }\Gamma _{\alpha \beta }^{'[{C}]i}\lambda _{\alpha }|{{\alpha }i}\rangle .}$

The computational cost

The dimension of the largest tensors in D is of the order ${\displaystyle {O}(M{\cdot }{\chi }^{2})}$; when constructing the ${\displaystyle \Theta _{\alpha \gamma }^{ij}}$ one makes the summation over ${\displaystyle \beta }$, ${\displaystyle {\it {m}}}$ and ${\displaystyle {\it {n}}}$ for each ${\displaystyle \gamma ,\alpha ,{\it {i,j}}}$, adding up to a total of ${\displaystyle {O}(M^{4}{\cdot }{\chi }^{3})}$ operations. The same holds for the formation of the elements ${\displaystyle \rho _{\gamma \gamma '}^{jj'}}$, or for computing the left-hand eigenvectors ${\displaystyle \lambda _{\beta }^{'}|{\Phi _{\beta }^{'[{\it {JC}}]}}\rangle }$, a maximum of ${\displaystyle {\it {O}}(M^{3}{\cdot }{\chi }^{3})}$, respectively ${\displaystyle {\it {O}}(M^{2}{\cdot }{\chi }^{3})}$ basic operations. In the case of qubits, ${\displaystyle {\it {M}}=2}$, hence its role is not very relevant for the order of magnitude of the number of basic operations, but in the case when the on-site dimension is higher than two it has a rather decisive contribution.

The Suzuki–Trotter expansion

The Suzuki–Trotter expansion of the first order (ST1) represents a general way of writing exponential operators:

${\displaystyle e^{(A+B)}=\lim _{n\rightarrow \infty }{\Bigl (}e^{\frac {A}{n}}e^{\frac {B}{n}}{\Bigr )}^{n}}$

or, equivalently

${\displaystyle e^{{\delta }(A+B)}=e^{{\delta }A}e^{{\delta }B}+{\it {O}}(\delta ^{2}).}$

The correction term vanishes in the limit ${\displaystyle \delta \rightarrow 0}$

For simulations of quantum dynamics it is useful to use operators that are unitary, conserving the norm (unlike power series expansions), and there's where the Trotter-Suzuki expansion comes in. In problems of quantum dynamics the unitarity of the operators in the ST expansion proves quite practical, since the error tends to concentrate in the overall phase, thus allowing us to faithfully compute expectation values and conserved quantities. Because the ST conserves the phase-space volume, it is also called a symplectic integrator.

The trick of the ST2 is to write the unitary operators ${\displaystyle e^{-iHt}}$ as:

${\displaystyle e^{-iH_{N}T}=[e^{-iH_{N}\delta }]^{T/{\delta }}=[e^{{\frac {\delta }{2}}F}e^{{\delta }G}e^{{\frac {\delta }{2}}F}]^{n}}$

where ${\displaystyle n={\frac {T}{\delta }}}$. The number ${\displaystyle {\it {n}}}$ is called the Trotter number.

Simulation of the time-evolution

The operators ${\displaystyle e^{{\frac {\delta }{2}}F}}$, ${\displaystyle e^{{\delta }G}}$ are easy to express, as:

${\displaystyle e^{{\frac {\delta }{2}}F}=\prod _{even\ \ l}e^{{\frac {\delta }{2}}F^{[l]}}}$

${\displaystyle e^{{\delta }G}=\prod _{odd\ \ l}e^{{\delta }G^{[l]}}}$

since any two operators ${\displaystyle F^{[l]}}$,${\displaystyle F^{[l']}}$ (respectively, ${\displaystyle G^{[l]}}$,${\displaystyle G^{[l']}}$) commute for ${\displaystyle l{\neq }l'}$ and an ST expansion of the first order keeps only the product of the exponentials, the approximation becoming, in this case, exact.

The time-evolution can be made according to

${\displaystyle |{{\tilde {\psi }}_{t+\delta }}\rangle =e^{-i{\frac {\delta }{2}}F}e^{{-i\delta }G}e^{{\frac {-i\delta }{2}}F}|{{\tilde {\psi }}_{t}}\rangle .}$

For each "time-step" ${\displaystyle \delta }$, ${\displaystyle e^{-i{\frac {\delta }{2}}F^{[l]}}}$ are applied successively to all odd sites, then ${\displaystyle e^{{-i\delta }G^{[l]}}}$ to the even ones, and ${\displaystyle e^{-i{\frac {\delta }{2}}F^{[l]}}}$ again to the odd ones; this is basically a sequence of TQG's, and it has been explained above how to update the decomposition ${\displaystyle {\it {D}}}$ when applying them.

Our goal is to make the time evolution of a state ${\displaystyle |{\psi _{0}}\rangle }$ for a time T, towards the state ${\displaystyle |{\psi _{T}}\rangle }$ using the n-particle Hamiltonian ${\displaystyle H_{n}}$.

It is rather troublesome, if at all possible, to construct the decomposition ${\displaystyle {\it {D}}}$ 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 ${\displaystyle \chi _{c}}$ 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 ${\displaystyle {\it {D}}}$ for a simple initial state, let us say, some product state ${\displaystyle |{\psi _{P}}\rangle }$, for which the decomposition is straightforward.

ii) relate ${\displaystyle |{\psi _{0}}\rangle }$ to the ground state ${\displaystyle |{\psi _{gr}}\rangle }$ of a Hamiltonian ${\displaystyle {\tilde {H}}}$ by a sufficiently local transformation Q (one that can be expressed as a product of TQGs, for example) ${\displaystyle |{\psi _{0}}\rangle =Q|{\psi _{gr}}\rangle }$

iii) make an imaginary-time evolution towards the ground state of the Hamiltonian ${\displaystyle {\tilde {H}}}$, ${\displaystyle |{\psi _{gr}}\rangle }$ according to:

${\displaystyle |{\psi _{gr}}\rangle =\lim _{\tau \rightarrow \infty }{\frac {e^{-{\tilde {H}}\tau }|{\psi _{P}}\rangle }{||e^{-{\tilde {H}}\tau }|{\psi _{P}}\rangle ||}},}$

or, alternatively, simulate an isentropic evolution using a time-dependent Hamiltonian, which interpolates between the Hamiltonian ${\displaystyle H_{1}}$, which has the product state ${\displaystyle |{\psi _{P}}\rangle }$ as its ground state, and the Hamiltonian ${\displaystyle {\tilde {H}}}$; 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 ${\displaystyle |{\psi _{0}}\rangle }$ towards ${\displaystyle |{\psi _{T}}\rangle }$ using the Hamiltonian ${\displaystyle H_{n}}$:

${\displaystyle |{\psi _{T}}\rangle =e^{-iH_{n}T}|{\psi _{0}}\rangle }$

Errors coming from the Suzuki–Trotter expansion

In the case of a Trotter approximation of ${\displaystyle {\it {p^{th}}}}$ order, the error is of order ${\displaystyle {\delta }^{p+1}}$. Taking into account ${\displaystyle n={\frac {T}{\delta }}}$ steps, the error after the time T is:

${\displaystyle \epsilon ={\frac {T}{\delta }}\delta ^{p+1}=T\delta ^{p}}$

The unapproximated state ${\displaystyle |{{\tilde {\psi }}_{Tr}}\rangle }$ is:

${\displaystyle |{{\tilde {\psi }}_{Tr}}\rangle ={\sqrt {1-{\epsilon }^{2}}}|{\psi _{Tr}}\rangle +{\epsilon }|{\psi _{Tr}^{\bot }}\rangle }$

where ${\displaystyle |{\psi _{Tr}}\rangle }$ is the state kept after the Trotter expansion and ${\displaystyle |{\psi _{Tr}^{\bot }}\rangle }$ accounts for the part that is neglected when doing the expansion.

The total error scales with time ${\displaystyle {\it {T}}}$ as:

${\displaystyle \epsilon ({\it {T}})=1-|\langle {\tilde {\psi _{Tr}}}|{\psi _{Tr}}\rangle |^{2}=1-1+\epsilon ^{2}=\epsilon ^{2}}$

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:

${\displaystyle \epsilon ({\it {D}})=1-\prod \limits _{n=1}^{N-1}(1-\epsilon _{n})}$

where ${\displaystyle \epsilon _{n}=\sum \limits _{\alpha =\chi _{c}}^{\chi }(\lambda _{\alpha }^{[n]})^{2}}$ is the sum of all the discarded eigenvalues of the reduced density matrix, at the bond ${\displaystyle {\it {n}}}$. The state ${\displaystyle |{\psi }\rangle }$ is, at a given bond ${\displaystyle {\it {n}}}$, described by the Schmidt decomposition:

${\displaystyle |{\psi }\rangle ={\sqrt {1-\epsilon _{n}}}|{\psi _{D}}\rangle +{\sqrt {\epsilon _{n}}}|{\psi _{D}^{\bot }}\rangle }$

where

${\displaystyle |{\psi _{D}}\rangle ={\frac {1}{\sqrt {1-\epsilon _{n}}}}\sum \limits _{{{\alpha }_{n}}=1}^{{\chi }_{c}}\lambda _{{\alpha }_{n}}^{[n]}|{\Phi _{\alpha _{n}}^{[1..n]}}\rangle |{\Phi _{\alpha _{n}}^{[n+1..N]}}\rangle }$

is the state kept after the truncation and

${\displaystyle |{\psi _{D}^{\bot }}\rangle ={\frac {1}{\sqrt {\epsilon _{n}}}}\sum \limits _{{{\alpha }_{n}}={\chi }_{c}}^{\chi }\lambda _{{\alpha }_{n}}^{[n]}|{\Phi _{\alpha _{n}}^{[1..n]}}\rangle |{\Phi _{\alpha _{n}}^{[n+1..N]}}\rangle }$

is the state formed by the eigenfunctions corresponding to the smallest, irrelevant Schmidt coefficients, which are neglected. Now, ${\displaystyle \langle \psi _{D}^{\bot }|\psi _{D}\rangle =0}$ 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:

${\displaystyle \epsilon _{n}=1-|\langle {\psi }|\psi _{D}\rangle |^{2}=\sum \limits _{\alpha =\chi _{c}}^{\chi }(\lambda _{\alpha }^{[n]})^{2}}$

After moving to the next bond, the state is, similarly:

${\displaystyle |{\psi _{D}}\rangle ={\sqrt {1-\epsilon _{n+1}}}|{{{\psi }'}_{D}}\rangle +{\sqrt {\epsilon _{n+1}}}|{{\psi '}_{D}^{\bot }}\rangle }$

The error, after the second truncation, is:

${\displaystyle \epsilon =1-|\langle {\psi }|\psi '_{D}\rangle |^{2}=1-(1-\epsilon _{n+1})|\langle {\psi }|\psi _{D}\rangle |^{2}=1-(1-\epsilon _{n+1})(1-\epsilon _{n})}$

and so on, as we move from bond to bond.

The second error source enfolded in the decomposition ${\displaystyle {\it {D}}}$ is more subtle and requires a little bit of calculation.

As we calculated before, the normalization constant after making the truncation at bond ${\displaystyle {\it {l}}}$ ${\displaystyle ([1..l]:[l+1..N])}$ is:

${\displaystyle R={\sum \limits _{{\alpha _{l}}=1}^{{\chi }_{c}}{|\lambda _{{\alpha }_{l}}^{[l]}|}^{2}}={1-\epsilon _{l}}}$

Now let us go to the bond ${\displaystyle {\it {l}}-1}$ and calculate the norm of the right-hand Schmidt vectors ${\displaystyle ||{\Phi _{\alpha _{l-1}}^{[l-1..N]}}||}$; taking into account the full Schmidt dimension, the norm is:

${\displaystyle n_{1}=1=\sum \limits _{\alpha _{l}=1}^{\chi _{c}}(c_{\alpha _{l-1}\alpha _{l}})^{2}(\lambda _{\alpha _{l}}^{[l]})^{2}+\sum \limits _{\alpha _{l}=\chi _{c}}^{\chi }(c_{\alpha _{l-1}\alpha _{l}})^{2}(\lambda _{\alpha _{l}}^{[l]})^{2}=S_{1}+S_{2}}$,

where ${\displaystyle (c_{\alpha _{l-1}\alpha _{l}})^{2}=\sum \limits _{i_{l}=1}^{d}(\Gamma _{\alpha _{l-1}\alpha _{l}}^{[l]i_{l}})^{*}\Gamma _{\alpha _{l-1}\alpha _{l}}^{[l]i_{l}}}$.

Taking into account the truncated space, the norm is:

${\displaystyle n_{2}=\sum \limits _{\alpha _{l}=1}^{\chi _{c}}(c_{\alpha _{l-1}\alpha _{l}})^{2}\cdot ({\lambda '}_{\alpha _{l}}^{[l]})^{2}=\sum \limits _{\alpha _{l}=1}^{\chi _{c}}(c_{\alpha _{l-1}\alpha _{l}})^{2}{\frac {(\lambda _{\alpha _{l}}^{[l]})^{2}}{R}}={\frac {S_{1}}{R}}}$

Taking the difference, ${\displaystyle \epsilon =n_{2}-n_{1}=n_{2}-1}$, we get:

${\displaystyle \epsilon ={\frac {S_{1}}{R}}-1\leq {\frac {1-R}{R}}={\frac {\epsilon _{l}}{1-\epsilon _{l}}}{\rightarrow }0\ \ as\ \ {\epsilon _{l}{\rightarrow }{0}}}$

Hence, when constructing the reduced density matrix, the trace of the matrix is multiplied by the factor:

${\displaystyle |\langle {\psi _{D}}|\psi _{D}\rangle |^{2}=1-{\frac {\epsilon _{l}}{1-\epsilon _{l}}}={\frac {1-2\epsilon _{l}}{1-\epsilon _{l}}}}$

The total truncation error

The total truncation error, considering both sources, is upper bounded by:

${\displaystyle \epsilon ({D})=1-\prod \limits _{n=1}^{N-1}(1-\epsilon _{n})\prod \limits _{n=1}^{N-1}{\frac {1-2\epsilon _{n}}{1-\epsilon _{n}}}=1-\prod \limits _{n=1}^{N-1}(1-2\epsilon _{n})}$

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.