A simplified approach for 3D tensor network simulations

(1)

MSc Physics & Astronomy, Track theoretical

physics

Master Thesis

A simplified approach for 3D tensor

network simulations

Author:

Supervisor:

P.C.G. Vlaar

dr. P.R. Corboz

Student number:

Second examiner:

10218521

dr. J. van Wezel

Number of credits:

Examination date:

60 ECTS

August 25, 2018

Institute for Theoretical Physics, Amsterdam

(2)

Title: A simplified approach for 3D tensor network simulations Author: P.C.G. Vlaar, p.c.g.vlaar@uva.nl, 10218521

Supervisor: dr. P.R. Corboz Second examiner: dr. J. van Wezel Examination date: August 25, 2018

Institute for Theoretical Physics University of Amsterdam

Science Park 904, 1098 XG Amsterdam uva.nl

(3)

3

Abstract

Tensor network methods have been very successful to accurately simulate one- and two-dimensional classical and quantum models, with White’s density-matrix renormalization group (DMRG) as a famous example [1]. A way to generalize tensor network methods to three dimensions would be highly desirable, especially for simulating fermionic and frustrated models which in many cases are hard or impossible to solve using other methods. Here a novel way of contracting three-dimensional tensor networks using clusters is introduced. These clusters consist of a certain number of tensors which are contracted exactly or with a higher accuracy as compared to the environment of the cluster, which is approximated effectively. The method is implemented on the square and cubic lattices and applied to the quantum Ising model with a transverse field and the Heisenberg model. The results are compared to corner transfer matrix (CTM) results and other reference values. Also, as a non-trivial example, a simulation is made for the SU(3) Heisenberg model. A three-sublattice striped state is found to be lower in energy than a two-sublattice striped state.

(4)

(5)

5

Populairwetenschappelijke samenvatting

Hoewel er veel bekend is over de fundamentele deeltjes waaruit materialen om ons heen zijn opgebouwd en de krachten die tussen deze deeltjes werken, blijkt het zeer moeilijk te zijn om met deze relatief simpele regels te bepalen hoe systemen die bestaan uit veel deeltjes zich gedragen. Dit is een soortgelijk probleem als het gedrag van mieren in een kolonie. Ondanks het feit dat een kolonie is opgebouwd uit simpele eenheden, de individuele mieren, kan het gedrag van de kolonie als geheel zeer complex zijn. Om de natuurkundige fenomenen die we zien toch te kunnen beschrijven wordt vaak gebruik gemaakt van versimpelde modellen waarin alleen de belangrijkste interacties tussen de deeltjes worden meegenomen. Vaak zijn zelfs deze versimpelde modellen echter te moeilijk om op te lossen.

Sinds de ontwikkeling van de computer zijn er veel algoritmen ontwikkeld die benaderingen kunnen geven van de oplossingen van een aantal van deze modellen. Eén groep van deze algoritmes is gebaseerd op zogenaamde tensor netwerken. Deze technieken zijn zeer succesvol om één dimen-sionale (waarbij de deeltjes zich op een lijn bevinden, zoals in een draad) en twee dimendimen-sionale (waarbij de deeltjes zich in een vlak bevinden) modellen te kunnen beschrijven. Hoewel veel van deze modellen interessant natuurkundig gedrag vertonen en ons veel kunnen ler en, staan drie dimensionale modellen dichter bij onze alledaagse realiteit. De huidige tensor netwerk technieken hebben echter te veel extra rekenkracht en geheugen nodig om deze direct toe te kunnen passen op deze modellen. In deze scriptie wordt er een nieuwe techniek ge¨ıntroduceerd waarmee deze simu-laties toch gedaan kunnen worden en er wordt gekeken naar de nauwkeurigheid van de resultaten ten opzichte van andere methoden.

(6)

(7)

Introduction

The 20th _{century has seen several important paradigm shifts in physics, among which is the}

discovery of the quantum world at the submicroscopic scale. Many breakthrough achievements have been made in identifying the fundamental particles that make up the universe and the interactions between them, with the standard model of particle physics as an important milestone. Despite the details of our understanding of these fundamental interactions, applying these laws to describe the behavior of macroscopic numbers of particles has turned out to be highly non-trivial. Many relatively simple models are still poorly understood to this day.

The development of computer technology has given us a new way of gaining more insight into these problems. However, exact evaluations, e.g. by diagonalizing the Hamiltonian, suffer from an exponential increase in the number of coefficients with the number of particles in the system. Therefore, it is necessary to use algorithms which give estimations. One of the main groups of algorithms is based on Monte Carlo (MC) sampling. These methods have been highly successful for bosonic and unfrustrated systems. Unfortunately, studies for fermionic and frustrated models suffer from the infamous negative sign problem, which significantly hinders the study of these models.

A breakthrough in the study of one-dimensional quantum systems was achieved by the density-matrix renormalization group algorithm (DMRG) [1]. The method is based on ordering the states in the exponentially large Hilbert space in such a way that only the states which are relevant for the low-energy physics are kept. Because of the success of the DMRG algorithm, and its underlying ansatz the matrix product state (MPS) [2, 3, 4, 5], generalizations of this method have been made to study two-dimensional systems, e.g. based on the projected entangled pair state (PEPS) [6]. A shared property of these techniques is that they are based on the contraction of networks of tensors, therefore they are collectively referred to as tensor network (TN) techniques. One of the advantages of these methods is that they do not suffer from the negative sign problem and therefore allow for the numerical study of fermionic and frustrated systems. Although these methods have been developed to study quantum many-body problems, connections are also found to other fields, for example quantum gravity [7].

Unlike MC methods, algorithms based on TN suffer from a steep growth in the computational cost when moving to higher dimensional networks. As most of the physical systems around us are three-dimensional, an extension of TN methods to three-dimensions would be highly desirable, especially for the study of fermionic and frustrated models.

A well-known property of quantum systems is that quantum fluctuations1 _{generally become}

less strong for higher-dimensional systems. One way in which this can be seen is by considering

1_{There are two types of fluctuations: quantum and thermal/classical fluctuations. Whereas thermal}

fluctua-tions only occur at non-zero temperatures, quantum fluctuafluctua-tions, which are caused by the Heisenberg uncertainty principle, also play a role at T = 0 [8]. In this work we will only study models at T = 0.

(10)

the monogamy property of entanglement. Consider a system consisting of three spin-1₂ particles A, B and C. If particles A and B are in a singlet state, and are therefore maximally entangled, A and B cannot be entangled with particle C. In a similar way a partial entanglement between particles A and B implies that there is a maximum on the entanglement that particles A or B can share with particle C [9]. Osborne & Verstraete showed that the same also holds for systems of n spin-1₂ particles [10]. The monogamy property implies that for higher-dimensional systems, where the entanglement must be shared among more particles, the entanglement between two individual particles will generally be smaller. Furthermore, for systems with a dimension larger than the so-called upper-critical dimension, their critical behavior does not depend on fluctuations anymore and they can be exactly described by mean-field theory [8]. In TN the accuracy with which quantum fluctuations are taken into account can be tuned by a parameter called the bond dimension D, which will be introduced more precisely later. Although the computational cost increases strongly for three-dimensional TN, the reduced importance of quantum fluctuations should make it possible to treat models using a lower D than their one- and two-dimensional counterparts. Although some studies have been done on three-dimensional TN, e.g. by Refs. [11, 12, 13, 14, 15], it remains a relatively unexplored area.

In this work a new cluster method will be presented to obtain estimations of expectation values for three-dimensional quantum systems. The method is based on contracting a certain number of tensors (i.e. a cluster) around the tensor(s) to which the operator is applied exactly or with a higher accuracy as compared to the environment, which is considered effectively. The approach is simple to implement and has a relatively small computational cost, especially for the simplest types of clusters. The accuracy of the method is evaluated for the quantum Ising model with a transverse field and the Heisenberg model in both two- and three-dimensions.

