Time evolution of quantum tensor

Pełen tekst

(1)TIME EVOLUTION OF QUANTUM TENSOR NETWORKS. Piotr Czarnik. A thesis submitted for the degree of Doctor of Philosophy Supervisor: prof. dr hab. Jacek Dziarmaga. Uniwersytet Jagielloński Instytut Fizyki im. Mariana Somluchowskiego. Kraków, May 2016.

(2) 2.

(3) Wydział Fizyki, Astronomii i Informatyki Stosowanej Uniwersytet Jagielloński. Oświadczenie Ja niżej podpisana/podpisany Piotr Czarnik doktorant Wydziału Fizyki, Astronomii i Informatyki Stosowanej Uniwersytetu Jagiellońskiego oświadczam, że przedłożona przeze mnie rozprawa doktorska pt. ,,Time Evolution of Quantum Tensor Networks’’ jest oryginalna i przedstawia wyniki badań wykonanych przeze mnie osobiście, pod kierunkiem prof. Jacka Dziarmagi. Pracę napisałem samodzielnie. Oświadczam, że moja rozprawa doktorska została opracowana zgodnie z Ustawą o prawie autorskim i prawach pokrewnych z dnia 4 lutego 1994 r. (Dziennik Ustaw 1994 nr 24 poz. 83 wraz z późniejszymi zmianami). Jestem świadomy, że niezgodność niniejszego oświadczenia z prawdą ujawniona w dowolnym czasie, niezależnie od skutków prawnych wynikających z ww. ustawy, może spowodować unieważnienie stopnia nabytego na podstawie tej rozprawy.. Kraków, dnia 01.06.2016 r.. ………………………………… podpis doktorantki/doktoranta.

(4) 4.

(5) Contents 0.1 0.2. Abstract. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Abstrakt. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 1 Quantum tensor networks 1.1 Area law for entanglement (ground states). . . . . . . . . . . 1.2 Graphical tensor notation. . . . . . . . . . . . . . . . . . . . . 1.3 Quantum tensor networks for ground state simulation. . . . . 1.4 Ground state simulations with MPS and PEPS. . . . . . . . . 1.5 Area law for thermal states. . . . . . . . . . . . . . . . . . . . 1.6 Quantum tensor networks for mixed states simulations. . . . 1.7 Thermal states of infinite systems. . . . . . . . . . . . . . . . 1.8 Quasi-exact tensor network representation of Gibbs operator.. . . . . . . . .. . . . . . . . .. . . . . . . . .. . . . . . . . .. . . . . . . . .. . . . . . . . .. . . . . . . . .. . . . . . . . .. . . . . . . . .. . . . . . . . .. . . . . . . . .. . . . . . . . .. . . . . . . . .. . . . . . . . .. 6 6 7 7 8 10 12 15 15 17 20. 2 Gibbs state simulation methods. 24 2.1 Imaginary time evolution of purification U (β). . . . . . . . . . . . . . . . . . . . . . . 24 2.2 Gibbs operator variational renormalization. . . . . . . . . . . . . . . . . . . . . . . . . 27 2.3 Comparison of the methods. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3 Benchmark results - summary. 3.1 General remarks. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 2D quantum Ising model. . . . . . . . . . . . . . . . . . . . . . . . 3.3 Isotropic 2D quantum compass model - frustrations of interactions. 3.4 Fermionic models. . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 The Hubbard model. . . . . . . . . . . . . . . . . . . . . . . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. . . . . .. 32 32 34 35 37 39. 4 Conclusions.. 41. 5 References. 42. 6 The articles.. 43. 5.

(6) 0.1. Abstract.. This PhD thesis consists of 5 articles. They introduce and benchmark novel methods of unbiased 2D Gibbs state simulation based on quantum tensor networks. The methods are suited for strongly correlated lattice models for which Quantum Monte Carlo suffers the sign problem. In the first chapter I present a short introduction to quantum tensor networks. The second chapter is a nontechnical introduction to the methods. It compares also the introduced methods. In the third chapter I summarize all benchmarks results presented in the articles and present additional benchmark results for the Hubbard model. The fourth chapter contains conclusions. The articles are included as the last chapter.. 0.2. Abstrakt.. Ta rozprawa doktorska składa się z 5 artykułów. Wprowadzają one i testują nowe metody symulacji stanów termicznych układów dwuwymiarowych. Metody te bazują na sieciach tensorowych. Ich celem jest symulacja silnie skorelowanych układów kwantowych dla których kwantowe Monte Carlo wykazuje problem znaku. Pierwszy rozdział rozprawy jest wprowadzeniem do sieci tensorowych. Drugi zawiera wprowadzenie do zaproponowanych metod. Trzeci rozdział podsumowuje wyniki uzyskane za pomocą tych metod. Czwarty rozdział zawiera wnioski. Artykuły zostały dołączone jako ostatni rozdział.. 6.

(7) 1 1.1. Quantum tensor networks Area law for entanglement (ground states).. Simulation of many-body quantum systems is a difficult problem because in their case Hilbert space H grows exponentially with number of degrees of freedom N dim H ∼ exp(N ).. (1). Many methods of quantum state simulation were proposed to deal with the problem. Nevertheless simulation of important classes of strongly correlated systems like high-TC superconductors or spin liquids still remains a challenge. Recently new approach to the problem - quantum tensor networks methods is intensively developed. I will quantum tensor networks briefly starting from methods suited for ground state simulation which are one of the best methods available for strongly correlated systems. I start from insight coming from quantum information theory giving theoretical motivation for them. Here I consider the case of lattice models. For lattice models the Hilbert space of a system H is constructed from Hilbert spaces of localized degrees of freedom. Here Hi consists of degrees of freedom localized on site i. The Hilbert space H is a tensor product: H = H1 ⊗ H2 ⊗ · · · ⊗ HN ,. (2). where N is a number of the lattice sites. To introduce the area law I need Von Neumann’s entanglement entropy S(ρA ) which quantifies entanglement between two parts of a system A, B such that H = HA ⊗ H B .. (3). I call A and B for which (3) is satisfied a system bipartition. An example of bipartition is shown in Fig. 1 I consider A which is a connected block of lattice sites. Than B is the rest of lattice. For a bipartition A, B and a state |Ψi: S(ρA ) = −Tr ρA log2 ρA ,. ρA = TrB |ΨihΨ|.. ρB = TrA |ΨihΨ|.. S(ρA ) = S(ρB ),. (4) (5). The entanglement entropy is nonnegative. It is minimal for a non-entangled product state: |Ψi = |ΨA i ⊗ |ΨB i =⇒ S(ρA ) = 0,. (6). and is maximized by a maximally entangled state: |Ψi =. d X 1 √ |iA i ⊗ |iB i, d i=1. d = min(dim HA , dim HB ) =⇒ S(ρA ) = log2 d,. (7). where {|iA i, i = 1, . . . , d} and {|iB i, i = 1, . . . , d} are sets of orthonormal states of A and B. The area law for entanglement is a property of some quantum states for which entanglement of a block of lattice sites A grows at most like its boundary surface area |∂A|: S(ρA ) ≤ O(|∂A|).. (8). Intuitive understanding of the area law can be obtained by considering a system with "entanglement correlation length" in which entanglement between degrees of freedom decreases rapidly for distances larger than the "entanglement correlation length". Then for blocks A significantly larger than the "entanglement correlation length", it can be expected that entanglement between degrees of freedom in the block A and degrees outside scales like the bloc boundary surface area |∂A|. There is a conjecture claiming that all ground states of so called local gapped Hamiltonians obey the area law [3]. Local Hamiltonians are defined as Hamiltonians with finite range of interactions and hoppings (I consider here only lattice model). Examples of such Hamiltonians are the Hubbard model or the Heisenberg model. A gapped Hamiltonians is defined as a Hamiltonian with an energy gap between the ground state and the first excited state. Proof of the conjecture was found in the cases of 1D Hamiltonians and quadratic bosonic and fermionic Hamiltonians [3]. It was also shown 7.

