© The Author(s) 2020
S.-J. Ran et al.Tensor Network ContractionsLecture Notes in Physics964https://doi.org/10.1007/978-3-030-34489-4_6

6. Quantum Entanglement Simulation Inspired by Tensor Network

Shi-Ju Ran1 , Emanuele Tirrito2, Cheng Peng3, Xi Chen4, Luca Tagliacozzo5, Gang Su6 and Maciej Lewenstein2
(1)
Department of Physics, Capital Normal University, Beijing, China
(2)
Quantum Optics Theory, Institute of Photonic Sciences, Castelldefels, Spain
(3)
Stanford Institute for Materials and Energy Sciences, SLAC and Stanford University, Menlo Park, CA, USA
(4)
School of Physical Sciences, University of Chinese Academy of Science, Beijing, China
(5)
Department of Quantum Physics and Astrophysics, University of Barcelona, Barcelona, Spain
(6)
Kavli Institute for Theoretical Sciences, University of Chinese Academy of Science, Beijing, China
 

Abstract

This chapter is focused on the quantum entanglement simulation approach (Ran et al., Phys Rev B 96:155120, 2017). The idea is to use few-body models embedded in the “entanglement bath” to mimic the properties of large and infinite-size quantum systems. The few-body models are dubbed as quantum entanglement simulators. Generally speaking, the QES approach consists of three stages: first, determine the effective interactions that give the infinite boundary condition (Ran et al., Phys Rev B 96:155120, 2017; Phien et al., Phys Rev B 86:245107, 2012) by the MPS/PEPS methods, such as iDMRG in one dimension or zero-loop scheme in two and three dimensions; second, construct the simulators by surrounding a finite-size cluster of the targeted model with the effective interactions; third, simulate the properties of the quantum entanglement simulator by the finite-size algorithms or quantum hardware, and extract the properties of the targeted model within the bulk.

6.1 Motivation and General Ideas

An impression one may have for the TN approaches of the quantum lattice models is that the algorithm (i.e., how to contract and/or truncate) will dramatically change when considering different models or lattices. This motivates us to look for more unified approaches. Considering that a huge number of problems can be transformed to TN contractions, one general question we may ask is: how can we reduce a non-deterministic polynomial hard TN contraction problem approximately to an effective one that can be computed exactly and efficiently by classical computers? We shall put some principles while considering this question: the effective problem should be as simple as possible, containing as few parameters to solve as possible. We would like to coin this principle for TN contractions as the ab initio optimization principle (AOP) of TN [1]. The term “ab inito” is taken here to pay respect to the famous ab inito principle approaches in condensed matter physics and quantum chemistry (see several recent reviews in [2–4]). Here, “ab inito” means to think from the very beginning, with least prior knowledge of or assumptions to the problems.

One progress achieved in the spirit of AOP is the TRD introduced in Sect. 5.​4. Considering the TN on an infinite square lattice, its contraction is reduced to a set of self-consistent eigenvalue equations that can be efficiently solved by classical computers. The variational parameters are just two tensors. One advantage of TRD is that it connects the TN algorithms (iDMRG , iTEBD , CTMRG ), which are previously considered to be quite different, in a unified picture.

Another progress made in the AOP spirit is called QES for simulating infinite-size physical models [1, 5, 6]. It is less dependent on the specific models; it also provides a natural way for designing quantum simulators and for hybridized-quantum-classical simulations of many-body systems. Hopefully in the future when people are able to readily realize the designed Hamiltonians on artificial quantum platforms, QES will enable us to design the Hamiltonians that will realize quantum many-body phenomena.

6.2 Simulating One-Dimensional Quantum Lattice Models

Let us firstly take the ground-state simulation of the infinite-size 1D quantum system as an example. The Hamiltonian is the summation of two-body nearest-neighbor terms, which reads $$\hat {H}_{Inf} = \sum _{n} \hat {H}_{n,n+1}$$ . The translational invariance is imposed. The first step is to choose a supercell (e.g., a finite part of the chain with $$\tilde {N}$$ sites). Then the Hamiltonian of the supercell is $$\hat {H}_{B} = \sum _{n=1}^{\tilde {N}} \hat {H}_{n,n+1}$$ , and the Hamiltonian connecting the supercell to the rest part is $$\hat {H}_{\partial } = \hat {H}_{n',n'+1}$$ (note the interactions are nearest neighbor).