In chapter 2 a review is given of several important concepts and algorithms for the theory of TN and several approximate contraction approaches for two-dimensional TN are discussed. This chapter concludes with a discussion of the cluster method that had been explored in this work. To clarify several aspects of the algorithms used as well as some limitations of them, chapter 3 gives some results for the one-dimensional Ising chain. In the next two chapters, 4 and 5, the cluster method is used for the Ising model with a transverse field and the Heisenberg model in two-and three-dimensions respectively, which are compared to the established corner-transfer matrix (CTM) contraction in 2D and several results obtained from other methods. A simulation of the SU(3) Heisenberg model, which cannot be evaluated using MC simulations, is presented next in chapter 6. A final discussion and possible future directions are given in chapter 7.

(11)

Chapter 2

Introduction to tensor networks

2.1 The matrix product state

The main challenge in computational many-body physics is how to efficiently perform simulations of large numbers of particles. To clarify this problem let us consider a spin-1₂ system on a lattice with N sites. In its most general form, the quantum state can be written as

|ψi = X

i1,i2,...,iN

Ci1,i2,...,iN|i1i ⊗ |i2i ⊗ . . . ⊗ |iNi (2.1)

where each |ini represents a two-dimensional local basis at each site. The number of coefficients

in Ci1,i2,...,iN grows exponentially with the number of sites as 2

N_{. Because of this rapid growth}

techniques like exact diagonalization of the Hamiltonian can in practice only be performed for relatively small system sizes of up-to several dozen sites. Macroscopic systems, on the other hand, consist of numbers of particles of the order of Avogadro’s number. Therefore, this way of representing the state clearly seems inefficient.

One observation which we can make is that all physical states only live on a tiny fraction of the total Hilbert space and can only reach a small section of it during their lifetime [4]. Furthermore, it can be shown that low-energy and ground states, which are of particular interest in condensed matter physics, occupy an even smaller section [4]. An efficient algorithm therefore only has to target this subsection of the Hilbert space, while omitting the degrees of freedom corresponding to unphysical or high-energy states.

To achieve such an algorithm, it is useful to rewrite the wave function in equation 2.1 to a matrix product state (MPS) [2, 4, 5]. In an MPS, matrices are assigned to each site in the lattice such that the product of these matrices describes the state

|ψi = X i1,...,iN α1,...,αN −1 Ai1 α1B i2 α1α2. . . P iN −1 αN −2αN −1Q iN αN −1|i1i ⊗ |i2i ⊗ . . . ⊗ |iN −1i ⊗ |iNi (2.2)

in which i1, . . . , iN correspond to the indices of the different sites and α1, . . . , αN −1represent the

connections between the different matrices and they are summed over. These are also referred to as the physical and virtual indices respectively. The virtual index runs from 1 to a value D, which is called the bond dimension. Note that by rewriting the state to an MPS we went from an exponential scaling with the number of sites to a polynomial scaling. In an exact MPS representation of a state, however, D will grow rapidly with the number of sites, such that we encounter a similar problem as the one that we started with. Luckily, it turns out that an efficient approximation of the state can be made by truncating D to a relatively small value.

(12)

(a) (b) (c)

....

L

(d) s vi Mij Ti1i2...iL

Figure 2.1: Graphical representation for tensor networks. Each shape represents a tensor and each leg an index of the tensor. (a) represents a scalar, (b) represents a vector, (c) represents a matrix and (d) represents a rank-L tensor.

The most famous example of an algorithm which makes use of the MPS representation is the density-matrix renormalization group (DMRG) [1]. Although DMRG was originally formulated using density matrices, the realization was later made that it makes use of an underlying MPS ansatz [3].

In the next section a useful diagrammatic notation for tensor networks will be introduced. After that we discuss the procedure to obtain an MPS from the general quantum state in equation 2.1, several properties of them and the way in which they can be used to e.g. determine expectation values. Also, the questions of why and when the MPS ansatz gives us an efficient representation of a quantum state will be elaborated upon.

2.1.1 Diagrammatic notation

A multiplication of many tensors tends to become complex quickly in the standard mathemat-ical notation. Fortunately, an alternative notation exists making use of diagrams which will be introduced here. The notation is loosely based on the tensor diagram notation introduced by Penrose [16].

A tensor network (TN) consists of a certain number of tensors that are contracted in a certain way. In the diagrammatic notation each tensor is represented by a shape with a number of legs, where each leg represents one of the indices of the tensor. Figure 2.1a-d give examples of a scalar, a vector, a matrix and a rank-L tensor respectively. These types of shapes will form the building blocks of our TNs. Contractions over one or more indices of two tensors is represented by a line between them. One of the simplest kinds of contractions we can perform is multiplying two matrices A and B together

Cαγ =

X

β

AαβBβγ (2.3)

where the contraction is performed over index β. The diagram for this contraction is given in figure 2.2a. The other indices, α and γ, which are not contracted over, are called open indices.

Of course, we can also draw more complicated diagrams, and this is where the notation becomes very useful. While a contraction as

F = X

α,β,γ,δ,,ζ

AαβBαγδCγDζEδζβ (2.4)

looks relatively complicated, it can be drawn as a simple diagram as depicted in figure 2.2b. The diagrammatic notation makes it much more transparent to see which contractions are performed. From a computational viewpoint, an important observation that should be made is that the order in which these contractions are performed often matters. As an example, we can take the TN of equation 2.4. For simplicity we assume all contractions have the same dimension D. Figure 2.3 shows two contraction sequences of this TN. Clearly the contraction sequence of order O(D4) is, in

(13)

2.1. THE MATRIX PRODUCT STATE 13 α β γ A B (a) A B E C D δ α β γ ζ (b)

Figure 2.2: Examples of contractions. (a) shows the matrix contraction of equation 2.3 and (b) shows the slightly more complicated contraction of equation 2.4 in the diagrammatic notation.