(8) B A. Figure 1: A bipartition of a square 5x5 lattice into the bloc of sites A and the rest of the lattice B. that ground states of quadratic gapless fermionic Hamiltonians obey the area law with logarithmic multiplicative correction [3]. For these Hamiltonians Von Neumann’s entanglement entropy S(L) of a D-dimensional hypercube with side length L obeys inequality: c1 LD−1 lnL ≤ S(L) ≤ c2 LD−1 ln2 L,. (9). where c1 , c2 are constants. Most of quantum states obey volume law for entanglement, which scales in their case for a bloc of sites as its volume. Furthermore in limit of an infinite system states with the area law become a subspace of measure zero in the Hilbert space [8]. These results motivated search for variational unbiased numerically tractable ansatzes for states with the area law. . Such ansatzes were found. They are based on parameterizations of states with the area law by contractions of tensors [4], [5]. These parameterizations are called quantum tensor networks or simply tensor networks.. 1.2. Graphical tensor notation.. I denote tensors and their contractions using graphical notation representing tensors by geometrical figures e.g spheres. Tensor indices are represented by lines coming out the geometrical figures. A contracted pair of tensor indices is represented by a line joining two geometrical figures (for a contraction of two tensors) or a loop attached to a geometrical figure (for a partial trace). The graphical notation usage is shown in Fig. 2.. 8.

(9) (a). (b). (c). (d). (e). (f). (g). Figure 2: The graphical tensor notation. In (a) a scalar s. In (b) a vector P va0. In00(c) a matrix Mab . In (d) a rank 5 tensor T . In (e) a matrix multiplication M = abcde ac P b Mab Mbc . In (f) a tensor P 0 0 0 contraction Tabcdef gh = i Tabcdi Tef . In (g) a partial trace T abc = ghi d Tabcdd .. 9.

(10) 1.3. Quantum tensor networks for ground state simulation.. For a quantum lattice model with N lattice sites probability amplitudes of a state |ψi are given by a rank N tensor ψi1 ,i2 ,...,iN . X H = H1 ⊗ H2 ⊗ · · · ⊗ HN ⇒ |ψi = ψi1 ,i2 ,...,iN |i1 i2 . . . iN i. (10) i1 ,i2 ,...,iN. Here I use a tensor product basis of the Hilbert space H: X |i1 i2 . . . iN i = |i1 i ⊗ |i2 i ⊗ · · · ⊗ |iN i,. (11). i1 i2 ...iN. where |i1 i, |i2 i, . . . , |iN i are basis vectors of H1 ,H2 ,. . . ,HN respectively. A representation of ψi1 i2 ...iN tensor by a contraction of smaller tensors is an example of a quantum tensor network. States obeying the area law can be constructed as tensor networks in following way. With each lattice site s a tensor As is associated. One of As indices is (called a physical index) becomes one of ψi1 ,i2 ,...,iN indices after all As tensors are contracted. The index is enumerates degrees of freedom at site s. All other indices of As are contracted to obtain ψi1 ,i2 ,...,iN (they are called virtual indices). Each of As virtual indices is associated with a different lattice bond connecting s with one of its nearest neighbor sites. Therefore rank of the tensor As - rs is determined by a number of nearest neighbors of site s - nns : rs = nns + 1. (12) To obtain ψi1 i2 ...iN pairs of virtual indices associated with the same lattice bond are contracted. Therefore a diagram denoting ψi1 i2 ...iN mirrors the lattice structure. For 1D lattice a tensor network obtained in such way is called a matrix product state. It is also sometimes called a MPS. A MPS is a contraction of rank 3 tensors. In the case of open boundary conditions rank 2 tensors appear at the first and the last site (see Fig. 3). For 2D square lattice pair entangled projected states (called also PEPS) are obtained. They are contractions of rank 5 tensors. Again in the case of open boundary conditions tensors with smaller ranks (rank 4 or rank 3 tensors) appear at the boundary (see Fig. 4). Other two dimensional lattices frequently can be mapped into a square lattice to enable their simulation with pair entangled pair operators. Another quantum tensor networks representing states obeying the area law can be proposed but usually numerical computations done with them are more costly than analogous computations done with MPS or PEPS tensor networks. Therefore MPS and PEPS tensor networks are the most popular ones. I will call As tensors MPS tenors or PEPS tensors respectively. Usually dimensions of all virtual indices are chosen to be the same. Then the virtual indices dimension is called bond dimension D. I will call both the introduced tensor network representations of states and state represented by them matrix product states and pair entangled projected states respectively. It can be easily shown that the area law is satisfied for MPS and PEPS quantum tensor networks. Here I present a proof for PEPS. I consider a system bipartition into degrees of freedom in a bloc of sites A and degrees of freedom in the rest of the lattice B. State |Ψi of the system is represented by a PEPS. I introduce two tensors created by a contraction of the PEPS tensors in the block - ψiAA ,k and. i1. i2. i3 i4 ψi1,i2,i3,i4,i5,i6. i5. i6. Figure 3: A Matrix Product State (a MPS) representing a state of a 1D lattice model with open boundary conditions. The lattice has 6 sites. The diagram represents tensor of probability amplitudes ψi1 i2 i3 i4 i5 i6 . The MPS tensors are represented by the circles. The lines connecting the circles represent the virtual indices. The remaining lines represent the physical indices. 10.

(11) Figure 4: A Pair Entangled Projected State (a PEPS) representing a state of a 2D lattice model with open boundary conditions. The lattice is a a square 4x4 lattice. the diagram represents a tensor of probability amplitudes ψP EP S . The PEPS tensors are represented by the spheres. As in Fig 3 the lines connecting the spheres represent the virtual indices and the remaining lines represent the physical indices. the PEPS tensors in the rest of the lattice ψiBB ,k (see Fig. 5). Here iA is an index which numerates all degrees of freedom in A, iB is an index which numerates all degrees of freedom in B, and k is an index which replaces all uncontracted virtual indices (its dimension is a product of dimensions of all replaced indices). For each value of k index ψiAA ,k and ψiBB ,k are probability amplitudes of unnonormalized states of A and B: |kiA and |kiB . Now the state |Ψi and the reduced density matrices ρA , ρB are: X X X |Ψi = |kiA ⊗ |kiB ρA = ρk,k0 |kiA A hk 0 | ρB = ρk,k0 |kiB B hk 0 |. (13) k. k,k0. k,k0. Rank of the reduced density matrices ρA , ρB - r is now bounded from above by dimension of k. Number of indices replaced by k is proportional to the length of A boundary |∂A|. In the case of the PEPS with the bond dimension D: r ≤ D|∂A| , so the entanglement entropy S(ρA ) is bounded from above: S(ρA ) ≤ log2 r = |∂A|log2 D.. (14). We obtain that all PEPS with the bond dimension D obey the area law by construction. Furthermore the upper limit of block entanglement (14) grows with the growing bond dimension D. This result suggests that in the case of strongly entangled state approximated by PEPS with growing D the better state approximation can be obtained. In the limit of D = 1 a PEPS state becomes a product state: |Ψi = |Ψ1 i ⊗ |Ψ2 i ⊗ · · · ⊗ |ΨN i. It was checked empirically that the bond dimension D necessary to obtain a good state approximation grows with growing entanglement (see e.g [13],[10]). In the case of small entanglement sometimes even the smallest nontrivial D = 2 suffices to obtain qualitatively correct results (see e. g. [14]).. 11.

(12) Figure 5: An example of tensors ψkiA and ψkiB appearing in the proof of the area law for PEPS. Here they correspond to the bipartition of the lattice into a 2x2 bloc of sites A at the bottom left and the rest of the lattice B. The index k replaces the uncontracted virtual indices k1 , k2 , k3 , k4 .. 1.4. Ground state simulations with MPS and PEPS.. It was shown that leading methods of ground state simulation in the case of 1D strongly correlated system (Density Matrix Renormalization Group - DMRG) are in fact methods of finding the optimal MPS ground state representation. This discovery led to further development of DMRG methods [1]. A variational ansatz for ground states of 2D infinite systems [9] based on PEPS tensor networks were proposed. The ansatz is one of the leading methods of unbiased ground state simulation for models with the notorious sign problem. For the t-J model and the Hubbard model (see Figs 6, 7) the PEPS ansatz gives the lowest variational energies beating the best variational Quantum Monte Carlo methods [11], [13]. Another class of models with the sign problem are spin models with frustrations of interactions. An example of the PEPS ansatz power in their case is simulation of Shastry-Sutherland model which led to solution of problem of SrCu2 (BO3 )2 magnetization plateaus (see Fig. 8).. 12.