Define the operator $$\hat {F}^{\partial }$$ as
$$\displaystyle \begin{aligned} \begin{array}{rcl} \hat{F}_{\partial} = \hat{I} - \tau \hat{H}_{\partial}, {} \end{array} \end{aligned} $$
(6.1)
with τ the Trotter-Suzuki step. This definition is to construct the Trotter-Suzuki decomposition [7, 8]. Instead of using the exponential form $$e^{-\tau \hat {H}}$$ , we equivalently chose to shift $$\hat {H}_{\partial }$$ for algorithmic consideration. The errors of these two ways concerning the ground state are at the same level ( $$\mathscr {O}(\tau ^2)$$ ). Introduce an ancillary index a and rewrite $$\hat {F}_{\partial }$$ as a sum of operators as
$$\displaystyle \begin{aligned} \begin{array}{rcl} \hat{F}_{\partial} = \sum_a \hat{F}_{L}(s)_a \hat{F}_{R}(s')_a, {} \end{array} \end{aligned} $$
(6.2)
where $$\hat {F}_{L}(s)_a$$ and $$\hat {F}_{R}(s')_a$$ are two sets of one-body operators (labeled by a) acting on the left and right one of the two spins (s and s′) associated with $$\hat {H}_{\partial }$$ , respectively (Fig. 6.1). Equation (6.2) can be easily achieved directly from the Hamiltonian or using eigenvalue decomposition. For example, for the Heisenberg interaction with $$\hat {H}_{\partial } = \sum _{\alpha =x,y,z} J_{\alpha } \hat {S}^{\alpha }(s)\hat {S}^{\alpha }(s')$$ with $$\hat {S}^{\alpha }(s)$$ the spin operators. We have $$\hat {F}_{\partial } = \sum _{a = 0,x,y,z} \tilde {J}_{a} \hat {S}^{a}(s) \hat {S}^{a}(s')$$ with $$\hat {S}^{0} = I$$ , $$\tilde {J}_0 = 1$$ , $$\tilde {J}_{a} = -\tau J_{a}$$ (α = x, y, z), hence we can define $$\hat {F}_{L}(s)_a = \sqrt {|\tilde {J}_{a}|} \hat {S}^{\alpha }(s)$$ and $$\hat {F}_{R}(s)_a = \text{sign}(\tilde {J}_{a}) \sqrt {|\tilde {J}_{a}|} \hat {S}^{\alpha }(s')$$ .
../images/489509_1_En_6_Chapter/489509_1_En_6_Fig1_HTML.png
Fig. 6.1

Graphical representations of Eqs. (6.1)–(6.4)

Construct the operator $$\hat {\mathscr {F}}(S)_{a'a}$$ , with $$S=(s_1,\cdots ,s_{\tilde {N}})$$ representing the physical spins inside the supercell, as
$$\displaystyle \begin{aligned} \begin{array}{rcl} \hat{\mathscr{F}}(S)_{a'a}= \hat{F}_{R}(s_1)_{a'}^{\dagger} \tilde{H}_B \hat{F}_{L}(s_{\tilde{N}})_a, {} \end{array} \end{aligned} $$
(6.3)
with $$\tilde {H}_B = \hat {I} - \tau \hat {H}_B$$ . $$\hat {F}_{R}(s_1)_{a'}^{\dagger }$$ and $$\hat {F}_{L}(s_{\tilde {N}})_{a}$$ act on the first and last sites of the supercell, respectively. One can see that $$\hat {\mathscr {F}}(S)_{a'a}$$ represents a set of operators labeled by two indexes (a′ and a) that act on the supercell.
In the language of TN , the coefficients of $$\hat {\mathscr {F}}(S)_{a'a}$$ in the local basis ( $$|S\rangle = |s_1\rangle \cdots |s_{\tilde {N}}\rangle $$ ) is a fourth-order cell tensor (Fig. 6.1) as
$$\displaystyle \begin{aligned} \begin{array}{rcl} T_{S'a'Sa} = \langle S'| \hat{\mathscr{F}}(S)_{a'a} | S\rangle. {} \end{array} \end{aligned} $$
(6.4)
On the left-hand side, the order of the indexes are arranged to be consistent with the definition in the TN algorithm introduced in the precious sections. T is the cell tensor, whose infinite copies form the TN of the imaginary-time evolution up to the first Trotter-Suzuki order. One may consider the second Trotter-Suzuki order by defining $$\hat {\mathscr {F}}(S)_{a'a}$$ as $$\hat {\mathscr {F}}(S)_{a'a} = (\hat {I} - \tau \hat {H}_B /2) \hat {F}_{R}(s_1)_{a'}^{\dagger } \hat {F}_{L}(s_{\tilde {N}})_a (\hat {I} - \tau \hat {H}_B /2)$$ . With the cell tensor T, the ground-state properties can be solved using the TN algorithms (e.g., TRD) introduced above. The ground state is given by the MPS given by Eq. (5.​59).
Let us consider one of the eigenvalue equations [also see Eq. (5.​57)] in the TRD
$$\displaystyle \begin{aligned} \begin{array}{rcl} \mathscr{H}_{S'b_1^{\prime}b_2^{\prime},Sb_1b_2} = \sum_{aa'} T_{S'a'Sa} v^{L}_{a'b_1^{\prime}b_1} v^R_{ab_2^{\prime}b_2}. {} \end{array} \end{aligned} $$
(6.5)
We define a new Hamiltonian $$\hat {\mathscr {H}}$$ by using $$\mathscr {H}_{S'b_1^{\prime }b_2^{\prime },Sb_1b_2}$$ as the coefficients
$$\displaystyle \begin{aligned} \begin{array}{rcl} \hat{\mathscr{H}} =\sum_{SS'} \sum_{b_1b_2b_1^{\prime}b_2^{\prime}} \mathscr{H}_{S'b_1^{\prime}b_2^{\prime},Sb_1b_2} |S'b_1^{\prime}b_2^{\prime}\rangle\langle Sb_1b_2|. {} \end{array} \end{aligned} $$
(6.6)
$$\hat {\mathscr {H}}$$ is the effective Hamiltonian in iDMRG [9–11] or the methods which represent the RG of Hilbert space by MPS [12, 13]. The indexes {b} are considered as virtual spins with basis {|b〉}. The virtual spins are called the entanglement bath sites in the QES .
By substituting with the cell tensor T [Eqs. (6.3) and (6.4)] inside the above equation, we have
$$\displaystyle \begin{aligned} \begin{array}{rcl} \hat{\mathscr{H}} = \hat{\mathscr{H}}_L \tilde{H}_B \hat{\mathscr{H}}_R, {} \end{array} \end{aligned} $$
(6.7)
where the Hamiltonians $$\hat {\mathscr {H}}_L$$ and $$\hat {\mathscr {H}}_R$$ locate on the boundaries of $$\hat {\mathscr {H}}$$ , whose coefficients satisfy
$$\displaystyle \begin{aligned} \begin{array}{rcl} \langle b_1^{\prime}s_{1}^{\prime}| \hat{\mathscr{H}}_L |b_1 s_{1}\rangle &\displaystyle =&\displaystyle \sum_{a} v^{L}_{ab_1^{\prime}b_1} \langle s_1^{\prime} |\hat{F}_{R}(s_1)_{a}^{\dagger} | s_1\rangle, \\ \langle s_{\tilde{N}}^{\prime} b_2^{\prime}| \hat{\mathscr{H}}_R |s_{\tilde{N}} b_2\rangle &\displaystyle =&\displaystyle \sum_a \langle s_{\tilde{N}}^{\prime}| \hat{F}_{L}(s_{\tilde{N}})_a^{\dagger} | s_{\tilde{N}}\rangle v^{R}_{ab_2b_2^{\prime}}. {} \end{array} \end{aligned} $$
(6.8)
$$\hat {\mathscr {H}}_L$$ and $$\hat {\mathscr {H}}_R$$ are just two-body Hamiltonians, of which each acts on the bath site and its neighboring physical site on the boundary of the bulk; they define the infinite boundary condition for simulating the time evolution of 1D quantum systems [14].
$$\hat {\mathscr {H}}_L$$ and $$\hat {\mathscr {H}}_R$$ can also be written in a shifted form as
$$\displaystyle \begin{aligned} \hat{\mathscr{H}}_{L(R)} = I-\tau \hat{H}_{L(R)}. {} \end{aligned} $$
(6.9)
This is because the tensor v L (and also v R) satisfies a special form [15] as
$$\displaystyle \begin{aligned} v^{L}_{0,bb'} = I_{bb'} - \tau Q_{bb'}, \end{aligned} $$
(6.10)
$$\displaystyle \begin{aligned} v^{L}_{a,bb'} = \tau^{a/2} (R^a)_{bb'} \ \ (a > 0), \end{aligned} $$
(6.11)
with Q and R two Hermitian matrices independent on τ. In other words, the MPS formed by infinite copies of v L or v R is a continuous MPS [16], which is known as the temporal MPS [17]. Therefore, $$\hat {H}_{L(R)}$$ is independent on τ, called the physical-bath Hamiltonian. Then $$\hat {\mathscr {H}}$$ can be written as the shift of a few-body Hamiltonian as $$\hat {\mathscr {H}} = I- \tau \hat {H}_{FB}$$ , where $$\hat {H}_{FB}$$ has the standard summation form as
$$\displaystyle \begin{aligned} \begin{array}{rcl} \hat{H}_{FB} = \hat{H}_L + \sum_{n=1}^{L} \hat{H}_{n,n+1} + \hat{H}_R. {} \end{array} \end{aligned} $$
(6.12)
For $$\hat {H}_{L}$$ and $$\hat {H}_{R}$$ with the bath dimension χ, the coefficient matrix of $$\hat {H}_{L(R)}$$ is (2χ × 2χ). Then $$\hat {H}_{L(R)}$$ can be generally expanded by $$\hat {S}^{\alpha _1} \otimes \hat {\mathscr {S}}^{\alpha _2}$$ with $$\{\hat {\mathscr {S}}\}$$ the generators of the SU(χ) group, and define the magnetic field and coupling constants associated to the entanglement bath
$$\displaystyle \begin{aligned} \begin{array}{rcl} \hat{H}_{L(R)} = \sum_{\alpha_1,\alpha_2} J^{\alpha_1\alpha_2}_{L(R)} \hat{\mathscr{S}}^{\alpha_1} \otimes \hat{S}^{\alpha_2}, {} \end{array} \end{aligned} $$
(6.13)
with $$\mathscr {S}$$ denoting the SU(χ) spin operators and $$\hat {S}$$ the operators of the physical spin (with the identity included).
Let us take the bond dimension χ = 2 as an example, and $$\hat {H}_{L(R)}$$ just gives the Hamiltonian between two spin-1∕2’s. Thus, it can be expanded by the spin (or Pauli) operators $$\hat {S}^{\alpha _1} \otimes \hat {S}^{\alpha _2}$$ as
$$\displaystyle \begin{aligned} \begin{array}{rcl} \hat{H}_{L(R)} = \sum_{\alpha_1,\alpha_2=0}^3 J^{\alpha_1\alpha_2}_{L(R)} \hat{S}^{\alpha_1} \otimes \hat{S}^{\alpha_2}, {} \end{array} \end{aligned} $$
(6.14)
where the spin-1∕2 operators are labeled as $$\hat {S}^0 = I$$ , $$\hat {S}^1 = \hat {S}^x$$ , $$\hat {S}^2 = \hat {S}^y$$ , and $$\hat {S}^3 = \hat {S}^z$$ . Then with α 1 ≠ 0 and α 2 ≠ 0, we have $$J^{\alpha _1\alpha _2}_{L(R)}$$ as the coupling constants, and $$J^{\alpha _1 0}_{L(R)}$$ and $$J^{0 \alpha _2}_{L(R)}$$ the magnetic fields on the first and second sites, respectively. $$J^{0 0}_{L(R)}$$ only provides a constant shift of the Hamiltonian which does not change the eigenstates.
As an example, we show the $$\hat {H}_{L}$$ and $$\hat {H}_{R}$$ for the infinite quantum Ising chain in a transverse field [18]. The original Hamiltonian reads $$\hat {H}_{Ising} = \sum _n \hat {S}^z_n \hat {S}^z_{n+1} -\alpha \sum _{n} \hat {S}^x_n$$ , and $$\hat {H}_{L}$$ and $$\hat {H}_{R}$$ satisfies
$$\displaystyle \begin{aligned} \begin{array}{rcl} \begin{aligned} &\displaystyle \hat{H}_{L} = J^L_{xz} \hat{S}^x_1 \hat{S}^z_2 + J^L_{zz} \hat{S}^z_1 \hat{S}^z_2 - h^L_{x} \hat{S}^x_1 - h^L_{z} \hat{S}^z_1 - \tilde{h}^L_{x} \hat{S}^x_2,\\ &\displaystyle \hat{H}_{R} = J^R_{zx} \hat{S}^z_{N-1} \hat{S}^x_N + J^R_{zz} \hat{S}^z_{N-1} \hat{S}^z_{N} - h^R_{x} \hat{S}^x_{N} - h^R_{z} \hat{S}^z_{N} - \tilde{h}^R_{x} \hat{S}^x_{N-1}. \end{aligned} {} \end{array} \end{aligned} $$
(6.15)
The coupling constants and magnetic fields depend on the transverse field α, as shown in Fig. 6.2. The calculation shows that except the Ising interactions and the transverse field that originally appear in the infinite model, the $$\hat {S}^x \hat {S}^z$$ coupling and a vertical field emerge in $$\hat {H}_{L}$$ and $$\hat {H}_{R}$$ . This is interesting, because the $$\hat {S}^x \hat {S}^z$$ interaction is the stabilizer on the open boundaries of the cluster state, a highly entangled state that has been widely used in quantum information sciences [19, 20]. More relations with the cluster state are to be further explored.
../images/489509_1_En_6_Chapter/489509_1_En_6_Fig2_HTML.png
Fig. 6.2

The α-dependence [18] of the coupling constants (left) and magnetic fields (right) of the few-body Hamiltonians (Eq. (6.15)). Reused from [18] with permission

The physical information of the infinite-size model can be extracted from the ground state of $$\hat {H}_{FB}$$ (denoted by |Ψ(Sb 1b 2)〉) by tracing over the entanglement-bath degrees of freedom. To this aim, we calculate the reduced density matrix of the bulk as
$$\displaystyle \begin{aligned} \begin{array}{rcl} \hat{\rho}(S) = \text{Tr}_{b_1b_2} |\varPsi(Sb_1b_2) \rangle \langle \varPsi(Sb_1b_2)|. {} \end{array} \end{aligned} $$
(6.16)
Note $$|\varPsi (Sb_1b_2)\rangle = \sum _{Sb_1b_2} \varPsi _{Sb_1b_2} |Sb_1b_2\rangle $$ with $$\varPsi _{Sb_1b_2}$$ the eigenvector of Eq. (6.5) or (5.​57). It is easy to see that $$\varPsi _{Sb_1b_2}$$ is the central tensor in the central-orthogonal MPS (Eq. (5.​59)), thus the $$\hat {\rho }(S)$$ is actually the reduced density matrix of the MPS. Since the MPS optimally gives the ground state of the infinite model, therefore, $$\hat {\rho }(S)$$ of the few-body ground state optimally gives the reduced density matrix of the original model.

In Eq. (6.12), the summation of the physical interactions is within the supercell that we choose to construct the cell tensor. To improve the accuracy to, e.g., capture longer correlations inside the bulk, one just needs to increase the supercell in $$\hat {H}_{FB}$$ . In other words, $$\hat {H}_L$$ and $$\hat {H}_R$$ are obtained by TRD from the supercell of a tolerable size $$\tilde {N}$$ , and $$\hat {H}_{FB}$$ is constructed with a larger bulk as $$\hat {H}_{FB} = \hat {H}_L + \sum _{n=1}^{\tilde {N}'} \hat {H}_{n,n+1} + \hat {H}_R$$ with $$\tilde {N}'>\tilde {N}$$ . Though $$\hat {H}_{FB}$$ becomes more expensive to solve, we have many well-established finite-size algorithms to compute its dominant eigenvector. We will show below that this way is extremely useful in higher dimensions.

6.3 Simulating Higher-Dimensional Quantum Systems

For (D  >  1)-dimensional quantum systems on, e.g., square lattice, one can use different update schemes to calculate the ground state. Here, we explain an alternative way by generalizing the above 1D simulation to higher dimensions [5]. The idea is to optimize the physical-bath Hamiltonians by the zero-loop approximation (simple update, see Sect. 5.​3), e.g., iDMRG on tree lattices [21, 22], and then construct the few-body Hamiltonian $$\hat {H}_{FB}$$ with larger bulks. The loops inside the bulk will be fully considered when solving the ground state of $$\hat {H}_{FB}$$ , thus the precision will be significantly improved compared with the zero-loop approximation.

The procedures are similar to those for 1D models. The first step is to contract the cell tensor, so that the ground-state simulation is transformed to a TN contraction problem. We choose the two sites connected by a parallel bond as the supercell, and construct the cell tensor that parametrizes the eigenvalue equations. The bulk interaction is simply the coupling between these two spins, i.e., $$\hat {H}_B = \hat {H}_{i,j}$$ , and the interaction between two neighboring supercells is the same, i.e., $$\hat {H}_{\partial } = \hat {H}_{i,j}$$ . By shifting $$\hat {H}_{\partial }$$ , we define $$\hat {F}_{\partial } = I-\tau \hat {H}_{\partial }$$ and decompose it as
$$\displaystyle \begin{aligned} \begin{array}{rcl} \hat{F}_{\partial} = \sum_{a} \hat{F}_{L}(s)_a \otimes \hat{F}_{R}(s')_a. {} \end{array} \end{aligned} $$
(6.17)
$$\hat {F}_{L}(s)_a$$ and $$\hat {F}_{R}(s')_a$$ are two sets of operators labeled by a that act on the two spins (s and s′) in the supercell, respectively (see the texts below Eq. (6.2) for more detail).
Define a set of operators by the product of the (shifted) bulk Hamiltonian with $$\hat {F}_{L}(s)_a$$ and $$\hat {F}_{R}(s)_a$$ (Fig. 6.3) as
$$\displaystyle \begin{aligned} \hat{\mathscr{F}}(S)_{a_1a_2a_3a_4}= \hat{F}_{R}(s)_{a_1} \hat{F}_{R}(s)_{a_2} \hat{F}_{L}(s')_{a_3} \hat{F}_{L}(s')_{a_4} \tilde{H}^B, {} \end{aligned} $$
(6.18)
with S = (s, s′) and $$\tilde {H}^B = I- \tau \hat {H}^B$$ . The cell tensor that defines the TN is given by the coefficients of $$\hat {\mathscr {F}}(S)_{a_1a_2a_3a_4}$$ as
$$\displaystyle \begin{aligned} \begin{array}{rcl} T_{S'Sa_1a_2a_3a_4} = \langle S'| \hat{\mathscr{F}}(S)_{a_1a_2a_3a_4} |S\rangle. {} \end{array} \end{aligned} $$
(6.19)
One can see that T has six bonds, of which two (S and S′) are physical and four (a 1, a 2, a 3, and a 4) are non-physical. For comparison, the tensor in the 1D quantum case has four bonds, where two are physical and two are non-physical [see Eq. (6.4)]. As discussed above in Sect. 4.​2, the ground-state simulation becomes the contraction of a cubic TN formed by infinite copies of T. Each layer of the cubic TN gives the operator $$\hat {\rho }(\tau ) = I-\tau \hat {H}$$ , which is a PEPO defined on a square lattice. Infinite layers of the PEPO $$\lim _{K \to \infty } \hat {\rho }(\tau )^K$$ give the cubic TN .
../images/489509_1_En_6_Chapter/489509_1_En_6_Fig3_HTML.png
Fig. 6.3

Graphical representation of the cell tensor for 2D quantum systems (Eq. (6.18))

The next step is to solve the SEEs of the zero-loop approximation. For the same model defined on the loopless Bethe lattice, the 3D TN is formed by infinite layers of PEPO $$\hat {\rho }_{Bethe}(\tau )$$ that is defined on the Bethe lattice. The cell tensor is defined exactly in the same way as Eq. (6.19). With the Bethe approximation, there are five variational tensors, which are Ψ (central tensor) and v [x] (x = 1, 2, 3, 4, boundary tensors). Meanwhile, we have five self-consistent equations that encodes the 3D TN $$\lim _{K \to \infty } \hat {\rho }_{Bethe}(\tau )^K$$ , which are given by five matrices as
$$\displaystyle \begin{aligned} \begin{array}{rcl} \mathscr{H}_{S'b_1^{\prime}b_2^{\prime}b_3^{\prime}b_4^{\prime},Sb_1b_2b_3b_4} = \sum_{a_1a_2a_3a_4} T_{S'Sa_1a_2a_3a_4} v^{[1]}_{a_1b_1 b_1^{\prime}} v^{[2]}_{a_2b_2 b_2^{\prime}} v^{[3]}_{a_3b_3 b_3^{\prime}} v^{[4]}_{a_4b_4 b_4^{\prime}}, {}\\ \end{array} \end{aligned} $$
(6.20)
$$\displaystyle \begin{aligned} \begin{array}{rcl} M^{[1]}_{a_1b_1b_1^{\prime},a_3b_3b_3^{\prime}} = \sum_{S'Sa_2a_4b_2b_2^{\prime}b_4b_4^{\prime}} T_{S'Sa_1a_2a_3a_4} A^{[1]\ast}_{S'b_1^{\prime}b_2^{\prime}b_3^{\prime}b_4^{\prime}} v^{[2]}_{a_2b_2 b_2^{\prime}} A^{[1]}_{Sb_1b_2b_3b_4} v^{[4]}_{a_4b_4 b_4^{\prime}},{}\\ \end{array} \end{aligned} $$
(6.21)
$$\displaystyle \begin{aligned} \begin{array}{rcl} M^{[2]}_{a_2b_2b_2^{\prime},a_4b_4b_4^{\prime}} = \sum_{S'Sa_1a_3b_1b_1^{\prime}b_3b_3^{\prime}} T_{S'Sa_1a_2a_3a_4} A^{[2]\ast}_{S'b_1^{\prime}b_2^{\prime}b_3^{\prime}b_4^{\prime}} v^{[1]}_{a_1b_1 b_1^{\prime}} A^{[2]}_{Sb_1b_2b_3b_4} v^{[3]}_{a_3b_3 b_3^{\prime}},{}\\ \end{array} \end{aligned} $$
(6.22)
$$\displaystyle \begin{aligned} \begin{array}{rcl} M^{[3]}_{a_1b_1b_1^{\prime},a_3b_3b_3^{\prime}} = \sum_{S'Sa_2a_4b_2b_2^{\prime}b_4b_4^{\prime}} T_{S'Sa_1a_2a_3a_4} A^{[3]\ast}_{S'b_1^{\prime}b_2^{\prime}b_3^{\prime}b_4^{\prime}} v^{[2]}_{a_2b_2 b_2^{\prime}} A^{[3]}_{Sb_1b_2b_3b_4} v^{[4]}_{a_4b_4 b_4^{\prime}},{}\\ \end{array} \end{aligned} $$
(6.23)
$$\displaystyle \begin{aligned} \begin{array}{rcl} M^{[4]}_{a_2b_2b_2^{\prime},a_4b_4b_4^{\prime}} = \sum_{S'Sa_1a_3b_1b_1^{\prime}b_3b_3^{\prime}} T_{S'Sa_1a_2a_3a_4} A^{[4]\ast}_{S'b_1^{\prime}b_2^{\prime}b_3^{\prime}b_4^{\prime}} v^{[1]}_{a_1b_1 b_1^{\prime}} A^{[4]}_{Sb_1b_2b_3b_4} v^{[3]}_{a_3b_3 b_3^{\prime}}.{}\\ \end{array} \end{aligned} $$
(6.24)
Equations (6.20) and (6.22) are illustrated in Fig. 6.4 as two examples. A [x] is an isometry obtained by the QR decomposition (or SVD ) of the central tensor Ψ referring to the x-th virtual bond b x. For example, for x = 2, we have (Fig. 6.4)
$$\displaystyle \begin{aligned} \begin{array}{rcl} \varPsi_{Sb_1b_2b_3b_4} = \sum_{b} A^{[2]}_{Sb_1bb_3b_4} R^{[2]}_{b b_2}. {} \end{array} \end{aligned} $$
(6.25)
../images/489509_1_En_6_Chapter/489509_1_En_6_Fig4_HTML.png
Fig. 6.4

The left figure is the graphic representations of $$\mathscr {H}_{S'b_1^{\prime }b_2^{\prime }b_3^{\prime }b_4^{\prime },Sb_1b_2b_3b_4}$$ in Eq. (6.20), and we take Eq. (6.22) from the self-consistent equations as an example shown in the middle. The QR decomposition in Eq. (6.25) is shown in the right figure, where the arrows indicate the direction of orthogonality of A [3] in Eq. (6.26)

A [2] is orthogonal, satisfying
$$\displaystyle \begin{aligned} \begin{array}{rcl} \sum_{Sb_1 b_3 b_4} A^{[2] \ast}_{Sb_1bb_3b_4} A^{[2]}_{Sb_1b'b_3b_4} = I_{bb'}. {} \end{array} \end{aligned} $$
(6.26)

The self-consistent equations can be solved recursively. By solving the leading eigenvector of $$\mathscr {H}$$ given by Eq. (6.20), we update the central tensor Ψ. Then according to Eq. (6.25), we decompose Ψ to obtain A [x], then update M [x] in Eqs. (6.21)–(6.24), and update each v [x] by M [x]v [x]. Repeat this process until all the five variational tensors converge. The algorithm is the generalized DMRG based on infinite tree PEPS [21, 22]. Each boundary tensor can be understood as the infinite environment of a tree branch, thus the original model is actually approximated at this stage by that defined on an Bethe lattice. Note that when only looking at the tree locally (from one site and its nearest neighbors), it looks the same to the original lattice. Thus, the loss of information is mainly long range, i.e., from the destruction of loops.

We can have a deeper understanding of the Bethe approximation with the help of rank-1 decomposition explained in Sect. 5.​3. Equations (6.21)–(6.24) encode a Bethe TN, whose contraction is written as $$Z_{Bethe} = \langle \tilde {\varPhi }| \hat {\rho }_{Bethe}(\tau ) | \tilde {\varPhi } \rangle $$ with $$\hat {\rho }_{Bethe}(\tau )$$ the PEPO of the Bethe model and $$| \tilde {\varPhi } \rangle $$ a tree iPEPS (Fig. 6.5). To see this, let us start with the local contraction (Fig. 6.5a) as
$$\displaystyle \begin{aligned} \begin{array}{rcl} Z_{Bethe} = \sum \varPsi^{\ast}_{S'b_1^{\prime}b_2^{\prime}b_3^{\prime}b_4^{\prime}} \varPsi_{Sb_1b_2b_3b_4} T_{S'Sa_1a_2a_3a_4} v^{[1]}_{a_1b_1 b_1^{\prime}} v^{[2]}_{a_2b_2 b_2^{\prime}} v^{[3]}_{a_3b_3 b_3^{\prime}} v^{[4]}_{a_4b_4 b_4^{\prime}}.\qquad {} \end{array} \end{aligned} $$
(6.27)
Then, each v [x] can be replaced by M [x]v [x] because we are at the fixed point of the eigenvalue equations. By repeating this substitution in a similar way as the rank-1 decomposition in Sect. 5.​3.​3, we will have the TN for Z Bethe, which is maximized at the fixed point (Fig. 6.5b). With the constraint $$\langle \tilde {\varPhi }| \tilde {\varPhi } \rangle =1$$ satisfied, $$| \tilde {\varPhi } \rangle $$ is the ground state of $$\hat {\rho }_{Bethe}(\tau )$$ .
../images/489509_1_En_6_Chapter/489509_1_En_6_Fig5_HTML.png
Fig. 6.5

The left figure shows the local contraction that encodes the infinite TN for simulating the 2D ground state. By substituting with the self-consistent equations, the TN representing $$\tilde {Z} = \langle \tilde {\varPhi }| \hat {\rho }_{Bethe}(\tau ) | \tilde {\varPhi } \rangle $$ can be reconstructed, with $$\hat {\rho }_{Bethe}(\tau )$$ the tree PEPO of the Bethe model and $$| \tilde {\varPhi } \rangle $$ a PEPS

Now, we constrain the growth so that the TN covers the infinite square lattice. Inevitably, some v [x]s will gather at the same site. The tensor product of these v [x]s in fact gives the optimal rank-1 approximation of the “correct” full-rank tensor here (Sect. 5.​3.​3). Suppose that one uses the full-rank tensor to replace its rank-1 version (the tensor product of four v [x]’s), one will have the PEPO of $$I- \tau \hat {H}$$ (with H the Hamiltonian on square lattice), and the tree iPEPS becomes the iPEPS defined on the square lattice. Compared with the NCD scheme that employs rank-1 decomposition explicitly to solve TN contraction, one difference here for updating iPEPS is that the “correct” tensor to be decomposed by rank-1 decomposition contains the variational tensor, thus is in fact unknown before the equations are solved. For this reason, we cannot use rank-1 decomposition directly. Another difference is that the constraint, i.e., the normalization of the tree iPEPS, should be fulfilled. By utilizing the iDMRG algorithm with the tree iPEPS, the rank-1 tensor is obtained without knowing the “correct” tensor, and meanwhile the constraints are satisfied. The zero-loop approximation of the ground state is thus given by the tree iPEPS.

The few-body Hamiltonian is constructed in a larger cluster, so that the error brought by zero-loop approximation can be reduced. Similar to the 1D case, we embed a larger cluster in the middle of the entanglement bath. The few-body Hamiltonian (Fig. 6.6) is written as
$$\displaystyle \begin{aligned} \begin{array}{rcl} \hat{\mathscr{H}} = \prod_{\langle n \in cluster, \alpha \in bath \rangle} \hat{\mathscr{H}}_{\partial}(n,\alpha) \prod_{\langle i,j \rangle \in cluster} [I-\tau \hat{H}(s_i,s_j)]. {} \end{array} \end{aligned} $$
(6.28)
$$\hat {\mathscr {H}}_{\partial }(n,\alpha )$$ is defined as the physical-bath Hamiltonian between the α-th bath site and the neighboring n-th physical site, and it is obtained by the corresponding boundary tensor v [x(α)] and $$\hat {F}_{L(R)}(s_n)$$ (Fig. 6.6) as
$$\displaystyle \begin{aligned} \begin{array}{rcl} \langle b_{\alpha}^{\prime}s_n^{\prime}| \hat{\mathscr{H}}_{\partial}(n,\alpha) |b_{\alpha}s_n\rangle = \sum_{a} v^{[x(\alpha)]}_{ab_{\alpha}^{\prime} b_{\alpha}} \langle s_n^{\prime} |\hat{F}_{L(R)}(s_n)_{a} | s_n\rangle. {} \end{array} \end{aligned} $$
(6.29)
Here, $$\hat {F}_{L(R)}(s_n)_{a}$$ is the operator defined in Eq. (6.17), and $$v^{[x(\alpha )]}_{ab_{\alpha }^{\prime } b_{\alpha }}$$ are the solutions of the SEEs given in Eqs. (6.20)–(6.24).
../images/489509_1_En_6_Chapter/489509_1_En_6_Fig6_HTML.png
Fig. 6.6

The left figure shows the few-body Hamiltonian $$\hat {\mathscr {H}}$$ in Eq. (6.28). The middle one shows the physical-bath Hamiltonian $$\hat {\mathscr {H}}_{\partial }$$ that gives the interaction between the corresponding physical and bath site. The right one illustrates the state ansatz for the infinite system. Note that the boundary of the cluster should be surrounded by $$\hat {\mathscr {H}}_{\partial }$$ ’s, and each $$\hat {\mathscr {H}}_{\partial }$$ corresponds to an infinite tree brunch in the state ansatz. For simplicity, we only illustrate four of the $$\hat {\mathscr {H}}_{\partial }$$ s and the corresponding brunches

$$\hat {\mathscr {H}}$$ in Eq. (6.28) can also be rewritten as the shift of a few-body Hamiltonian $$\hat {H}_{FB}$$ , i.e., $$\hat {\mathscr {H}} = I-\tau \hat {H}_{FB}$$ . We have $$\hat {H}_{FB}$$ possessing the standard summation form as
$$\displaystyle \begin{aligned} \begin{array}{rcl} \hat{H}_{FB} = \sum_{\langle i,j \rangle \in cluster} \hat{H}(s_i,s_j) + \sum_{\langle n \in cluster,\alpha \in bath \rangle} \hat{H}_{PB}(n,\alpha), {} \end{array} \end{aligned} $$
(6.30)
with $$\hat {\mathscr {H}}_{\partial }(n,\alpha ) = I-\tau \hat {H}_{PB}(s_n,b_{\alpha })$$ . This equations gives a general form of the few-body Hamiltonian: the first term contains all the physical interactions inside the cluster, and the second contains all physical-bath interactions $$\hat {H}_{PB}(s_n,b_{\alpha })$$ . $$\hat {\mathscr {H}}$$ can be solved by any finite-size algorithms, such as exact diagonalization, QMC , DMRG [9, 23, 24], or finite-size PEPS [25–27] algorithms. The error from the rank-1 decomposition will be reduced since the loops inside the cluster will be fully considered.
Similar to the 1D cases, the ground-state properties can be extracted by the reduced density matrix $$\hat {\rho }(S)$$ after tracing over the entanglement-bath degrees of freedom. We have $$\hat {\rho }(S) = \text{Tr}_{/(S)} |\varPhi \rangle \langle \varPhi |$$ (with |Φ〉 the ground state of the infinite model) that well approximate by
$$\displaystyle \begin{aligned} \begin{array}{rcl} \hat{\rho}(S) \simeq \sum_{SS'b_1b_2 \cdots} \varPsi^{\ast}_{S'b_1b_2 \cdots} \varPsi_{Sb_1b_2 \cdots} |S \rangle \langle S'|, {} \end{array} \end{aligned} $$
(6.31)
with $$\varPsi _{Sb_1b_2 \cdots } $$ the coefficients of the ground state of $$\hat {H}_{FB}$$ .

Figure 6.6 illustrates the ground state ansatz behind the few-body model. The cluster in the center is entangled with the surrounding infinite tree brunches through the entanglement-bath degrees of freedom. Note that solving Eq. (6.20) in Stage one is equivalent to solving Eq. (6.28) by choose the cluster as one supercell.

Some benchmark results of simulating 2D and 3D spin models can be found in Ref. [5]. For the ground state of Heisenberg model on honeycomb lattice, results of the magnetization and bond energy show that the few-body model of 18 physical and 12 bath sites suffers only a small finite-effect of O(10−3). For the ground state of 3D Heisenberg model on cubic lattice, the discrepancy of the energy per site is O(10−3) between the few-body model of 8 physical plus 24 bath sites and the model of 1000 sites by QMC . The quantum phase transition of the quantum Ising model on cubic lattice can also be accurately captured by such a few-body model, including determining the critical field and the critical exponent of the magnetization.

6.4 Quantum Entanglement Simulation by Tensor Network: Summary

Below, we summarize the QES approach for quantum many-body systems with few-body models [1, 5, 6]. The QES contains three stages (Fig. 6.7) in general. The first stage is to optimize the physical-bath interactions by classical computations. The algorithm can be iDMRG in one dimension or the zero-loop schemes in higher dimensions. The second stage is to construct the few-body model by embedding a finite-size cluster in the entanglement bath, and simulate the ground state of this few-body model. One can employ any well-established finite-size algorithms by classical computations, or build the quantum simulators according to the few-body Hamiltonian. The third stage is to extract physical information by tracing over all bath degrees of freedom. The QES approach has been generalized to finite-temperature simulations for one-, two-, and three-dimensional quantum lattice models [6].
../images/489509_1_En_6_Chapter/489509_1_En_6_Fig7_HTML.png
Fig. 6.7

The “ab initio optimization principle” to simulate quantum many-body systems

As to the classical computations, one will have a high flexibility to balance between the computational complexity and accuracy, according to the required precision and the computational resources at hand. On the one hand, thanks to the zero-loop approximation, one can avoid the conventional finite-size effects faced by the previous exact diagonalization, QMC , or DMRG algorithms with the standard finite-size models. In the QES , the size of the few-body model is finite, but the actual size is infinite as the size of the defective TN (see Sect. 5.​3.​3). The approximation is that the loops beyond the supercell are destroyed in the manner of the rank-1 approximation, so that the TN can be computed efficiently by classical computation. On the other hand, the error from the destruction of the loops can be reduced in the second stage by considering a cluster larger than the supercell. It is important that the second stage would introduce no improvement if no larger loops were contained in the enlarged cluster. From this point of view, we have no “finite-size” but “finite-loop” effects. In addition, this “loop” scheme explains why we can flexibly change the size of the cluster in stage two: which is just to restore the rank-1 tensors inside the chosen cluster with the full tensors.

The relations among other algorithms are illustrated in Fig. 6.8 by taking certain limits of the computational parameters. The simplest situation is to take the dimension of the bath sites dim(b) = 1, and then $$\hat {\mathscr {H}}_{\partial }$$ can be written as a linear combination of spin operators (and identity). Thus in this case, v [x] simply plays the role of a classical mean field. If one only uses the bath calculation of the first stage to obtain the ground-state properties, the algorithm will be reduced to the zero-loop schemes such as tree DMRG and simple update of iPEPS . By choosing a large cluster and dim(b) = 1, the DMRG simulation in stage two becomes equivalent to the standard DMRG for solving the cluster in a mean field. By taking proper supercell, cluster, algorithms, and other computational parameters, the QES approach can outperform others.
../images/489509_1_En_6_Chapter/489509_1_En_6_Fig8_HTML.png
Fig. 6.8

Relations to the algorithms (PEPS, DMRG, and ED) for the ground-state simulations of 2D and 3D Hamiltonian. The corresponding computational set-ups in the first (bath calculation) and second (solving the few-body Hamiltonian) stages are given above and under the arrows, respectively. Reused from [5] with permission

The QES approach with classical computations can be categorized as a cluster update scheme (see Sect. 4.​3) in the sense of classical computations. Compared with the “traditional” cluster update schemes [26, 28–30], there exist some essential differences. The “traditional” cluster update schemes use the super-orthogonal spectra to approximate the environment of the iPEPS. The central idea of QES is different, which is to give an effective finite-size Hamiltonian; the environment is mimicked by the physical-bath Hamiltonians instead of some spectra.

In addition, it is possible to use full update in the first stage to optimize the interactions related to the entanglement bath. For example, one may use TRD (iDMRG , iTEBD , or CTMRG ) to compute the environment tensors, instead of the zero-loop schemes. This idea has not been realized yet, but it can be foreseen that the interactions among the bath sites will appear in $$\hat {H}_{FB}$$ . Surely the computation will become much more expensive. It is not clear yet how the performance would be.

The idea of “bath” has been utilized in many approaches and gained tremendous successes. The general idea is to mimic the target model of high complexity by a simpler model embedded in a bath. The physics of the target model can be extracted by integrating over the bath degrees of freedom. The approximations are reflected by the underlying effective model. Table 6.1 shows the effective models of two recognized methods (DFT and dynamic mean-field theory (DMFT) [31]) and the QES . An essential difference is that the effective models of the former two methods are of single-particle or mean-field approximations, and the effective model of the QES is strongly correlated.
Table 6.1

The effective models under several bath-related methods: density functional theory (DFT , also known as the ab initio calculations), dynamical mean-field theory (DMFT) , and QES

Methods

DFT

DMFT

QES

Effective models

Tight binding model

Single impurity model

Interacting few-body model

The QES allow for quantum simulations of infinite-size many-body systems by realizing the few-body models on the quantum platforms. There are several unique advantages. The first one concerns the size. One of the main challenges to build a quantum simulator is to access a large size. In this scheme, a few-body model of only O(10) sites already shows a high accuracy with the error ∼O(10−3) [1, 5]. Such sizes are accessible by the current platforms. Secondly, the interactions in the few-body model are simple. The bulk just contains the interactions of the original physical model. The physical-bath interactions are only two-body and nearest neighbor. But there exist several challenges. Firstly, the physical-bath interaction for simulating, e.g., spin-1∕2 models, is between a spin-1∕2 and a higher spin. This may require the realization of the interactions between SU(N) spins, which is difficult but possible with current experimental techniques [32–35]. The second challenge concerns the non-standard form in the physical-bath interaction, such as the $$\hat {S}^x\hat {S}^z$$ coupling in $$\hat {H}_{FB}$$ for simulating quantum Ising chain [see Eq. (6.15)] [18]. With the experimental realization of the few-body models, the numerical simulations of many-body systems will not only be useful to study natural materials. It would become possible to firstly study the many-body phenomena by numerics, and then realize, control, and even utilize these many-body phenomena in the bulk of small quantum devices.

The QES Hamiltonian was shown to also mimics the thermodynamics [6]. The finite-temperature information is extracted from the reduced density matrix
$$\displaystyle \begin{aligned} \hat{\rho}_R = \text{Tr}_{\text{bath}} \hat{\rho}, {} \end{aligned} $$
(6.32)
with $$\hat {\rho }=e^{- \hat {\mathscr {H}}_{\mathscr {F}\mathscr {B}}/\mathrm {T}}$$ the density matrix of the QES at the temperature T and Trbath the trace over the degrees of freedom of the bath sites. $$\hat {\rho }_R$$ mimics the reduced density matrix of infinite-size system that traces over everything except the bulk. This idea has been used to simulate the quantum models in one, two, and three dimensions. The QES shows good accuracy at all temperatures, where relatively large error appears near the critical/crossover temperature.

One can readily check the consistency with the ground-state QES . When the ground state is unique, the density matrix is defined as $$\hat {\rho } = |\varPsi \rangle \langle \varPsi |$$ with |Ψ〉 the ground state of the QES. In this case, Eqs. (6.32) and (6.16) are equivalent. With degenerate ground states, the equivalence should still hold when the spontaneous symmetry breaking occurs. With the symmetry preserved, it is an open question how the ground-state degeneracy affects the QES, where at zero temperature we have $$\hat {\rho } = \sum _a^{\mathscr {D}} |\varPsi _a \rangle \langle \varPsi _a|/ \mathscr {D}$$ with {|Ψ a〉} the degenerate ground states and $$\mathscr {D}$$ the degeneracy.

Creative Commons

Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.