O(D6₎ O (D₄ ) O(D4₎ O(D 4) O(D3₎ _O(D2₎

Figure 2.3: Different contraction orders often result in different costs of contraction. The upper contraction sequence has a cost of order O(D6_{), while the lower contraction sequence only has a cost of order O(D}4_).

terms of computational cost, preferable to the sequence of order O(D6). When implementing TN algorithms, defining the computationally cheapest contraction sequence is therefore an important consideration to make as it can have a significant impact on the efficiency of the algorithm. In this work extensive use has been made of the Netcon and NCON functions developed in Refs. [17] and [18] which can determine the optimal contraction order and perform the contractions of our tensor networks respectively.

2.1.2 Obtaining the MPS representation

Let us go back to the problem of finding an MPS from a general quantum state. This section will closely follow the discussion given in Ref. [5]. The key manipulation that is necessary to obtain the MPS form is the singular value decomposition (SVD) from linear algebra. By making use of the SVD any rectangular matrix M of dimensions NA× NB can also be expressed as

M = U λV† (2.5)

in which U and V are two matrices with orthonormal columns of dimensions NA× min (NA, NB)

and NB × min (NA, NB) respectively and λ is a real positive diagonal matrix of dimensions

min (NA, NB) × min (NA, NB). The entries in λ are also called singular values. We will order

the singular values in λ in such a way that they are decreasing, i.e. λ1> λ2> . . . > λN.

To obtain the MPS form we will iteratively split off one of the physical legs of Ci1,i2,...,iN from

(14)

into a matrix Mi1,(i2...iN). The SVD gives Mi1,(i2...iN)= X α1 Ui1,α1λ [1] α1(V †₎ α1,(i2...iN). (2.6)

To split off the next leg we reorganize λ[1]α1(V

†₎

α1,(i2...iN) = M(α1i2),(i3...iN) and perform another

SVD Ci1,i2,...,iN = X α1 Ui1,α1M(α1i2),(i3...iN)= X α1,α2 Ui1,α1U(α1i2),α2λ [2] α2(V †₎ α2,(i3...iN). (2.7)

This process of splitting off one of the legs is continued until only one physical leg remains on the matrix on the right. In the end we obtain

Ci1,i2,...,iN = X α1,α2,...,αN −1 Ui1,α1U(α1i2),α2. . . U(αN −2iN −1),αN −1λ [N −1] αN −1(V †₎ αN −1,iN = X α1,α2,...,αN −1 A[1]i1 α1 A [2]i2 α1,α2. . . A [N −1]iN −1 αN −2,αN −1λ [N −1] αN −1(V †₎ αN −1,iN (2.8)

where in the last step the U matrices were reshaped to rank-3 tensors (except for the very first U matrix which was reshaped to a rank-2 tensor).

For our further discussion it will turn out to be useful to rewrite equation 2.8 to a form where the singular value spectra are reintroduced on the bonds, a notation introduced by Vidal [19]. This can be done by decomposing each A[l]il

αl−1,αl= λ

[l−1] αl−1Γ

[l]il

αl−1,αl, in this way equation 2.8 can be

rewritten as Ci1,i2,...,iN = X α1,α2,...,αN −1 Γ[1]i1 α1 λ [1] α1Γ [2]i2 α1,α2λ [2] α2. . . λ [N −1] αN −1Γ [N ]iN αN −1 (2.9)

where we removed the dummy matrix λ[0]α0 and we also replaced (V

†₎iN

αN −1 = Γ

[N ]iN

αN −1. The full

procedure can also be displayed in the diagrammatic notation as shown in figure 2.4.

Now let us perform the SVD procedure up to site n and pick the singular value matrix between site n and n + 1 to rewrite our state as

where we have collected all terms to the left of λ[n]αn in |uαniL and all terms to the right of λ [n] αn

in |vαniR. Note that the right-hand side of equation 2.10 is the Schmidt decomposition of the

state, which has several useful properties. One of them is that when we assume λ[n]αn consists of

D0 singular values, the state can be approximated up to D < D0 by

Because of the ordering we applied to the singular values before, this corresponds to discarding the smallest singular values. It can be shown that the difference between the original and approximate state

|ψi − | ˜ψi

corresponds to the Frobenius norm, i.e. it gives the optimal approximation of the state that can be made for a certain D [5]. This approximation will be essential for our further work as it will allow us to significantly reduce the computational efforts that need to be made in order to still find a good approximation to a quantum state in many cases.

(15)

2.1. THE MATRIX PRODUCT STATE 15 . . . i1 i2 iN −1iN C . . . α1 i1 i2 iN −1iN A[1] C . . . . . . α1 α2 i1 i2 i3 iN −1iN A[1] A[2] C . . . α1 α2 αN −2 αN −1 i1 i2 iN −1 iN A[1] _A[2] _A[N −1] _λ[N −1]_V† . . . i1 i2 iN −1 iN Γ[1] λ[1] Γ [2] λ[2] λ[N −2]Γ [N −1] λ[N −1]Γ [N ]

Figure 2.4: A general quantum state can be decomposed into an MPS by iteratively splitting off the physical legs by performing SVDs. At the final step, the singular value matrices are put back on the bonds, such that they are more easily accessible.

. . . . . . A[1] (A[1]₎† A[2] (A[2]₎† A[N −1] (A[N −1]₎† λ[N −1]V† λ[N −1]∗_V . . . . . . A[2] (A[2]₎† A[N −1] (A[N −1]₎† λ[N −1]V† λ[N −1]∗_V . . . A[N −1] (A[N −1]₎† λ[N −1]V† λ[N −1]∗_V λ[N −1] V† λ[N −1]∗ _V λ[N −1] λ[N −1]∗

Figure 2.5: Determining the norm of a left-canonical MPS. When contracted from the left the product forms unitaries from the left side in such a way that only the product of the singular value matrices on the last bond remains. If the MPS is normalized this product should be equal to one.

We have now described how an MPS form can be obtained from any quantum state of the form of equation 2.1, however this is not how an MPS is obtained in practice. For the problems considered in this work an MPS ansatz is assumed, consisting of random matrices, which is then brought to the ground state. Several algorithms exist to do this, which will be discussed in section 2.3.

2.1.3 The canonical form

The particular form of the MPS we obtained in equation 2.8 has an appealing property. To see this, let us calculate its norm. To determine the norm, we have to multiply the MPS by its complex conjugate, as depicted in figure 2.5. The property of orthonormality for the U and V matrices implies that U†_{U = I and V}†_{V = I. In our calculation this means that the tensors A}[1]

up to A[N −1] _{will form unitaries with their complex conjugates when contracted from the left,}

such that only a relatively simple contraction remains. This particular form of the MPS is called left-canonical. The V†V pair on the right vanishes in a similar way, such that only the two singular value matrices are left.

(16)

O A[1] M[2] B[3] B[4] (A[1]₎†_(M[2]₎†_(B[3]₎†_(B[4]₎† O M[2] (M[2]₎†

Figure 2.6: Determining the expectation value of a local operator of an MPS representing four sites in mixed canonical form around site two.

explained in the previous section in reverse from the final physical leg Ci1i2...iN = X αN −1 U(i1...iN −1),αN −1λ [N −1] αN −1(V †₎ αN −1,iN = X αN −2,αN −1 U(i1...iN −2),αN −2λ [N −2] αN −2(V †₎ αN −2,(iN −1αN −1)(V †₎ αN −1,iN = X α1,...,αN −1 Ui1,α1λ [1] α1(V †₎ α1,(i2α2). . . (V †₎ αN −1,iN = X α1,...,αN −1 Ui1,α1λ [1] α1B [2]i2 α1α2. . . B [N ]iN αN −1 (2.12)

where in the last step the V†_{matrices were reshaped to rank-3 tensors (except for the last matrix}

which was reshaped to a rank-2 tensor) as in the case of the left canonical MPS. Like the case of the left-canonical form, the right-canonical MPS of equation 2.12 can also be re-expressed in the notation of Ref. [19] by decomposing each B[l]il

αl−1,αl= Γ

[l]il

αl−1αlλ

[l] αl.

It is important to emphasize that although the matrices corresponding to each physical index in equation 2.12 are not the same as in equation 2.8, the states they represent are identical. The particular form of the MPS we obtain therefore strongly depends on the technicalities of the decomposition, in other words the virtual indices have introduced a gauge freedom. It is this property of the MPS that makes it possible for us to transform an MPS into a left- or right-canonical form.

There is another gauge that is particularly useful to determine expectation values of an MPS: the mixed-canonical form. This form can be obtained from a left- or right-canonical MPS by performing the SVD process in reverse until the site of interest is reached. In this form the tensors to the left of the site of interest form unitaries as in a left-canonical MPS and the tensors to the right of the site form unitaries as in a right-canonical MPS. Take, for example, an MPS consisting of four sites and say we want to determine the expectation value of a local operator at site two

hψ| O2|ψi = X α1,α2,i1,i2 M[2]i2 α1α2O i2i02 2 (M [2]i0 2 α1α2) †_. _(2.13)

as depicted in figure 2.6. This contraction reduces to the contraction of solely some local tensors. In the Vidal notation the conditions for the left- and right-canonical form can be reformulated to [5] X αl−1,il (A[l]il αl−1α0l )†A[l]il αl−1αl= X αl−1,il (Γ[l]il αl−1α0l )†λ[l−1]_α_l−1λ[l−1]_α_l−1Γ[l]il αl−1αl= δα0lαl (2.14) X αl,il B[l]il αl−1αl(B [l]il α0 l−1αl) †₌X αl,il Γ[l]il αl−1αlλ [l] αlλ [l] αl(Γ [l]il αl−1α0l )† = δα0 l−1αl−1 (2.15)

(17)

2.1. THE MATRIX PRODUCT STATE 17 A[l] (A[l]₎† = λ[l−1] λ[l−1]∗ Γ[l] (Γ[l]₎† = δ B[l] (B[l]₎† = Γ[l] (Γ[l]₎† λ[l] λ[l]∗ = δ

Figure 2.7: The left-canonical and the right-canonical conditions in the notation of Ref. [19].

. . . .

λB ΓA λA ΓB λB ΓA λA ΓB λB ΓA

Figure 2.8: An iMPS consisting of two sites repeated infinitely many times.

respectively, as graphically depicted in figure 2.7. When an MPS obeys both conditions 2.14 and 2.15 for all sites, the MPS is in a special mixed-canonical form which is also simply called the canonical form. Note that when the MPS is in a canonical form it also automatically forms a Schmidt decomposition on each bond. Therefore, a truncation performed on an MPS in this form will be optimal.

2.1.4 Infinite matrix product states

MPS can be used for systems of different sizes both with open and periodic boundary conditions1_.

Often, we want to make estimates of properties of systems extrapolated in the thermodynamic limit. One approach would be to consider different sizes of MPS such that we can extrapolate the results to an infinite chain. However, a more accurate path to take would be to probe the thermodynamic limit directly. In a translationally invariant system this can be achieved by taking one (or several) tensors and repeating them infinitely many times, see figure 2.8. An MPS of this form is called an infinite matrix product state (iMPS).

Local operators on iMPS can be evaluated in a similar manner as finite MPS when they are brought in the canonical form of equations 2.14 and 2.15. Take for example a two-site operator that we want to evaluate on an iMPS that is in the canonical form, as depicted in figure 2.9. The conditions 2.14 and 2.15 reduce this contraction to a simple form. Also, the SVD truncation discussed before will be optimal if the iMPS is in the canonical form, as in the finite case. Every iMPS can be brought into the canonical form by using e.g. the procedure explained in Ref. [20]. Note that it is not necessary for an iMPS to be in the canonical form when we want to compute expectation values, as we can also determine the left- and right extremal eigenvectors of the iMPS.

2.1.5 Area law of entanglement entropy

We have already seen that, when considering one bond of an MPS, the optimal approximation up to a certain bond dimension D is given by truncating the singular value spectrum. Some models can even be expressed exactly in terms of an MPS with finite D as is exemplified by the Affleck-Kennedy-Lieb-Tasaki Hamiltonian for which the solution can be represented by an MPS of D = 2 [5].

1_{Although MPS with periodic boundary conditions are more challenging, because it is not possible to bring}

(18)

. . . . λBΓAλAΓBλBΓAλAΓBλBΓAλAΓBλB λBΓ†_AλAΓ†_BλBΓ†_AλAΓ†_BλBΓ†_AλAΓ†_BλB O . . . . λBΓAλAΓBλB λBΓ†_AλAΓ†_BλB O

Figure 2.9: The expectation value of an operator on an iMPS can be reduced to a relatively simple contraction when it obeys the left- and right canonical conditions of equations 2.14 and 2.15.

To determine the accuracy with which a quantum state can be described using an MPS with a certain bond dimension one in principle has to know how quickly the singular value spectrum of the exact state decays. It turns out that the Schmidt decomposition from equation 2.10 has another useful property in this context. The singular value spectrum can be directly related to the entanglement entropy between the subsystems L and R. This can be shown by considering the reduced density matrix of, for example, subsystem L

ˆ ρL= TrR|ψi hψ| = X αn (λ[n]αn) 2_|u αniLhuαn|L (2.16)

inserting this into the formula for the von Neumann entropy shows that the entanglement entropy between subsystem L and R is

SL|R= − Tr ˆρLlog2ρˆL= − X αn (λ[n]_α n) 2_log 2(λ [n] αn) 2_. _(2.17)

Using this relationship, we see the entanglement entropy can be used to get a rough understanding of the singular value spectrum [5]. Therefore, for systems for which the states have a small entanglement entropy, it is likely that the MPS ansatz allows for an efficient simulation of the state.2

It turns out that many ground states and lowly-exited states have an interesting scaling behav-ior in their entropy. Whereas the entanglement entropy of random states scales with the volume of the subsystem under consideration [21], the entanglement entropy of many ground- and lowly ex-cited states scales with the boundary of the subsystem. This also implies that their entanglement entropy is much lower than that of random states. Even though a general prove for the existence of area laws is lacking they have been shown to hold for several important classes of models, such as for one-dimensional gapped Hamiltonians with a non-degenerate ground state as well as for several higher-dimensional gapped models, and they are expected to hold more generally for gapped and local Hamiltonians in any dimension [21]. It turns out that MPS also reproduce an area law in 1D and therefore they are especially well suited to describe such states.3 Area laws are known to break down in some cases or they can have a prefactor, as is the case for many critical models which are characterized by a logarithmic correction.

2.2 Projected entangled pair states

Because of the successes of MPS in one dimension, a lot of effort has been focused on extending them to two or even higher dimensions. An initial proposal was to use the same MPS description,

2_{The correspondence between a rapidly decaying singular value spectrum and a low entanglement does not hold}

necessarily, however in most practical cases it does [5].

3_{To be precise MPS satisfy S ≤ 2 log} 2D.

(19)

2.2. PROJECTED ENTANGLED PAIR STATES 19

(a) (b)

Figure 2.10: Proposals to generalize MPS to higher dimensions. (a) shows the snake MPS and (b) shows the projected entangled pair state (PEPS).

but instead of aligning the tensors in a one-dimensional chain they would be laid over a two-dimensional grid in a snake-like way (see figure 2.10a) [6]. One of the downsides of these algorithms is that all entanglement in the horizontal direction is carried over only one bond. But more importantly, these MPSs, unfortunately, do not obey the two-dimensional area law which scales as S ∝ L for the systems that we want to study and not as S ∝ const. as is deeply-rooted in the definition of MPS. Small widths, however, can be considered because they do not obey the general two-dimensional area law.

A state that does reproduce the area law in two dimensions is the projected entangled pair state (PEPS) as proposed in Ref. [6]. This state is obtained by adding additional legs in order to connect the tensors with its nearest-neighbors, as is shown in figure 2.10b. More technically, PEPS can be thought of in terms of maximally entangled pairs [6, 16]. In this picture the physical sites are thought of as consisting of pairs of two virtual particles, see figure 2.11a. Each virtual particle is maximally entangled (i.e. it has a bond) with its neighbor

|ωni = D

X

αn=1

|αnαni . (2.18)

The tensors at each site then represent projections of the pairs of two virtual particles back to a physical particle A[n]= d X in=1 D X αn−1,αn=1 A[n]in αn−1αn|ini hαn−1αn| (2.19)

such that the total wavefunction is given by

|ψi =A[1]⊗ A[2]_{⊗ . . . ⊗ A}[N ]_(|ω

1i ⊗ |ω2i ⊗ . . . ⊗ |ωN −1i) . (2.20)

When going to two dimensions, instead of two virtual particles at each site, four are introduced as depicted in figure 2.11b. The projectors then become

A[n]= d X in=1 D X αn−1,αn,βn−1,βn=1 A[n]in αn−1,αn,βn−1,βn|ini hαn−1αnβn−1βn| . (2.21)

This explains the name projected entangled pair states. At this point it is clear that there is nothing that prevents us from considering three or even higher dimensions and MPS turn out to be a special case of PEPS in one dimension (when using the term PEPS, we mean the two or

(20)

(a) _(b)

Figure 2.11: Virtual entangled pair state picture for (a) MPS and for (b) PEPS

higher dimensional variant from this point on). It is possible to represent other 2D geometries, such as the honeycomb or triangular lattices, using PEPS as well. However, for simplicity we will limit ourselves to the case of square and cubic geometries.

A practical difference between MPS and PEPS is that the latter cannot be contracted in an exact way. Also, a canonical form, which would allow for an exact contraction, does not exist for PEPS. In order to contract PEPS to determine expectation values we therefore must rely on approximate methods. Several of these methods will be described in section 2.4.

Similar to the case of MPS, states can be directly probed in the thermodynamic limit using PEPS by choosing several tensors and repeating these infinitely many times. These states are called infinite PEPS or iPEPS [22].

2.3 Ground state algorithms

In a previous section we saw that, because of the area law of entanglement entropy, MPS and PEPS are especially well suited to describe ground states. However, we have not yet touched upon the issue of how to obtain an MPS or PEPS that approximately represents the ground state of a given Hamiltonian. Note that MPS and PEPS are variational ans¨atze, i.e. the energy always forms an upper bound to the true ground state energy. The ground state can therefore be approximated by minimizing the energy

E = hψ| ˆH |ψi

hψ|ψi . (2.22)

In this section we will describe the two main methods through which an initial state can be brought to an approximate ground state. First, variational optimization will be briefly described, and we will discuss what the complexities with this method are in higher dimensions. After that we will discuss methods based on imaginary time evolution, which have been used in this work, in more detail.

2.3.1 Variational optimization

In its simplest form, variational optimization methods approximate the ground state by sweeping over the individual tensors and optimizing them along the way until convergence is reached. The optimization for a certain tensor Aiis achieved by taking equation 2.22 and solving

∂ ∂A†_i

h

(21)

2.3. GROUND STATE ALGORITHMS 21

which leads to a generalized eigenvalue problem

HiAi− µNiAi= 0 (2.24)

where µ is the Lagrange multiplier and Hi and Ni represent hψ| ˆH |ψi and hψ|ψi respectively

without the tensors at site i.

To evaluate Hi and Ni for MPS we can make use of the methods described before to efficiently

evaluate expectation values. When moving to two-dimensional PEPS we encounter the problem that expectation values cannot be evaluated exactly anymore. Approximate methods to determine expectation values exist for two-dimensional systems and they can be used to perform variational optimization on PEPS and even on iPEPS [4, 24].

Although these methods can produce very accurate results, their major downside is that they tend to become computationally expensive. This is caused by the requirement to contract the network often during the algorithm, especially for the determination of Hi which requires us to

perform a summation over every term in the Hamiltonian. For our purposes both the lack of an established algorithm to determine expectation values in three-dimensions and their computational complexity do not make this method a good starting point. Fortunately, there is another option: imaginary time evolution.

2.3.2 Imaginary time evolution

Another commonly used method to prepare a ground state is by imaginary time evolution. In this method an initial state |φi is projected onto the ground state |ψ0i by evolving it in imaginary

time. Let us write our initial state in the Hamiltonian basis |φi =P

nan|φni, such that ˆH |φni =

En|φni. Then perform an imaginary time evolution on this state

where E0 corresponds to the ground state energy and we have written ∆En = En− E0. So for

β → ∞ we indeed obtain the ground state. Imaginary time evolution works especially well for systems with an energy gap between the ground state and the first-exited state, as can be seen from this equation.

Let us rewrite e−β ˆH _{= (e}−τ ˆH₎M _{with β = τ M . To apply this method to our MPS the time}

evolution operator first has to be decomposed into local gates. This can be done by using the Trotter-Suzuki decomposition. Assume the Hamiltonian can be decomposed as ˆH = PN

i=1Hˆi,

where the ˆHi can, for example, describe nearest-neighbor interaction terms. In this case the

first-order Trotter-Suzuki approximation is given by [5]

e−τ ˆH = e−τ ˆH1_e−τ ˆH2_{. . . e}−τ ˆHN_{+ O(τ}2₎ _(2.26)

where the error is caused by the non-commutativity of the different Hamiltonian terms. τ is also called the Trotter parameter. Applying the terms, which from now on will be referred to as Trotter gates, in a different order will give a different Trotter error, however generally the computational cost will increase. In practice a useful ordering is given by the second-order Trotter-Suzuki decomposition

(22)

λB ΓA λA ΓB λB ΓA λA ΓB λB ΓA λA ΓB λB

Figure 2.12: By making use of the translational invariance of iMPS the Hamiltonian terms are applied in parallel in iTEBD. ΓA λA ΓB λB λB e−τ ˆHA,B SVD U λ0_A V λB(λB)−1U λ0A V(λB)−1λB λB Γ0A λ0A Γ 0 B λB

Figure 2.13: The steps taking place in the iTEBD algorithm. First the Trotter gate is applied, then the SVD is taken. Because of the application of the gate the bond dimension has grown, which therefore must be truncated. Finally, identities are inserted on the outer bonds such that the singular value spectra are reintroduced, and the updated tensors are obtained. Next, we shift to the other bond and the process is iterated.

which only gives a minor increase in computational cost and is the decomposition used in this work.

By assuming a translationally invariant system that is described by a Hamiltonian with local or nearest-neighbor interactions, the imaginary time evolution algorithm can be extended to probe systems in the thermodynamic limit. The assumption that the Hamiltonian terms are nearest-neighbor allows us to apply them in parallel, as depicted in figure 2.12. This algorithm is called infinite time-evolving block decimation (iTEBD) [19, 25, 26]. The computational cost of these algorithms is actually lower than is the case for finite systems, because a smaller number of sites has to be updated during each sweep.

The update step that is used in iTEBD is depicted in figure 2.13 and makes use of an SVD to obtain the new evolved tensors. After the application of the Trotter gate the bond dimension of the bond under consideration grows to D0 and is brought back to D by truncating the singular value spectrum.

A thing that could worry us is that the imaginary time-evolution operator is non-unitary. For the truncation to be optimal the iMPS should therefore first be brought back to the canonical form. However, the algorithm gives comparable results when this step is skipped. For small τ the Trotter gates are close to identities which means they bring the iMPS to an approximate canonical form. The error introduced by the non-unitarity of the operator turns out to be partially compensated by the SVD truncation and it vanishes when the Trotter step is taken smaller during the evolution [20].

Generalizing imaginary time evolution to higher dimensions

We encounter a problem when we want to generalize these algorithms to higher dimensions. When considering the singular value spectrum on a bond, it is no longer related to a Schmidt

(23)

decomposi-2.3. GROUND STATE ALGORITHMS 23 ΓA λA ΓB QA RAλARB QB SVD QA QB QA U λ 0 A V QB Γ 0 A λ 0 A Γ 0 B

Figure 2.14: The reduced tensor update for a 2D PEPS. The physical legs are split off from the rest of the tensors before the gate is applied. After applying the gate and performing an SVD the tensors with the physical legs are contracted giving the updated tensors.

tion as it is for MPS. Therefore, approximating a PEPS by truncating the singular value spectrum of the bond that has been updated no longer gives an optimal approximation. To make a proper truncation that minimizes

|ψi − | ˜ψi

a full contraction of the other tensors of the PEPS (also called the environment of the bond) should be performed. This approach is called full update (FU) [6, 22].

The FU method has a significant downside. Since the computational cost of calculating the environment is high, the efficiency of the algorithm decreases significantly. For this reason, com-putationally cheaper, although less accurate, approaches have been developed. The simplest of these is the simple update (SU), which gives a direct translation of the iTEBD algorithm and applies it to higher dimensions [27]. This means that, just like in iTEBD, only the singular values adjacent to the sites to which the Trotter gate has been applied are taken into account in the truncation. Because of the lower computational cost, the SU can reach higher D than the FU approach. However, this does not necessarily imply that the results from the SU will be more accurate. At a certain point the error that is caused by the inaccurate representation of the envi-ronment in SU will start to dominate and the results will converge to a value that deviates from the exact one. However, despite this rough approach, the algorithm does often produce reasonable results in systems that are characterized by a limited correlation length.

Several algorithms have been developed that lie in between SU and FU, giving an intermediate trade-off between accuracy and computational complexity. One interesting category is formed by methods performing updates on clusters, as e.g. described in Refs. [28] and [29].

Reduced tensor update

An approach that can be used to reduce the computational cost of the SU, which is especially useful to reach higher bond dimensions, is the reduced tensor update. The procedure is graphically displayed in figure 2.14. Before applying the gate, the physical legs are separated from the rest of the tensor by means of a QR-decomposition. This is followed by a contraction of the gate and an SVD as usual. Finally, the tensor containing the physical leg is contracted with the remaining tensor again and the process continues to another bond.

(24)

. .. _... _... _... . .. _... _... _... χ D D χD . .. _... _... _... χ

Figure 2.15: The MPS-MPO contraction as explained in the text. The outer edge of the PEPS is multiplied with the next layer and a truncation is performed to a truncation parameter χ, with which the accuracy of the contraction can be tuned. The process is repeated until only two layers are left, which can be contracted exactly.

2.4 Contracting tensor networks in two and three

dimen-sions

2.4.1 Two-dimensional contractions

In sections 2.1.3 and 2.1.4 it was shown that finite MPS and iMPS can be contracted efficiently and that the canonical form reduces the contraction to a simple form. One of the main difficulties in simulating higher-dimensional systems using PEPS is that contractions can only be performed for small system sizes in practice. A solution could be given by a canonical form for PEPS, however such a form has not been formulated so far and it is unknown whether such a form exists at all.

Several schemes have been developed to make approximate contractions of PEPS anyway, like the MPS-MPO contractions which were originally proposed in Ref. [6], corner-transfer matrix (CTM) (e.g. Ref. [30, 31]), tensor renormalization group methods (TRG) [32] and many variations of these. In the next paragraphs the general ideas behind these three methods will be described. It is important to note that even though PEPS is a variational ansatz the approximate contractions can still give an energy that is lower than the ground state energy.

MPS-MPO contractions

When looking back at the imaginary time evolution of section 2.3.2 and the scheme that was discussed there, one may realize that the network that is created by the time evolution gates actually forms a two-dimensional TN. One way to contract the TN created by PEPS is therefore to use a similar method [6]. Here we call these methods MPS-MPO4contractions. In this scheme one layer of the PEPS is multiplied with the next layer. After this step a truncation is performed to a bond dimension χ to counteract the growth of the bonds, e.g. based on a scheme similar to DMRG. By increasing χ the accuracy of the contraction can be increased. This process is repeated until only two layers are left, which can be contracted exactly. One step of the procedure is shown in figure 2.15.

The contraction is not limited to the vertical direction and also contractions in the horizontal or diagonal direction can be performed [23].

Corner-transfer matrix contractions

Corner-transfer matrix (CTM) methods form another way of contracting two-dimensional TN. In this method the environment around a certain site of interest is approximated by corner- and

4_{MPO stands for matrix product operator. It is possible to decompose a Hamiltonian, or any operator, acting}

on an MPS to a form in which one tensor per site is obtained. MPOs form very interesting objects in their own right and we refer the reader to Refs. [4] and [5] for a more detailed discussion.

(25)

2.4. CONTRACTING TENSOR NETWORKS IN TWO AND THREE DIMENSIONS 25 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. _... _... . .. . .. . .. _... _... . .. . .. CLU TU CRU TR CRD TD CLD TL χ _D

Figure 2.16: In the CTM method the environment around the site of interest is approximated by corner-and edge tensors, each representing a part of the TN. The bond dimension connecting the corner- corner-and edge tensors is often called χ. By varying χ the accuracy of the environment approximation can be tuned.

edge tensors which form approximative contractions of their respective regions of the TN, see figure 2.16.

Different variants of this algorithm exist. Here we will discuss the directional CTM algorithm for iPEPS as discussed in Ref. [30]. For simplicity we assume that the network only consists of one type of tensor. The approach will be too start from a simple 3×3 cluster and to grow the network in the horizontal and vertical directions until convergence of the environment tensors is achieved. A left-move is depicted in figure 2.17. We start by growing in the horizontal direction by introducing a new column of tensors. When the new column is contracted with the environment tensors on the left side the bond dimension in the vertical direction will grow. This growth can be counteracted by introducing projectors Z†_{Z ≈ I on both bonds, where Z is isometric. These projectors project} the bond down to a lower-dimensional subspace, effectively performing a truncation on the bond. It is important to choose the projectors in a proper way such that they act in the relevant subspace of the bond. Different ways exist to determine the proper projectors. In the literature the bond dimension that is reached after applying the projectors is often called χ. Varying χ allows us to tune the accuracy of the contraction. After this step the process continues in the other directions and is iterated until convergence is reached.

Just like for the MPS-MPO contraction, the CTM method can also be used as a way to contract the 2D TN of a 1D imaginary time evolution [31].

Tensor Renormalization Group

In the tensor renormalization group (TRG) a coarse-graining of the tensors is performed. TRG was originally proposed by Levin and Nave [32]. Since we are working with square and cubic lattices we will only explain the method for square lattices, but it can easily be generalized to other lattice types. The method starts by performing a decomposition of every tensor of the form

Tlurd=

X

i

SlidSuri (2.28)

which can be done by performing an SVD on every tensor and absorbing the square of the singular values in each of the S tensors. These steps are displayed in figure 2.18. This is followed by a coarse graining step in which four S tensors are combined into a new tensor reducing the total

(26)

CLU TU CRU TR CRD TD CLD TL CLU TU TU CRU TR CRD TD TD CLD TL TU CRU TR CRD TD I I TU CRU TR CRD TD Z Z† Z Z† ˜ CLU TU CRU TR CRD TD ˜ CLD ˜ TL

Figure 2.17: A left-move of the one-directional CTM method from Ref. [30]. First an expansion is done in one direction. Because of the growth in the bond dimension approximate identities are introduced in the form of projectors. After absorbing the projectors, the updated corners and edges are obtained. The absorption process is also done in the right-, up- and down-directions after which the full procedure is iterated until convergence is reached.

(27)

2.4. CONTRACTING TENSOR NETWORKS IN TWO AND THREE DIMENSIONS 27 T U V s U V √ s √ s l d r u i S S l r u d

Figure 2.18: In the first step of TRG a decomposition of all tensors is performed using an SVD. By absorbing the square roots of the singular values into the unitary matrices, the S tensors of equation 2.28 are obtained. To prevent the growth of the bond dimension the singular values are truncated.

. . . . . . . . . . . . . .. _... _... _... . . . . . . . . . . . . . .. _... _... _...

Figure 2.19: The TRG method for the cubic lattice for one iteration. In the first step an SVD is performed on all tensors (in the way depicted in figure 2.18) which gives the second lattice. The dashed lines indicate in what direction the SVD is performed. After this four of the resulting tensors are combined into one new tensor and the process is iterated until convergence is reached.

number of tensors by a factor of four. These steps are repeated until convergence is reached. One full step of the method is shown in figure 2.19. To prevent an exponential growth of the bond dimension a truncation is performed on the singular value spectra that are obtained in the first step.

The truncation in TRG, however, only involves local tensors and therefore it is not the optimal truncation we can perform globally. To make a globally optimal truncation the environment of the tensors should be taken into account. A method to do this is the second renormalization group (SRG) [33]. By using this method, the results of TRG can be dramatically improved.

Several other methods have been proposed to improve TRG and SRG, such as the higher-order tensor renormalization group (HOTRG) and the higher-higher-order second renormalization group (HOSRG) which make use of the higher-order singular value decomposition (HOSVD)5to further improve the truncations performed in TRG and SRG respectively [12]. Another method worth mentioning is tensor network renormalization (TNR) [34]. The TRG and the derived methods mentioned before do not give accurate results near critical points, because short-range entangle-ment is not completely removed during the coarse-graining steps. This problem is addressed in the TNR method, producing proper scale-invariance at the critical point and giving better results near criticality.

(28)

2.4.2 Three-dimensional contractions

The contractions mentioned above are the most commonly used contractions for two-dimensional PEPS. However, we encounter some problems when we want to generalize these schemes to three-dimensional networks. On the cubic lattice the number of legs on each tensor increases to seven, this implies that the amount of memory required to store the tensors, and more importantly the amount of memory and computation time required to perform manipulations on them, increases. This makes the methods mentioned above harder to perform, if at all. Also, to perform simulations of three-dimensional quantum systems the contraction becomes four-dimensional which is even more challenging.

Some adaptations to the contraction algorithms described above have been made and applied to 3D. For example, Refs. [11] and [12] have done this for TRG. In Ref. [11] a variant of TRG is proposed in which, instead of performing the truncation before the coarse-graining step, the truncation is performed afterwards. As an example, for a 3D contraction they provide results for the transverse field Ising model. In Ref. [12] the HOTRG is extended to 3D and it is shown to provide accurate results for the 3D classical Ising model.

Also, several alternative schemes have been proposed for 3D systems. In the string-bond states ansatz one-dimensional MPS-like strings are laid over a two- or three-dimensional lattice according to a certain pattern, e.g. along the rows and columns of the lattice, along its diagonals, in loops or a combination of these [13]. Expectation values are evaluated using MC sampling. Although the method produces larger errors than traditional PEPS algorithms because of the inability of string-bond states to describe correlations between sites that are not connected by strings properly, the decrease in computational cost makes it able to reach higher dimensional systems

Other schemes, for example, make use of weighted graph states [14], which do not depend on the geometry of the system and can therefore easily be extended to 3D, or the ab-initio op-timization principle (AOP) approach for TN [15, 35] in which a small few-body system is put in an entanglement bath which simulates the interactions between the few-body cluster and it environment.

2.4.3 Cluster contractions

In this thesis we want to look at an alternative way of approximating the contraction of a 3D tensor network. In our method, only a certain number of tensors (called a cluster) close to the sites that are being evaluated are contracted exactly. As we have seen in sections 2.1.3 and 2.1.4, an MPS in the canonical form can be contracted exactly by contracting singular value matrices at the edges of the sites under evaluation. In our method a similar contraction at the edges of the cluster is done to approximate the environment around the cluster in an effective way. In the simplest cluster that can be chosen, only the tensor(s) to which the operator is (are) applied are considered, as displayed in figure 2.20 for the cases of a one- and two-site operator in 2D (these clusters will later be called the 1 × 1 and 1 × 2 clusters respectively). Note that the environment approximation in the 1 × 2 cluster is similar to the effective environment used in the SU and uses the assumption that the PEPS can be contracted as if it were in a canonical form. The inaccuracy in this approximation is caused by the presence of loops in the PEPS. When the loops, however, are long as compared to the correlation length (i.e. l ξ with l the length of the loops) the contribution of the loops is expected to become small. Therefore, the evaluation is expected to be more accurate for PEPSs with small correlation lengths. Furthermore, by increasing the size of the cluster the contraction is expected to become more accurate as the loops inside the cluster are taken into account more accurately and the loops outside the cluster generally become longer, as shown in figure 2.21 for a 2D 3 × 4 cluster.

Other works making use of clusters in 2D exist, but they mainly focus on the update procedure. By increasing the size of the environment taken into account during the truncation procedure it

(29)

2.4. CONTRACTING TENSOR NETWORKS IN TWO AND THREE DIMENSIONS 29

(a) (b)

Figure 2.20: The simplest clusters for (a) one- and (b) two-operator evaluation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. _... _... _... _... _... _... _... . .. _... _... _... _... _... _... _...

Figure 2.21: Loops connecting the tensors at the boundary of the clusters generally become longer for larger clusters.

is possible to mediate between the SU and FU. In Ref. [28] an update method is proposed where imaginary time evolution is performed on rectangular clusters. In Ref. [29] a cluster update procedure on PEPS is described in which a row of the PEPS and a certain number of rows δ surrounding this row are approximated in detail. In this work we will solely focus on the effect of clusters on the evaluation of expectation values and not on their effect in the update procedure.

In chapters 4 and 5 the method will be benchmarked against the transverse field Ising model and the Heisenberg model in two- and three-dimensional quantum systems respectively. We will also examine the effect of taking different shapes of the cluster.

(30)

(31)

Chapter 3

Ising chain in a transverse field

This chapter discusses several results obtained by the iTEBD method in one dimension in order to show the capabilities of the methods and to get a better understanding of the limitations of the algorithm. The model that will be considered here is the Ising chain in a transverse field which is given by the Hamiltonian

ˆ H = −JX hi,ji ˆ σizσˆjz− h X i ˆ σx_i (3.1) with ˆσx

i and ˆσzi the x and z Pauli-spin matrices and with J and h two parameters that can be used

to tune the strength of the respective interactions. The model has a quantum phase transition1

from a ferromagnetic phase due to the first term to a paramagnetic phase due to the second term taking place at the critical point J = h. The transverse Ising model is one of the simplest models describing a quantum phase transition, therefore it is commonly used as a benchmark model. The model has several experimental realizations, such as in KH2PO4 for which it was originally

formulated [36].

One of the advantages of using the transverse field Ising model is that the Hamiltonian can be exactly diagonalized (see also Appendix A). The ground state energy per site is given by

E0= − 1 2π Z π −π dqpJ2_{+ 2J h cos(q) + h}2_, _(3.2)

which gives us the opportunity to compare the iTEBD algorithm with an exact result. From here on we will assume J = 1.

Simulations

Simulations have been run for several randomly chosen initial MPS for D = 2. To get a faster convergence these MPS have been used as initial states for the D = 4 case, which have been used again as initial state for D = 8 continuing to D = 55. Several points for which the accuracy was not satisfactory were run again from random initial conditions. During the evolution the Trotter parameter was gradually reduced from 0.01 to 0.01/28_.

Results and discussion

Figure 3.1 shows the relative error of the energy compared to equation 3.2 for different D. We see the MPS with small bond dimensions already give a good approximation of the ground state

1_{Quantum phase transitions are fully quantum in nature and take place upon the variation of a non-thermal}

parameter. In contrast to classical phase transitions, quantum phase transitions only take place at zero tempera-ture [8].

(32)

Figure 3.1: Relative error of the energy obtained for the transverse Ising model with iTEBD as compared to equation 3.1. The deviations near 10−14− 10−15

are mainly caused by the machine precision of the computer. The lines only form a guide to the eyes.

Figure 3.2: The order parameter for the transverse Ising model obtained from iTEBD around the critical point. The exact critical point is located at hc= 1.

energy away from the critical point. Near the critical point the MPS description starts to get less accurate. This can also be seen from the order parameter in figure 3.2. Particularly for smaller bond dimensions the critical point seems to be shifting away from hc = 1 and the transition

broadens. This behavior near critical points is common for MPS and it is also said that a pseudo-critical point is formed. For higher bond dimensions the pseudo-pseudo-critical point again approaches the exact critical point of the model. As shown in Ref. [37], it is possible to formulate scaling relations for the position of the pseudo-critical point depending on the bond dimension.

To get a better understanding of what is happening near the critical point we will look at the entropy and the correlation length. The entropy is obtained from equation 2.17 and is depicted for our simulation in figure 3.3a. Away from the critical point, where the states typically have a small entanglement entropy, it can be seen that the MPS of smaller bond dimensions and the ones with higher bond dimensions describe a similar amount of entropy. When we move closer to the critical point the entanglement entropy becomes stronger. However, the MPS, which has a finite bond dimension, can only describe a limited amount of entanglement which is especially apparent for the MPS with smaller bond dimensions.

A similar behavior can be observed when considering the correlation length as can be seen in figure 3.3b. The correlation length is obtained using a method described in Ref. [16]. For MPS correlations are known to decay exponentially, which means that a typical length scale for this decay can be determined.2 _{Assuming translational invariance of our MPS the effective correlation}

length is found to be [16] ξ−1= − ln λ2 λ1 (3.3)

2_{Non-critical gapped 1D ground states also have correlation functions which decay exponentially, giving another}

(33)

33

(a) (b)

Figure 3.3: The (a) entanglement entropy and (b) the correlation length around the critical point of the iTEBD simulation for the Ising chain with a transverse field.

where λ1and λ2 are the largest and second-to-largest eigenvalues of the so-called transfer matrix,

which is defined by T_α[n] 1α01α2α02 = d X in=1 A[n]in α1α2(A [n]in α0 1α 0 2 )†. (3.4)

Near the critical point we see that the effective correlation length grows strongly, indicating that long-ranged correlations become more important. At the critical point itself the correlation length will grow to infinity. The MPS description can only take a limited amount of this into account giving a reason for the decrease in accuracy in this region. Also, the broadening of the transition of the order parameter can be understood from this. For finite-size systems it is known that phase transitions spread out because of the inability of the correlation length to diverge [38]. The cut-off that is used for the bond dimension has a similar effect in our system, i.e. it induces a cut-off on the diverging correlation length. To study critical properties using MPS finite-size scaling analysis is commonly used.

(34)

(35)

Chapter 4

Results for two-dimensional

quantum systems

To determine whether the cluster contractions are able to give reasonable results and to see how the errors of the clusters behave for different cluster sizes and shapes, we have benchmarked them against the Ising model with a transverse field (from here on referred to as simply the Ising model) and the Heisenberg model.

The two-dimensional Ising model is the direct generalization of the one-dimensional model described in chapter 3. A useful characteristic of the Ising model is that by tuning the magnetic field parameter h the effective correlation length that is present in the ground state can be tuned, as we have already seen for the 1D case in figure 3.3b. In this way it is possible to get an idea of the accuracy of the cluster contraction for states with varying correlation lengths.

For the spin-1₂ Heisenberg model the Hamiltonian is given by ˆ

H = JX

hi,ji

ˆ

Si· ˆSj (4.1)

where J < 0 corresponds to the ferromagnetic model and J > 0 to the antiferromagnetic model. While the ground state of the ferromagnetic model corresponds to the alignment of the different spins, the ground state of the antiferromagnetic model turns out to be less trivial. In the classical Heisenberg model, the ground state corresponds to the N´eel state, in which the spins on neighboring sites are anti-parallel. Introducing quantum mechanics will cause the ground state to start showing deviations from the N´eel ordering. For both models no exact results are known in 2D, therefore the cluster contractions have been compared to several quantum Monte Carlo (QMC) and series expansion (SE) results. Also, the simulations have been directly compared to contractions with the CTM. For the CTM contraction a variant of the method described in section 2.4.1 is used in which the projectors are determined in the same way as in Ref. [39]. Appendix B gives the values that have been used for χ.

In the following sections several clusters that have been tried will be discussed.

Simulations

The ground state is approximated by the SU scheme for both models. Simulations have been run for several randomly chosen initial MPS and the ones that gave the lowest energy were selected. In contrast to the simulations that have been performed in 1D, the simulations here were done independently, i.e. without inserting the result obtained from simulations with lower D into ones for higher D. The Trotter parameter was gradually reduced from 0.01 to 0.01/28 during the time evolution.

(36)

Figure 4.1: The 3 × 3 cluster. The square roots of the singular value spectra have been absorbed into each of the tensors. For simplicity the cluster contractions will be represented by the graph on the right-hand side in the following, in which the gray tensor indicates the site to which the operator is applied.

Besides the error in the evaluation of the PEPS due to the cluster method, there are several sources of error in the simulations themselves. One of them is the error due to finite D, which formed the main source of error in the 1D case. Another error is due to the sub-optimality of the truncation scheme of the SU.

4.1 Square and rectangular clusters

The first type of clusters that will be described are the square and rectangular clusters. The square clusters of l × l tensors have been used to evaluate one-site operators and the rectangular clusters of l × (l + 1) tensors have been used to evaluate two-site operators. The 1 × 1 and 1 × 2 clusters have already been introduced in section 2.4.3. When going one step beyond the simplest clusters we obtain the 3 × 3 (see figure 4.1 where a simplified notation is introduced that will be used in the following sections) and 3 × 4 clusters. The largest clusters that have been implemented are the 7 × 8 cluster for the rectangular case and the 9 × 9 cluster for the square case. The leading contraction- and memory costs for our implementation of the clusters is given in appendix C.

First the error in the evaluation of the PEPS will be discussed. Figure 4.2 gives the relative error of the energy as compared to the CTM contraction for the Heisenberg model and for several field strengths of the Ising model. A significant decrease in the error can be observed for larger cluster sizes. The error is larger for states that encode more long-distance correlations, i.e. states with a field that is close to the critical point as well as states with a larger value for the bond dimension.

In figure 4.3 the relative errors for the order parameter and the staggered magnetization are given for the Ising and Heisenberg models respectively. Both quantities were evaluated with square clusters as well as with the rectangular clusters by inserting σi⊗ I with σithe Pauli spin matrices.

Similar trends are seen as for the energy, although in this case the error decreases monotonically with cluster size.

As mentioned before, there are other sources of error which enter the simulations due to the SU procedure and the effect of finite D. To get a better understanding of the influence of these errors in the overall method, the results have been compared to values from the literature for both models.

Figure 4.4 shows the energy per site and magnetization parallel and orthogonal to the magnetic field for the Ising model for D = 2. The evaluations have been done for the smallest and largest rectangular- and square clusters respectively. A rough estimate of the critical point is found by looking for the point where the order parameter vanishes. This point is found around hc ≈ 3.3.

Simulations performed for D = 3 − 4 show no significant shifts in this value. For comparison, table 4.1 gives several values for the critical point from the literature. The SE, cluster MC and fast-full update (FFU) results differ significantly from our result. When comparing our value to the SU result in Ref. [11] for D = 4, however, our result does agree reasonably. The error for

A simplified approach for 3D tensor network simulations

MSc Physics & Astronomy, Track theoretical

physics

Master Thesis

A simplified approach for 3D tensor

network simulations

Author:

Supervisor:

P.C.G. Vlaar

dr. P.R. Corboz

Student number:

Second examiner:

10218521

dr. J. van Wezel

Number of credits:

Examination date:

60 ECTS

August 25, 2018

Abstract

Populairwetenschappelijke samenvatting

Contents

Chapter 1

Introduction

Chapter 2

Introduction to tensor networks

2.1

The matrix product state

2.1.1

Diagrammatic notation

2.1.2

Obtaining the MPS representation

2.1.3

The canonical form

2.1.4

Infinite matrix product states

2.1.5

Area law of entanglement entropy

2.2

Projected entangled pair states

2.3

Ground state algorithms

2.3.1

Variational optimization

2.3.2

Imaginary time evolution

2.4

Contracting tensor networks in two and three

dimen-sions

2.4.1

Two-dimensional contractions

2.4.2

Three-dimensional contractions

2.4.3

Cluster contractions

Chapter 3

Ising chain in a transverse field

Chapter 4

Results for two-dimensional

quantum systems

4.1

Square and rectangular clusters