(13) Figure 6: Variational energies of ground state candidates obtained by PEPS simulation for the t-J model plotted versus the inverse bond dimension 1/D. Favorable for high-Tc superconductivity values of parameters J/t = 0.4 and the hole density δ = 0.12 were chosen. The PEPS results for an infinite system were compared with the best variational Monte Carlo energies obtained by Fixed Node Monte Carlo (FNMC - the horizontal lines). FNMC gives results for a finite cluster of N lattice sites. The PEPS variational energies beat the ones obtained by FNMC. The FNMC energies additionally grow with growing N . PEPS results for different ground state candidates are shown. The figure is taken from P. Corboz, T. M. Rice, M. Troyer, Phys. Rev. Lett. 113, 046402 (2014).. Figure 7: Variational energies of a ground state candidate obtained by PEPS simulation for the Hubbard model (again favorable for high-Tc superconductivity values of parameters U/t = 0.875 and the electron density n = 0.875 were chosen). The energies are plotted versus the inverse bond dimension 1/D (the squares) and another more precise results convergence measure w (the dots). PEPS results (for an infinite system) again were compared with the best Monte Carlo variational energies (obtained by Fixed Node Monte Carlo) beating the Monte Carlo results. The results were also compared with other non-variational methods whose systematic errors are not controlled, so their results can’t be considered as reliable reference points. The figure is taken from P. Corboz, Phys. Rev. B 93, 045116 (2016). 13.

(14) Figure 8: SrCu2 (BO3 )2 magnetization in external magnetic field. This compound is a real world realization of a spin Shustry-Sutherland model with frustrations of interactions. Comparison of results obtained by PEPS ground state simulations with experimental results and the other numerical methods shows power of PEPS simulations. PEPS results obtained for the model led to identification of spin orders realized in magnetization plateaus and to an analytical model explaining the realized spin orders. The problem of nature of the plateaus had remained unsolved by 30 years. The figure is taken from Y. H. Matsuda, N. Abe, S. Takeyama, H. Kageyama, P. Corboz, A. Honecker, S. R. Manmana, G. R. Foltin, K. P. Schmidt, and F. Mila, Phys. Rev. Lett. 111, 137204 (2013).. 14.

(15) 1.5. Area law for thermal states.. Density matrix of a thermal state (Gibbs state) ρ(β) equals: ρ(β) =. 1 −βH e . Z. (15). Here β is the inverse temperature T : β=. 1 , T. (16). and Z is the partition function: Z = Tr e−βH .. (17). Quantum thermal states of local Hamiltonians obey area law for mutual information [6]. For mixed states measure of quantum correlations between two system parts A, B is mutual information I(A : B). For a mixed state with the density matrix ρ: I(A : B) = S(ρA ) + S(ρB ) − S(ρ),. (18). where ρA , ρB are reduced density matrices of A, B and S(ρ) = −Tr ρ log2 ρ.. (19). It was shown that Gibbs states of local Hamiltonians obey area law for the mutual information analogous to the area law for pure states [6]: I(A : B) ≤ cβ|∂A|,. (20). where |∂A| is the block boundary surface area, and c is a constant dependent on the Hamiltonian. Although mutual information does not vanish in the case of classical thermal states (for which it also obeys the area law I(A : B) ≤ log2 d |∂A|, (21) where d is a number of degrees of freedom per lattice site), quantum systems in general have much larger mutual information than classical ones. Mixed states obeying the area law again can be parametrized by tensor networks.. 1.6. Quantum tensor networks for mixed states simulations.. For lattice models density matrix elements create a rank 2N tensor ρi1 ,i2 ,...,iN ,j1 ,j2 ,...jN : X ρ= ρi1 ,i2 ,...,iN ,j1 ,j2 ,...,jN |i1 , i2 , . . . , iN ihj1 , j2 , . . . jN |.. (22). i1 ,i2 ,...,iN ,j1 ,j2 ,...,jN. MPS and PEPS tensor networks can be easily generalized to obtain tensor networks representing ρi1 ,i2 ,...,iN ,j1 ,j2 ,...,jN . An additional physical index js has to be added to MPS tensors (or PEPS tensors) As . The new index becomes one of ρi1 ,i2 ,...,iN ,j1 ,j2 ,...,jN indices after the contraction. The generalization gives me tensor networks called Matrix Product Operators (MPO) and Pair Entangled Projected Operators (PEPO) (see Figs 9, 10). By analogy with PEPS and MPS tensors networks I use PEPO and MPO to call both the tensor networks and the states represented by them. I call tensors As in their case PEPO tensors or MPO tensors. For these new tensor networks dimension of their virtual indices is called the bond dimension D as in the case of PEPS and MPS tensor networks.. 15.

(16) j1. j2. j3. j4. j5. j6. i1. i2 i3 i4 i5 ρi1,i2,i3,i4,i5,i6, j1, j2, j3, j4, j5, j6. i6. Figure 9: A Matrix Product Operator (a MPO) representing a mixed state of a 1D lattice model with open boundary conditions. The lattice has 6 sites. The diagram represents the density matrix tensor ρi1 ,i2 ,i3 ,i4 ,i5 ,i6 ,j1 ,j2 ,j3 ,j4 ,j5 ,j6 . The MPO tensors are represented by the circles.. Figure 10: A Pair Entangled Projected Operator (a PEPO) representing a mixed state of a 2D lattice model with open boundary conditions. The lattice is a square 4x4 lattice. The diagram again represents the density matrix tensor ρP EP O . The PEPO tensors are represented by the spheres.. 16.

(17) Pair Entangled Projected Operators obey the area law by construction [6] I(A : B) ≤ 2log2 D |∂A|.. (23). Comparison of (23) with (21) shows that even with very small D they can represent Gibbs states with non-classical correlations. It was shown [7] that Gibbs states of local quantum spin lattice models can be approximated with any precision by PEPO with number of parameters scaling only polynomially with number of lattice sites N . In the case of a generic local quantum lattice model earlier the sub-exponential scaling had been obtained [8]. These results proof that pair entangled projected operators are indeed efficient parametrization of Gibbs states of local Hamiltonians encouraging their usage as variational numerical ansatzes for such thermal states. Search for numerical ansatzes based on PEPO tensor networks is subject of this thesis.. 1.7. Thermal states of infinite systems.. One can define a PEPO tensor network for an infinite lattice in the same way as for a finite lattice. In general case with each lattice site a different tensor may be associated leading to the infinite number of parameters. But in the case of a PEPO with translational symmetry the number of parameters is finite. Such PEPO can be used to simulate state of an infinite system with the translational symmetry numerically. In the simplest case when the simulated state is translationally invariant one can use a translationally invariant PEPO which is obtained by associating with each lattice site a copy of the same PEPO tensor A. In the case of an infinite system I will search for a tensor network representation of the Gibbs operator e−βH instead of tensor network representation of the density matrix ρ(β) to avoid problem of finding the normalization. Furthermore I will represent the Gibbs operator using its purification U (β): e−βH = U (β)U (β)† . (24) U (β) becomes a pure state when its right index is treated as an index numerating states of an ancillary Hilbert space (which is just a copy of the system Hilbert space), therefore it is called the thermal state purification. Representation (24) is used to guarantee hermitianity of numerically found e−βH . PEPO representation of U (β) is shown in Fig. 11(a). U (β)† PEPO representation can be easily found. Elements of the tensor U (β)† are obtained from elements of the U (β) tensor: ∗ U (β)†...,ik−1 ,ik ,ik+1 ,...,jk−1 ,jk ,jk+1 ,... = U (β)...,jk−1 ,jk ,jk+1 ,...,ik−1 ,ik ,ik+1 ,... , (25) where ∗ denotes complex conjugation. To satisfy (25) the new PEPO tensor A˜ should be chosen as follows: ∗ A˜i,j,t,l,b,r = Aj,i,t,l,b,r . (26) Here t, l, b, r are the virtual indices and i, j are the physical indices. The new PEPO is shown in Fig. 11(c). A tensor network representation of the Gibbs operator U (β)U (β)† is shown finally in Fig. 11(e). A PEPO representation of the Gibbs operator e−βH can be obtained from the U (β) and U (β)† PEPO representations. A contraction of PEPO tensors A and A˜ associated with the same lattice site (see Fig. 12(b)) gives the new PEPO tensor B. Both tensor network representations of the Gibbs operator are shown in Fig. 12(a). The Gibbs operator e−βH is an infinite rank tensor. But the partition function Tr e−βH is a scalar. It can be written as a contraction of rank 8 tensors a obtained by taking a partial trace of the B PEPO tensor (the physical indices are traced out) - see Figs 13(a),(d). Another scalar is Tr e−βH O, where O is an observable. In the simplest case when O acts on 1 site I can represent Tr e−βH O by a contraction of a tensors and one o tensors obtained as a contraction of B tensor with the observable O - see Figs 13(b),(e). The tensor network representations of Tr e−βH and Tr e−βH O differ only at one site. I call a tensor network obtained as a contraction of all a tensor expect one associated with the site at which O acts an 1-site environment Ea (see Fig. 13(c)). Now Tr e−βH and Tr e−βH O can be written as: Tr e−βH = Tr Ea a, (27) Tr e−βH O = Tr Ea o, 17. (28).

(18) (a). (b). (c). (d). (e). Figure 11: Tensor network representation of a translationally invariant Gibbs state. The system is infinite. In (a) the translationally invariant PEPO representing U (β). The PEPO is a contraction of copies of the same PEPO tensor Ai,j,t,l,b,r (shown in (b)). The copies of Ai,j,t,l,b,r are represented by the green spheres. The dashed lines denote that only a part of the infinite diagram is shown. In (c) the PEPO representing U (β)† . This PEPO is the contraction of copies of the PEPO tensor A˜i,j,t,l,b,r (shown in (d)). A˜ elements are defined in (26). In (e) a tensor network representation (24) of the Gibbs operator e−βH . Pairs of A, A˜ physical indices contracted to obtain e−βH are denoted by the red lines. Physical indices of A, A˜ which become indices of e−βH are denoted by the blue lines. 18.

(19) (a). (b). Figure 12: In (b) the new Bi,j,td ,ld ,bd ,rd is obtained from a rank 10 tenP PEPO tensor ˜ 0 0 0 0 sor Ci,j,t,t0 ,l,l0 ,b,b0 ,r,r0 = A A by replacing pairs of the virtual indices k i,k,t,l,b,r k,j,t ,l ,b ,r (t, t0 ), (b, b0 ), (r, r0 ), (l, l0 ) associated with the same lattice bond by new virtual indices td , ld , rd , bd . The bond dimension of the new PEPO equals D2 , where D is the old PEPO bond dimension. The new PEPO virtual indices are denoted by the thick lines. In (a) the new PEPO representation of e−βH and its old tensor network representation are shown.. 19.

(20) (a). (b). (c). (d). (e). Figure 13: In (d) (the left panel) the partition function Tr e−βH is represented with a tensor obtained as a partial trace of the B PEPO tensor over the physical indices as shown in (a). In (e) (the left panel) the tensor network representation of Tr e−βH O is shown using o tensor introduced in (b). Here o is an 1-site observable. In (c) a 1-site environment Ea is defined as a contraction of all copies of a appearing in the Tr e−βH representation expect a copy corresponding to a site on which O acts. In (d) and (e) (the right panels) Tr e−βH and Tr e−βH O are represented using the environment Ea . where Tr Ea a and Tr Ea o just mean contractions of a, o tensors with the environments giving the desired scalars (see Figs 13(d),(e)). So having Ea I can compute straightforwardly expected values of all 1-site observables: Tr e−βH O Tr Ea o hOi = = . (29) −βH Tr e Tr Ea a Actually to compute hOi it is enough to know Ea up to its normalization. It turns out that there exist numerical method which enables controlled approximation of the environment Ea up to its normalization. It is called CTMRG (Corner Transfer Matrix Renormalization Group) [22, 23]. Expected values of 2-site observables also can be obtained with CTMRG. Their computation require computation of two site environments. The two site environment is a contraction of a tensors appearing in the Tr e−βH tensor network representation expect of two of them associated with sites at which the two site observable acts.. 1.8. Quasi-exact tensor network representation of Gibbs operator.. A Gibbs operator e−βH can be approximated in a controlled way by a 3D tensor network. This tensor network is much more difficult to work numerically with then a PEPO tensor network because of larger computational complexity of tensor contractions. Nevertheless it is a starting point for one of the 20.

(21) proposed methods of Gibbs state simulation. The 3D tensor network is obtained by Suzuki-Trotter decomposition and one of elementary matrix decompositions (the singular value decomposition). I will introduce Suzuki-Trotter decomposition in the case of a Hamiltonian H which is a sum of 1 2 two non-commuting classical Hamiltonians Hcl , Hcl with the nearest neighbor interactions: 1 2 H = Hcl + Hcl ,. 1 2 [Hcl , Hcl ] 6= 0.. (30). A classical Hamiltonian Hcl is a sum of commuting terms Hk : X Hcl = Hk , ∀k,l [Hk , Hl ] = 0.. (31). k. For H Suzuki-Trotter decomposition gives an approximation of U (β) = e−(β/2)H by exponents of the 1 2 classical Hamiltonians Hcl and Hcl . The quality of the approximation is controlled by a small time step dβ. Using second order Suzuki-Trotter decomposition I obtain: U (β) = e−βH/2 = UST (dβ)UST (dβ) . . . UST (dβ) +O(dβ), {z } |. 1. 2. 1. UST (dβ) = e−dβHcl /4 e−dβHcl /2 e−dβHcl /4 .. β/dβ terms. (32) An exponent of a classical Hamiltonian usually can be exactly represented by a PEPO with the bond dimension D small enough to be numerically tractable. I consider here a classical translationally invariant Hamiltonian Hcl which is a sum of commuting nearest neighbor interactions Hi,j . Here i, j are indices of sites on which Hi,j acts. The translational invariance is obtained if each term Hi,j is a copy of the same rank 4 tensor Hi1 ,j1 ,i2 ,j2 . An exponent e−cHcl is decomposed as a product of Hij exponents: Y e−cHcl = e−cHij , (33) <ij> −cHij. as shown in Fig. 14(b). Here c is a constant. e ≡ Mi1 ,i2 ,j1 ,j2 is a rank 4 tensor. Replacing two indices i1 , i2 associated with site i by one index id and replacing indices j1 , j2 associated with site j by one index jd I obtain a matrix Mid ,jd . The singular value decomposition of Mid ,jd gives: Mid ,jd =. Dcl X l=1. Uid ,l Sl,l Vjd ,l =. Dcl X. Lid ,l Rjd ,l ,. Lid ,l = Uid ,l. p. Sl,l ,. Rjd ,l = Vjd ,l. p. Sl,l .. (34). l=1. Here U † U = 1Dcl , V † V = 1Dcl , Sl,l > 0 and Dcl ≤ d2 ,. (35). where d is a number of degrees of freedom per a lattice site. Going back from id to i1 , i2 indices and from jd to j1 , j2 indices I obtain a decomposition of Mi1 ,i2 ,j1 ,j2 X Mi1 ,i2 ,j1 ,j2 = Li1 ,i2 ,l Rj1 ,j2 ,l . (36) l. Further I consider only the most symmetric case of U = V leading to translational invariance of the searched PEPO representation. Then L = R ≡ P and X Mi1 ,i2 ,j1 ,j2 = Pi1 ,i2 ,l Pj1 ,j2 ,l (37) l. (see Fig. 14(a)). Using (37) I obtain a new tensor network representation of e−cHcl shown in 14(b). Replacing the contraction of all tensors P associated with the same lattice site I obtain the searched PEPO tensor Acl (see Fig 14(c)). Dcl becomes the new PEPO bond dimension. Bound (35) usually ensures that the obtained PEPO is numerically tractable. 1 2 Using the PEPO representations of the classical Hamiltonians Hcl , Hcl and (32) I obtain the 3D tensor network representation of U (β) which is a contraction of 3β/dβ PEPO layers (see Fig. 15).. 21.

(22) (a). (b). (c). Figure 14: In (a) the decomposition (37) of e−cHij ≡ Mi1 ,i2 ,j1 ,j2 (the green tube) into two P tensors (the green spheres). In (b) e−cHcl decomposition into Mi1 ,i2 ,j1 ,j2 tensors (see (33)) (the left panel) and its further decomposition into P tensors (the right panel). Here the sliced tubes and the dashed lines denote that only a part of the infinite diagram is shown. In (c) construction of the searched Acl PEPO tensor.. 22.

(23) Figure 15: The 3D tensor network approximate representation of U (β). Each PEPO layer represents one exponent of the classical Hamiltonian in (32). The dashed blue lines mean that only the part of the whole diagram is shown. The red lines should be contracted to obtain e−βH .. 23.

(24) 2 2.1. Gibbs state simulation methods. Imaginary time evolution of purification U (β).. In the first three articles PEPO representation of the purification U (β) was obtained by its imaginary time evolution. I will consider here this method for translationally and rotationally invariant Gibbs operator. Furthermore I will consider Hamiltonian (30) for which PEPO representations of exponents 1 2 of classical Hamiltonians Hclass and Hclass are translationally and rotationally invariant. PEPO representation of U (β+dβ) is obtained here step by step starting from PEPO representation of U (β): 1 ˜1 (β + dβ) = e−(dβ/4)Hcl U U (β), (38) 2 ˜2 (β + dβ) = e−(dβ/2)Hcl ˜1 (β + dβ), U U. U (β + dβ) = e. 1 −(dβ/4)Hcl. ˜2 (β + dβ). U. (39) (40). In the first step (38) is modified to: 1 ˜1 (dβ) = e−(dβ/4)Hcl U .. (41). In each step after the first one (41) a new tensor network is obtained (see Fig. 16(a)). It is a contraction of rank 10 tensors (see Fig. 16(b)) Anonren . I denote it by U nonren . The rank 10 tensor Anonren is a contraction of two PEPO tensors A1 and A2 with the bond dimensions D1 and D2 . It can be in principle replaced by a new PEPO tensor analogically to Fig. 12(b), but the new bond dimension would be D1 D2 so a renormalization scheme replacing Anonren by a PEPO tensor Aren with the bond dimension Dren < D1 D2 is necessary to prevent exponential growth of the bond dimension. The renormalized PEPO representation U ren is obtained from U nonren after at each lattice bond a projector P = XX † is inserted (as shown in 17(a)). Here X is an isometry preserving distance in a Dren -dimensional subspace so P is an orthogonal projection on a Dren -dimensional subspace. Contracting the isometries X with tensors Anonren I obtain the renormalized PEPO tensor Aren see Fig. 17(b). The projectors P are chosen to minimize discrepancy between norms of U (β + dβ) and U (β + dβ). Such variational problem has to be solved in approximate way as the discrepancy is a function of infinite number of the projector P copies. To solve the problem I initialize the projector P (usually by P from one of the earlier iterations). I calculate the projector environment EP defined as norm of U nonren with one copy of the projector P removed. ||U nonren || = Tr U nonren (U nonren )† = Tr P EP . (42) In the case of step (40) EP is the partition function Z with one copy of the projector P removed. The environment EP can be computed straightforwardly using two-site environment mentioned in sub-chapter 1.7. For the case of translationally and rotationally invariant U nonren the same EP is obtained for all choices of the removed P copy. I diagonalize symmetric part of EP : EP + EPT =. DX 1 D2. Vi Λi Vi†. i < j ⇒ |Λi | > |Λj |.. (43). i=1. I choose the updated projector P as P =. D ren X. Vi Vi† .. (44). i=1. This choice of the updated P minimizes Tr P EP − Tr EP , which in the case of EP independent from P measures how much norm of U nonren (β + dβ) is distorted by insertion of the projector P . I compute projector environment EP for the updated projector P and again update the projector P . This procedure is repeated until convergence measured by convergence of observables is reached to guarantee self-consistency of the obtained solution. Simulation of the Gibbs state for a given β requires for the proposed method 3 ∗ β/dβ renormalizations (actually usually in few first steps one may straightforwardly replace Anonren (β + dβ) by a PEPO tensor without truncation, because in their case D1 D2 is smaller than the final desired bond dimension D, so the number of renormalizations is slightly smaller). 24.

(25) (a). (b). Figure 16: In (a) construction of the tensor network representation of U nonren obtained in one of steps of the imaginary time evolution (38), (39, (40). The representation is a contraction of copies of the rank 10 Anonren tensor shown in (b). Anonren is a contraction of two PEPO tensors. The orange 1 2 one comes from PEPO representation of one of the exponents of the classical Hamiltonians Hcl , Hcl . ˜ ˜ The green one comes from PEPO representation of one of U (β), U1 (β + dβ) or U2 (β + dβ) tensors. The red indices of U nonren should be contracted to obtain e−βH. 25.

(26) (a). (b). (c). Figure 17: In (a) copies of the projector P (shown in (c)) are inserted at all lattice bonds (in the left part) into the tensor network representation of U nonren to obtain its approximate PEPO representation U ren (shown in the right part). The projector P is a contraction of two isometries X represented by the orange cones as shown in (c). Actually P is a rank 4 tensor so strictly speaking it becomes projector when pairs of its left and right indices are replaced by one index. I call a rank 3 X tensors isometries analogically. In (b) the isometries X are used to truncate Anonren tensor giving the new PEPS tensor Aren with the truncated bond dimension Dren . In the figure indices with the bond dimension Dren are represented by the orange lines. It is easy to see that the renormalization preserves translational and rotational invariance of U nonren .. 26.

(27) The simplified version of the renormalization procedure was proposed at the beginning. The simplified version doesn’t guarantee self-consistency and has larger computational complexity so it is inferior to the presented above version. In the simplified version the updated projector P is computed in the first iteration from the projector environment EP as before but EP is computed for a trivial projector P = 1D1 D2 . Now EP computation requires working with a PEPO with large bond dimension D1 D2 . This PEPO comes from replacing Anonren by a PEPO tensor analogically to Fig. 12(b). In the earlier procedure the renormalized PEPO is used to compute EP . Furthermore in the simplified update there are no additional iterations of the P update which are necessary to ensure self-consistency. The imaginary time evolution of the purification with the simplified update was proposed in the first article and generalized to fermionic case in the second one. It was developed to the self-consistent version in the third article.. 2.2. Gibbs operator variational renormalization.. Another way to find PEPO representation of a Gibbs operator is its variational renormalization. I start from Gibbs operator representation introduced in sub-chapter 1.8. Here I consider the case of classical Hamiltonian (30) with nearest neighbor interactions. Furthermore I assume that the U (β) representation which is renormalized is translationally and rotationally invariant. The U (β) representation (see Fig. 15) consists of columns of PEPO tensors. One can obtain from it straightforwardly a PEPO representation by replacing all virtual indices associated with the same lattice bond by one new virtual index. After the replacement the column is turned into a PEPO tensor AST with the bond dimension DST = 23β/dβ . (45) . The obtained bond dimension DST is huge enough to forbid any numerical use of the obtained PEPO as in typical applications β/dβ > 100. A renormalization scheme has to be proposed which transforms the tensor network representation from Fig. 15 into a PEPO with the bond dimension small enough to be numerically tractable. For a sake of clarity I replace a contraction of three neighboring tensors A1cl , A2cl and A1cl associated with the same UST (dβ) in (32) by a PEPO tensor A0 with the bond dimension 8 as shown in Fig. 18(a). After this replacement U (β) is represented by Nβ = β/dβ layers of PEPO tensors T0 contracted as in Fig. 15. To renormalize the Gibbs operator representation from Fig. 15. I insert a projector Plarge at each set of virtual bonds associated with the smae lattice bond as shown in Fig. 18(d). The projector Plarge is an orthogonal projection on a Dren dimensional subspace given by an isometry Xlarge . † Plarge = Xlarge Xlarge. (46). The isometry Xlarge applied to the column of tensors T0 turns it into a PEPS tensor with the truncated bond dimension Dren as shown in Fig. 18(c). The isometry Xlarge is a huge tensor with DST ∗ Dren elements. So it can’t be represented numerically without additional constraints. The additional constraint is an assumption that the isometry Xlarge is a contraction of smaller isometries (see Fig.19(b)). The small isometries are rank 3 tensors truncating two neighboring (along the axis of imaginary time β - orthogonal to the PEPO planes in Fig. 15) indices into a one new index. The are arranged in layers. I have in total nβ = log2 Nβ layers. Final renormalized PEPO U (β) representation is obtained step by step. At k-th step the layer of isometries Xk acts on U (β) representation from the previous on (with Nβ /2k−1 PEPO layers and with PEPO tensors Ak−1 ) giving new U (β) representation with Nβ /2k PEPO layers and with PEPO tensors Ak (see Fig.19(a),(c)). The bond dimension of all Ak tensors equals Dren . The presented here Xlarge structure is called an example of a Tree Tensor Network (TTN). To preserve rotational and translational invariance of U (β) representation I use a copy of the same projector Plarge at each lattice bond. Elements of the isometries Xk are variational parameters. It can be shown that partition function Z representation constructed from T0 tensors is invariant under translations along imaginary time β axis Therefore I can choose all isometries Xk in the same layer as a copies of the same tensor. Now number of variational parameters (nβ − 1) ∗ D3 + D2 ∗ D0 scales logarithmically with β. All isometries Wk are chosen to minimally distort the partition function Z. The partition function Z is a function of an infinite number of copies of isometries Xk . Therefore the variational problem 27.

(28) (a). (b). (c). (d). Figure 18: Construction of the renormalized PEPO representation of U (4dβ) by variational Gibbs operator renormalization is presented. In (a) a contraction of three neighboring tensors A1cl , A2cl and A1cl associated with the same UST (dβ) in (32) is replaced by a PEPO tensor T0 . In (d) the renormalized PEPO representation is obtained starting from the unrenormalized one by insertion of projectors Plarge (46) at all "lattice bonds". In (c) columns of 4 T0 tensors (shown in (b)) are renormalized by contraction with isometries Xlarge represented by the orange cones giving the renormalized PEPO tensors T2 with the bond dimension Dren . The orange lines represent indices with dimension Dren .. 28.

(29) (a). X1 X2 X1. ≡ Xlarge. X3 X1 X2 X1 (b). (c). Figure 19: Construction of the renormalized PEPO representation of U (8dβ). The TTN structure of Xlarge is taken into account. Xlarge is a contraction of smaller isometries X1 , X2 , X3 as shown in (b). To obtain the final PEPO representation of U (8dβ) the isometries Xn are applied step by step. In each step Xn isometries renormalize a contraction of two Tm−1 PEPO tensors giving a new Tm PEPO tensor as shown in (a). In (c) initial U (8dβ) representation with T0 PEPO tensors and 3 renormalized representations obtained in the consecutive steps. Again the orange lines have dimension Dren . 29.

(30) is solved approximately. At the beginning isometries Xk are chosen randomly (another choices of isometries speeding up convergence of the variational update can also be proposed). For the initial choice of isometries Xk one can compute straightforwardly the PEPO tensor Tnβ . Using 1-site environment Ea one can compute straightforwardly environments of all isometries Xk . The environment of the isometry Xk - EXk is defined as the partition function Z with one copy of the isometry Xk removed: Z = Tr Xk EXk . (47) Here Tr means a contraction of the rank 3 tensor Xk with its environment EXk giving the partition function Z. Right now Xk is updated. The updated Xk is given by the Singular Value Decomposition (SVD) of EXk : EXk = U SV † ⇒ Xk = V U † . (48) This choice of Xk maximizes Z if EXk is independent from all isometries Xk . After all Xk are updated the new tensor Tnβ is obtained. The update procedure: Xk. →. Tnβ. →. EXk. →. Xk .. is repeated until the convergence of observables is reached. Why Xk was chosen to maximize the partition function Z? To answer this question one need to consider the environment of the projector P large - EP large . It can be shown that for the translationally and rotationally invariant representation of the Gibbs operator EP large is non-negatively defined. Therefore for EP large treated as independent from the isometries Xk each choice of isometries Xk decreases the partition function Z or doesn’t change it.. 2.3. Comparison of the methods.. The main advantage of Gibbs operator variational renormalization is that it optimize the whole projector using the partition function Tr e−βH as a figure of merit. The U (β) representation obtained by the imaginary time evolution is a function of many projectors. Each of them optimizes fidelity of U (β 0 ) representation for β 0 ≤ β. I am interested only in quality of e−βH approximation for given the bond dimension D. There is no reason that the imaginary time evolution representation is as good from this point of view as the representation directly targeting Tr e−βH . Indeed it turns out that the imaginary time evolution requires significantly larger bond dimension D to obtain the same quality of approximation. This point is the very important one as computational complexity scales in both cases as D6 . Additionally it can be shown that for Gibbs operator variational renormalization number of the most costly procedures in one iteration of the renormalization procedure scales logarithmically with β. For the imaginary time evolution number of such procedures in total scales at least linearly with β as number of renormalization steps scales linearly with β. Therefore if in a low temperature phase the self consistent update converges sufficiently fast (as it’s the case of low temperature phases of 2D quantum Ising model and 2D isotropic compass model) significantly smaller number of the most costly calculations is needed in th case of Gibbs operator variational renormalization, leading to further simulation acceleration. Furthermore in the case of low temperature phase laying below a critical point imaginary time evolution requires Gibbs state simulation close to the critical point contrary to Gibbs operator variational renormalization . Speed of convergence of CTMRG procedure (which is the most costly procedure for both methods) decreases with growing correlation length which at the critical point becomes infinite. Therefore another significant simulation time reduction can be obtained here. Superiority of Gibbs operator variational renormalization was demonstrated numerically in the quantum Ising model (see Fig. 21).. 30.

(31) Figure 20: Comparison of an order parameter hZi close to the critical point in 2D Ising model obtained by variational Gibbs operator renormalization (the green curve) with the results obtained by the imaginary time evolution for the same bond dimension D = 6. Results for the imaginary time evolution were obtained with (the red curve) and without (the black curve) a bias term changing the phase transition into a crossover introduced to speed up simulations. The figure is taken from the fourth article.. 31.

(32) 3 3.1. Benchmark results - summary. General remarks.. Presented above methods are unbiased because they do not make assumptions about physical nature of simulated states. They share this property with Quantum Monte Carlo methods. Unlike Monte Carlo methods they don’t suffer the sign problem for system with frustrations or fermionic systems. Their computational complexity depends primarily on the bond dimension D depending in turn on amount of entanglement in the simulated state. Quality of obtained results is controlled by their convergence in D. Here I present two examples of the convergence analysis. In the case of 2D quantum Ising model the order parameter close to the critical point is converged for D = 6 (see Fig. 21). For 2D isotropic compass model more sophisticated analysis is necessary. The order parameter and its susceptibility are not strictly convergence in D but regime in which they scale linearly with measure of convergence (introduced in an appendix to the fifth article) is obtained for D ≥ 8 (see Figs 22, 23).. Figure 21: Convergence of the order parameter hZi in the bond dimension D for 2D quantum Ising model close to the critical point. The results are converged for D = 6. The figure is taken from the fourth article.. Figure 22: Convergence of the order parameter Q for the compass model close to the critical point. Convergence measure ez was introduced in an appendix to the fifth article. The data convergence is approximated well by linear function making possible extrapolation to ez = 0 (the limit of the large bond dimension D) . The figure is taken from the fifth article. 32.

(33) Figure 23: Convergence of the order parameter susceptibility χ for the compass model close to the critical point. The data convergence is again approximated well by linear function. The figure is taken from the fifth article.. 33.

(34) 3.2. 2D quantum Ising model.. 2D quantum Ising model is a spin 1/2 lattice model. The model is defined by a Hamiltonian X X H=− Zi Zj − h Xi , <i,j>. (49). i. for a square lattice. Here Z, X are Pauli matrices. The model has a nearest-neighbor ferromagnetic coupling and a transverse magnetic field with strength h. Finite temperature second order phase transition from paramagnetic to ferromagnetic phase occurs for h ≤ h0 ≈ 3.04 [17]. The model is used as a standard benchmark for new numerical methods. Neighborhood of the critical point at h = 2/3 · h0 was simulated as a benchmark of all proposed methods (the results were presented in the first, the third and the fourth paper). The results are consistent with the second order phase transition. Close to the critical point states with very large correlation lengths ξ were obtained (up to ξ = 290) (see Fig. 24). Furthermore very close to the critical point critical behavior of the order parameter and correlators in intermediate distance range was obtained (see Fig. 25).. Figure 24: Ferromagnetic correlation function CR = hZx,y Zx+R,y i − hZx,y ihZx+R,y i in the ferromagnetic phase close to the critical point at h = 2/3h0 . Results were obtained for the bond dimension D = 6. The correlation length CR diverges as the critical point is approached. The largest correlation length obtained is ξ = 290. The figure is taken from the fourth article.. Figure 25: In (A), the order parameter (spontaneous magnetization) hZi near the critical point. The results for the largest simulated bond dimension D = 6 (the red line) can be fitted well by critical 0 behavior hZi ∝ (β − βc )β (the dashed blue line). In (B) the ferromagnetic correlation function CR = hZx,y Zx+R,y i − hZx,y ihZx+R,y i at βc . The results for the largest simulated bond dimension D = 6 (the red line) again can be fitted well by critical behavior CR ∝ R−η (the dashed brown line) in intermediate distance range 1 ≤ R ≤ 30 . The figure is taken from the fourth article. The results were obtained at h = 2/3h0 .. 34.

(35) 3.3. Isotropic 2D quantum compass model - frustrations of interactions.. 2D quantum compass model is a lattice spin 1/2 model given by a Hamiltonian H=−. 1+AX 1−AX Xj Xj+ea − Zj Zj+eb , 4 4 j j. (50). for square lattice. Here X,Z are Pauli matrices, and ea (eb ) is unit vectors along horizontal (vertical) lattice axis. The model has ferromagnetic nearest-neighbor XX coupling along horizontal bonds and ferromagnetic nearest-neighbor ZZ coupling along vertical bonds. The model is a canonical example of frustrations of interactions. A second order phase transition to a low temperature nematic phase occurs in the model for the isotropic case A = 0. The nematic order parameter is [18] Q = hXj Xj+ea i − hZj Zj+eb i.. (51). By adding small anisotropy one can obtain the order parameter susceptibility χ=. dQ

(36)

(37) .

(38) dA A=0. (52). The nematic phase is described as correlated Ising chains. The sign error free Quantum Monte Carlo (QMC) simulations of the model can be performed for proper choice of the boundary conditions (e. g. twisted or periodic boundary conditions) [16]. They were used to estimate the critical temperature TC . Nevertheless estimation of TC with QMC turned to be a challenging task because of anomalous scalings occurring for intermediate cluster sizes [16]. It was claimed that twisted boundary conditions removes the problem [16] giving correct estimate TC = 0.0585(3). To estimate TC the order parameter Q and the order parameter susceptibility χ in the symmetric phase were simulated close to the critical point. The obtained results can be fitted well by critical behavior (see Figs 26, 27). These fits were used to obtain two estimates of TC . The order parameter fit gives TC = 0.0609 and the susceptibility fit gives TC = 0.0602. Both estimates were combined into an estimate Tc = 0.0606. Mine estimate differ from the QMC estimate by 3.5%. Accuracy of this estimate is limited by extrapolations to limit of large bond dimension D (see sub-chapter 3.1) and β − βc = 0. On the other side accuracy of the QMC estimate is limited by extrapolation to limit of infinite system. In the fifth paper it was concluded that regardless what source of the discrepancy is both estimates agree well. Furthermore possibility of correlator simulation was demonstrated (see Fig. 28). As a result correlation lengths were estimated at β = 17.2 (ξ = 40(2) along the chains and ξ = 6.9(4) across the chains). Phase diagrams of models with frustrations become controversial when unbiased results are not available because of the QMC sign problem. This benchmark suggest that the introduced method can provide such results for frustrated spin models. Work on its application to the Kitaev-Heisenberg model is in progress [21].. 35.

(39) 0. Figure 26: The order parameter Q (the points) fitted with critical behavior Q = Q0 (Tc − T )β (the blue curve) and its error bars (the red and green curves). The error bars estimate uncertainty of extrapolation to limit of the large D (see sub-chapter 3.1). Panel (B) is a log-log plot of Q as a function of distance from the fitted critical βc . The figure is taken from the fifth article.. Figure 27: The order parameter susceptibility χ (the points) fitted with critical behavior χ = χ0 (T − Tc )γ (the blue curve) and its error bars (the red and green curves). Panel (B) is a log-log plot of χ as a function of distance from the fitted critical βc . The figure is taken from the fifth article.. 36.

(40) Figure 28: Correlators in the nematic phase of the compass model simulated with D = 15 at β = 17.2. In panel (A) hXm Xm+d i correlation function along axis parallel to the chains. In panel (B) hXm Xm+d i correlation function along axis orthogonal to the chains. Both plots are semi-log plots. The figure is taken from the fifth article.. 3.4. Fermionic models.. The proposed methods of Gibbs operator simulation can be generalized to fermionic lattice models. Their computational complexity is than the same function of the PEPO bond dimension D as in the case of spin models. Such generalization of purification imaginary time evolution was proposed in the second article. The Gibbs operator variational renormalization also can be properly generalized. Its generalized fermionic version was briefly sketched in an appendix to the fifth paper. In the second paper a second order phase transition in a model of interacting spinless fermions X X H=− (ci c†j + cj c†i ) − g (ni − 1/2)(nj − 1/2), (53) <i,j>. <i,j>. (here ci is a fermionic annihilation operator on site i and ni = c†i ci is an occupation number operator on site i) was simulated as a proof of principle of the proposed algorithm (see Fig. 29). Also benchmark results for a non-interacting spinless fermion model X X H=− (ci c†j + cj c†i ) − δ ci cj + c†j c†i (54) <i,j>. <i,j>. were presented in the second paper. This model can be solved analytically. Ground state of the model doesn’t obey the area law (9) so at low temperatures its simulation becomes challenging. Comparison of obtained results with the analytical solution was preformed. The results from the paper are outdated. Here I present unpublished benchmark results (for δ = 0) obtained by the Gibbs operator variational renormalization (see Fig. 30).. 37.

(41) 1 0,9. n. 0,8 0,7 -3. b=10 -4 b=10 -5 b=10. 0,6 0,5. 0,15. 0,2. 0.175. 0.225. β Figure 29: The fermion density n close to the second order phase transition in the interacting spinless fermion modelPwith g = 10. The density n converges to critical behavior as a small bias term Hbias = −2gb i ni changing the phase transition into a crossover goes to zero. The figure is taken from the second article.. 0.2. hc†i cj i. 0.198 0.196 analytical D=8 D = 10 D = 12. 0.194 0.192 2.5. 3. 3.5. 4.. 4.5. β Figure 30: Nearest neighbor hopping hc†i cj i for the free fermion model with δ = 0. Results obtained by the variational Gibbs operator renormalization are compared with the analytical solution. The results converge fast to the analytical solution with the growing bond dimension D.. 38.

(42) 3.5. The Hubbard model.. I present here also the benchmark results for the Hubbard model. The model is given by square lattice Hamiltonian: X X X H = −t ni,↑ ni,↓ + µ ni , (55) c†i,σ cj,σ + h.c. + U i. <i,j>,σ. i. where σ is a spin orientation (σ =↑, ↓), ci,σ is annihilation operator for an electron with spin σ on site i, ni,σ is occupation number operator for an electron with spin σ on site i and ni = ni,↑ + ni,↓ . µ is chemical potential used to control electron density. The model is a candidate for high-Tc superconductivity model, which is expected to appear close to U/t = 8 and n = 0.875. The model was simulated at high T by Determinant Quantum Monte Carlo (DQMC) [24] and by Dynamical Cluster Approximation (DCA) [25]. DCA is a variant of cluster DMFT methods using QMC impurity solver. In both cases the sign problem appears becoming more severe with growing β and cluster size N . It prevents both methods from obtaining conclusive results for T /t < 0.25 close to U/t = 8 and n = 0.875 [24, 25]. For high T the model was also simulated by high temperature series expansion methods which example is Numerically Linked Cluster Expansion (NLCE) [26]. Strongly correlated case of U/t = 8 was simulated. In figure 31 results at two temperatures T /t = 1, 0.5 and electron densities n = 1, 0.875, 0.85 are collected and compared with the most recent dynamical cluster approximation (DCA) results from Ref. [28, 27]. Both methods have refinement parameters: the bond dimension D in the tensor network and cluster size N in DCA. The DCA results can be extrapolated with N → ∞. The extrapolated values are marked in the figure with the dashed horizontal lines. For the simulated temperatures DCA sign problem is at most mild [27]. Furthermore DCA results for n = 1 and n = 0.85, T /t = 1 were shown to agree with NLCE and DQMC results [28, 27], so they are reliable reference points. With increasing bond dimension D mine results approach DCA results. For the largest simulated D = 20 the energies (the double occupancies) differ from the DCA extrapolated results by a relative error smaller than 1.9% (3.6%).. 39.

(43) Figure 31: Results for energy E and double occupancy hni,↑ ni,↓ i from variational Gibbs operator renormalization (circles) versus extrapolated to N = ∞ results from dynamical cluster approximation (lines) for U/t = 8. D ≤ 20 is the bond dimension of the PEPO.. 40.

(44) 4. Conclusions.. Novel method of Gibbs operator simulation for 2D infinite lattice models is presented in the thesis. It is based on variational renormalization of the Gibbs operator tensor network representation obtained by Suzuki-Trotter decomposition. After the renormalization numerically tractable approximation of the Gibbs operator by a Pair Entangled Projected Operator (PEPO) is obtained preserving the area law for mutual information. The presented numerical method is unbiased like Quantum Monte Carlo (QMC) methods and doesn’t suffer the sign problem unlike QMC. The method was benchmarked for a model with frustrations of interactions (2D isotropic compass model) for which QMC results without the sign problem exist. It was also benchmarked for an analytically solvable free fermion model challenging to tensor network methods. A benchmark for strongly correlated fermionic system was performed in the Hubbard model by comparison with more standard unbiased methods. The results obtained in all cases agree well with the reference points. The benchmark results encourage the method application to models with the sign problem and its further development. The first article introduces the simplified version of outdated purification imaginary time evolution. The second generalizes it to fermionic systems. The second generalizes it to fermionic systems. The third introduces improved self-consistent and faster version of purification imaginary time evolution. The fourth introduces Gibbs operator variational renormalization which is superior to all variants of purification imaginary time evolution. The fifth benchmarks the new method for 2D isotropic compass model.. 41.

(45) 5. References. [1] U. Schollwöck, Ann. of Phys. 326, 96 (2011). [2] F. Verstraete and J. I. Cirac, arXiv:cond-mat/0407066 (2004); [3] J. Eisert, M. Cramer, M.B. Plenio, Rev. Mod. Phys. 82, 277 (2010); [4] R. Orus, Ann. of Phys. 349, 117 (2014). [5] F. Verstraete, J.I. Cirac, V. Murg, Adv. Phys. 57,143 (2008). [6] N. Schuch, M. M. Wolf, F. Verstraete, J. I. Cirac, Phys. Rev. Lett. 100, 070502 (2008). [7] A. Molnár, N. Schuch, F. Verstraete, J. I. Cirac, Phys. Rev. B 91, 045138 (2015). [8] M. B. Hastings, Phys. Rev. B 73, 085115 (2006). [9] J. Jordan, R. Orús, G. Vidal, F. Verstraete, and J. I. Cirac, Phys. Rev. Lett. 101, 250602 (2008); [10] P. Corboz i G. Vidal, Phys. Rev. B 80, 165129 (2009); P. Corboz, R. Orús, B. Bauer, i G.Vidal, Phys. Rev. B 81, 165104 (2010). [11] P. Corboz, S. R. White, G. Vidal, and M. Troyer, Phys. Rev. B 84, 041108 (2011); P. Corboz, T. M. Rice, M. Troyer, Phys. Rev. Lett. 113, 046402 (2014). [12] P. Corboz, Phys. Rev. B 93, 045116 (2016). [13] B. Keimer, S. A. Kivelson, M. R. Norman, S. Uchida and J. Zaanen, Nature 518, 179 (2015). [14] R. Orus, A. C. Doherty, G. Vidal; Phys. Rev. Lett. 102, 077203 (2009). [15] Y. H. Matsuda, N. Abe, S. Takeyama, H. Kageyama, P. Corboz, A. Honecker, S. R. Manmana, G. R. Foltin, K. P. Schmidt, and F. Mila, Phys. Rev. Lett. 111, 137204 (2013); P. Corboz, F. Mila, Phys. Rev. Lett. 113, 147203 (2014). [16] S. Wenzel, W. Janke, A. M. Läuchli, Phys. Rev. E 81, 066702 (2010). [17] H. Rieger and N. Kawashima, Eur. Phys. J. B 9, 233 (1999); H. W. J. Blote and Y. Deng, Phys. Rev. E 66, 066110 (2002). [18] Z. Nussinov and J. van den Brink, Rev. Mod. Phys. 87, 1 (2015). [19] R. Orús and G. Vidal, Phys. Rev. B 80, 094403 (2009); [20] P. Czarnik, M. Rams and J. Dziarmaga, in preparation. [21] A. Francuz and J. Dziarmaga, in preparation. [22] T. Nishino and K. Okunishi, J. Phys. Soc. Jpn. 67, 3066 (1998); [23] R. Orús and G. Vidal, Phys. Rev. B 80, 094403 (2009). [24] R. T. Scalettar, G. C. Batrouni, Understanding Quantum Phase Transitions, ed. L. Carr, (CRC Press, Boca Raton, 2010) 265. [25] D. J. Scalapino, Handbook of High Temperature Superconductivity, ed. J. R. Schriefer, (Springer, 2007). [26] E. Khatami and M. Rigol, Rhys. Rev. A 84 053611 (2011). [27] J. P. F. LeBlanc and E. Gull, Rhys. Rev. B 88 155108 (2013). [28] J. P. F. LeBlanc et al., Phys. Rev. X 5, 041041 (2015) .. 42.

(46) 6. Articles.. 43.

(47) PHYSICAL REVIEW B 86, 245101 (2012). Projected entangled pair states at finite temperature: Imaginary time evolution with ancillas Piotr Czarnik,1 Lukasz Cincio,2 and Jacek Dziarmaga1 1. Instytut Fizyki Uniwersytetu Jagiello´nskiego and Centre for Complex Systems Research, ul. Reymonta 4, 30-059 Krak´ow, Poland 2 Perimeter Institute for Theoretical Physics, Waterloo, Ontario, Canada N2L 2Y5 (Received 5 September 2012; revised manuscript received 13 November 2012; published 5 December 2012) A projected entangled pair state (PEPS) with ancillas is evolved in imaginary time. This tensor network represents a thermal state of a two-dimensional (2D) lattice quantum system. A finite-temperature phase diagram of the 2D quantum Ising model in a transverse field is obtained as a benchmark application. DOI: 10.1103/PhysRevB.86.245101. PACS number(s): 75.10.Jm, 03.65.Ud, 03.67.−a, 05.30.Fk. I. INTRODUCTION. Quantum tensor networks are a competitive tool to study strongly correlated quantum systems on a lattice. They originate from the density matrix renormalization group1 — an algorithm to minimize the energy of a matrix product state (MPS) ansatz in one dimension. In recent years, the MPS was generalized to a two-dimensional (2D) “tensor product state” better known as a projected entangled pair state (PEPS).2 Another type of tensor network is the multiscale entanglement renormalization ansatz (MERA)3 that is, in some respects, a refined version of the real space renormalization group. Both PEPS and MERA can be applied to strongly correlated fermions in two dimensions4 because they do not suffer from the notorious fermionic sign problem. This makes them a powerful tool to attack some of the hardest problems in strongly correlated electronic systems, including the enigmatic high-temperature superconductivity.5 Indeed, PEPS has already provided initial results for the ground-state energy of the t-J model6 that can compete with the best variational Monte Carlo results.7 In contrast to the ground state, thermal states have been explored mainly with the MPS in one dimension,8,9 but they are more interesting in two dimensions, where they can undergo finite-temperature phase transitions. In two dimensions, thermal states were represented by tensor product states and contracted with the help of the higher-order singular value decomposition in Ref. 10. A similar projected entangled-pair operator (PEPO) ansatz was proposed in Ref. 11. In this paper, we follow a different route. In a way that can be easily generalized to two dimensions, the MPS can be extended to finite temperature by appending each lattice site with an ancilla.8 A thermal state is obtained by imaginary time evolution of a pure state in the enlarged Hilbert space, starting from infinite temperature. Unfortunately, in contrast to one dimension, where the time evolution of a MPS can be simulated accurately and efficiently, in two dimensions the time evolution of PEPS appears to be a hard problem. It requires accurate computation of a tensor environment that is often hard to approximate accurately and reliably. The aim of this work is to overcome this problem.. The enlarged Hilbert space is spanned by states s |is ,as , where the product runs over lattice sites s. The state of spins at infinite temperature, ρ(β = 0) = s ( S1 S−1 i=0 |is is |) ∝ 1, is obtained from a pure state in the enlarged Hilbert space,. where. 1098-0121/2012/86(24)/245101(6). (1). S−1 1 |ψ(0) = √ |is ,ia S s i=0. (2). is a product of maximally entangled states of every spin with its ancilla. The state ρ(β) ∝ e−βH at finite β is obtained from |ψ(β) ∝ e− 2 βH |ψ(0) ≡ U (β)|ψ(0) 1. after imaginary time evolution for time β with. (3). 1 H. 2. III. PEPS. In the quantum Ising model with spin- 12 that we consider in the rest of this paper, the translational invariance is not broken and a unit cell encloses only one lattice site. Therefore, for an efficient simulation of the time evolution, we represent |ψ(β) by a translationally invariant PEPS with the same tensor Aia trbl (β) at every site. Here i,a = 0, . . . ,S − 1 are the spin and ancilla indices, respectively, S = 2, and t,r,b,l = 0, . . . ,D − 1 are the bond indices to contract the tensor with similar tensors at the nearest-neighbor sites; see Fig. 1(a). The ansatz is |ψ(β) = A [{is ,as }] |is ,as ≡ |ψA . (4) {is ,as }. s. Here the sum runs over all indices is ,as at all sites s. The amplitude A is the tensor contraction in Fig. 1(b). The initial state (2) can be represented by a tensor ia Aia trbl = δ δt0 δr0 δb0 δl0 .. (5). D = 1 is the minimal bond dimension sufficient to represent the initial state. IV. QUANTUM ISING MODEL IN TWO DIMENSIONS. We proceed with H=− Zs Zs − h Xs ≡ HZZ + HX .. II. THERMAL STATES. We consider spins on an infinite square lattice with a Hamiltonian H. Every spin has S states i = 0, . . . ,S − 1 and is accompanied by an ancilla with states a = 0, . . . ,S − 1.. ρ(0) = trancillas |ψ(0)ψ(0)|,. s,s . (6). s. Here Z,X are Pauli matrices. The model has a ferromagnetic phase with nonzero spontaneous magnetization Z for small. 245101-1. ©2012 American Physical Society.