Thermodynamics of the Shastry-Sutherland model using tensor network methods

(1)

Thermodynamics of the

Shastry-Sutherland model

using tensor network methods

Kim William Torre

Master’s thesis

Institute for Theoretical Physics University of Amsterdam

Master of Science Theoretical Physics

Supervisor dr. Philippe Corboz Second reviewer dr. Edan Lerner

(2)

(3)

Abstract

The Shastry-Sutherland model is a frustrated two-dimensional spin 1/2 system that with the right tuning can effectively describe the quasi two-dimensional spin-gap ma-terial SrCu2(BO3)2. Large efforts have been dedicated in the last years to describe

its phase diagram, and to numerically study it. It is particularly challenging because accurate Monte Carlo simulations are lacking for this model due to the negative sign problem. Recently, progresses have been achieved through the use of tensor network methods, which are immune to the negative sign problem and are capable to perform largely unbiased simulations. In this thesis the latter model is simulated in order to obtain data at finite temperature using a novel approach based on different tech-niques from the field of tensor network methods. In doing so, benefits and limits of the various applied methods are studied.

(4)

Chapter 1 Introduction

1.1 Spin systems: Shastry-Sutherland model

One of the primary goals of theoretical physics is to provide models which are as simple as possible, but at the same time capable to capture various complex physical phenomena. In condensed matter physics, as well as in statistical mechanics, spin systems often serve in this role. Indeed, as many classical spin models have been a central tool to shed light on thermal phase transitions and critical phenomena, so are various quantum spin models now contributing to broaden our understanding of exotic quantum many-body states and quantum phase transitions. Moreover, even though they are not always intended to describe all the details of specific real-life experiments, spin systems have also proven, in some cases, capable to reproduce properties of real materials [1].

Indeed, a striking example of this, and main object of interest of this thesis, is the Shastry-Sutherland model, which is a frustrated spin system that for a particular tuning can effectively describe the two-dimensional spin-gap material SrCu2(BO3)2

[2] [3] [4]. Large efforts have been made over the past three decades to describe its phase diagram, but still many questions remain open today. Moreover, since the model is frustrated, accurate Quantum Monte Carlo simulations are not in general feasible due to the negative sign problem [5], and only recently has been possible to overcome this issue, largely thanks to novel numerical techniques that rely on the so called tensor network methods [6].

(6)

1.2 Tensor network methods

In numerical studies of quantum many-body systems, one of the most difficult chal-lenge is the exponential growth of the Hilbert space of quantum states. In general this feature prevents to simulate systems large enough to well approximate the ther-modynamic limit. Despite this, it is often the case that only a thin portion of the Hilbert space of quantum theories contains physical relevant states.

Indeed, tensor network methods have proven to be an invaluable set of techniques to address the curse of dimensionality, by exploiting the above mentioned fact [7] [8]. Tensors are mathematical objects that generalize the idea of multilinear maps, and a tensor network is simply a collection of these, connected by bonds which represent contractions between them. "Tensor network methods" is the term used to refer to the entire set of relative techniques, usually employed in quantum information science, condensed matter physics and computer science. These methods come with an intuitive diagrammatic language, that dates back to at least the early 1970s by Roger Penrose [9], and have seen many developments and extensions since their first appearance in the context of density matrix renormalization group (DMRG) [10]. The main reasons of their success is the fact that these methods allow to perform largely unbiased simulations, which moreover are also unaffected by the negative sign problem. Therefore, thanks to the latter, various classes of quantum systems can now be simulated more efficiently and studied in greater detail, and this has opened new avenues for a greater understanding of certain physical systems.

1.3 Outline

In quantum statistical physics the main object of interest is the partition function

Z = T r(ˆρ) = T r(e−β ˆH), (1.1)

where β is the inverse temperature and ˆH is the Hamiltonian of the underlying theory.

(7)

Therefore, thermodynamic properties can be studied if we are able to compute the density matrix operator ˆρ, or at least its trace, for system-size large enough to well approximate the thermodynamic limit.

This thesis focuses on the application of novel methods belonging to the above men-tioned class of techniques to study finite temperature properties of the Shastry-Sutherland model. For this purpose, first the density matrix operator of the system will be represented as a three-dimensional tensor network, and then a coarse-graining procedure will be performed on the latter network in order to obtain an approximate contraction of it, and therefore thermal expectation values such as energy, specific heat and magnetic susceptibility. The material is organized as follows.

In chapter two the model itself will be presented, illustrating its phase diagram and features. In chapter three it will be shown that, in the case of local Hamiltonians, is possible to build a tensor network that, up to some systematically controlled error, represents the density matrix operator [11] [12]. Then, a way to perform an approximate contraction of this large network will be presented, which is equivalent to take the trace of ˆρ [13] [14]. Results of the simulations, in which the above method was applied to the Shastry-Sutherland model and, as benchmarks, Heisenberg and transverse Ising model, are shown and discussed in chapters four, while in chapter five conclusions will be drawn regarding the applications, efficiency and limits of the used methods.

(8)

Chapter 2 Shastry-Sutherland model

2.1 Definition

Shastry-Sutherland model, also referred to as the orthogonal dimer model, is a frus-trated two dimensional spin 1/2 system named after its creators [2], and defined by the Hamiltonian ˆ H = J X <i,j> ˆ Si· ˆSj + JD X <<i,j>> ˆ Si· ˆSj, (2.1)

in which the first sum is performed over nearest neighbour sites while the second over next-nearest neighbours, as shown in Figure 2.1.

Figure 2.1: Interactions geometry of intra-dimer coupling (JD) and inter-dimer

coup-ling (J ) for the Shastry-Sutherland model

Spins are said to be frustrated if there exist competing interactions which cannot be satisfied simultaneously, and as a consequence of this the ground state of such systems often exhibits a degeneracy on the classical level and the magnetic order is

(9)

(a) (b) (c)

Figure 2.2: Three different regimes of the coupling constants and relative low energy classical configurations, where arrows indicate z-component of spin, blue bonds rep-resent negative energies and red ones positive. In (a) intra-dimer exchange prevails, in (b) the antiferromagnetic Néel phase emerges while in (c) the system is highly frustrated due to the competition between J and JD.

suppressed [15]. Indeed, in the Shastry-Sutherland model the frustration is caused by direct competition between intra-dimer coupling JD and inter-dimer coupling J ,

as we can see from the example in Figure 2.2.

Unexpectedly, after ten years from the work of Shastry and Sutherland, the model, for J/JD ≈ 0.63, has been experimentally realized in SrCu2(BO3)2, which has a

layered structure, with each Cu2+ _{ion carrying spin 1/2 [3], [4]. The material has}

exactly the same geometry up to only weak non-Heisenberg interactions, and it naturally lies in the so called dimer-product phase of the Shastry-Sutherland model.

2.2 Ground state phase diagram

Shastry-Sutherland model ground state presents three different phases, each living in a portion of the space of parameters J and JD. The most trivial is the already

mentioned dimer-product phase, for J/JD < 0.675, for which analytical results are

available. In order to illustrate these, consider the model for J = 0. With this choice is trivial to see that the Hamiltonian reduces to the one of decoupled dimers having a ground state made by a product of singlets Q _√1

(10)

exact ground state also for small values of J 6= 0.

For high values of J/JD, instead we find the well known antiferromagnetic Néel

or-der, and in between these two established regimes exists an intermediate phase that exhibits plaquette long-range order, as shown in Figure 2.3. The latter has been subject of many studies in the last years, in which other candidates like columnar-dimer order were proposed, and only recently [16] it has been possible to numerically confirm the plaquette nature of this intermediate phase.

Figure 2.3: Shastry-Sutherland model ground state phase diagram as a function of nearest-neighbor coupling J (JD = 1). The width of a bond is proportional to

the magnitude of the bond energy, and full (dashed) lines correspond to negative (positive) energies, Corboz and Mila [16]

(11)

Chapter 3 Method

3.1 Tensor network diagrammatic

As already mentioned, the methods applied are based on tensor network techniques, and their main purpose is to represent the density matrix operator as a tensor net-work, and subsequently compute an approximate contraction of it obtaining ther-modynamic observables as results.

In this section the diagrammatic notation used in tensor network methods will be quickly reviewed.

As shown in Figure 3.1, shapes, usually spherical, with legs are used to represent tensors. Thus, a vector is simply a shape with one leg, while a matrix would have two legs and so on.

(12)

A connected index between two tensors denotes a sum over the index, as depic-ted in Figure 3.2, or in other word, a contraction.

Figure 3.2: Tensor contraction using tensor network notation. Note that, usually triangular shapes represent left or right unitary matrices, used as isometries to rep-resent tensors in basis were they can be more efficiently truncated.

The great usefulness of this notation comes from the fact that we can easily handle contractions among tensors with many indices. Moreover, as quantum mechanics is formalized in terms of vectors and operators, the tensor network notation can be also used as an alternative way to work with quantum many body problems.

Indeed we can define a Matrix Product State [10], MPS, as the tensor network representation of a physical vector in the Hilbert space of a one-dimensional quantum theory with local Hamiltonian, and similarly a Matrix Product Operator, MPO. These concepts have been extended to two-dimensional systems, giving birth to Projected Entangled Pair States [17], PEPS. A graphical overview of the latter is given in Figure 3.3, in which vertical legs represent physical degrees of freedom and horizontal ones are virtual levels. The degeneracy of the latter levels is usually referred as "bond dimension".

(a) Matrix Product State (MPS) (b) Projected Entangled Pair State (PEPS)

Figure 3.3: At the left: wave function representation for a one-dimensional spin system with local Hamiltonian in terms of an MPS, where ik goes over local basis

states and Ψi1i2...iN is a rank-N tensor containing the expansion coefficients. At the

(13)

3.2 Singular value decomposition

All numerical methods used in this project rely on a very powerful set of techniques for treating matrices that are either singular or at least numerically very close to singular. This set of techniques, known as singular value decomposition, or SVD, is capable to diagnose and, in some cases, solve this problem, in the sense of giving an exploitable numerical answer. In this section it will be briefly shown an outline of the relevant theory regarding SVD.

This technique is based on the following linear algebra theorem [18] :

Any M × N matrix A can be written as the product of an M × N column-unitary matrix U , an N × N diagonal matrix s with positive or zero elements (singular values), and the conjugate transpose of a N × N unitary matrix V ,

A = U ·s · V† (3.1) Aab= N −1 X i,j=0 Uai sij Vbj∗ = N −1 X k=0 Uak sk Vbk∗ (3.2)

where in last equation we exploited the diagonal form of s. Therefore, for matrix V holds

V · V†=1 (3.3) N −1 X i=0 Vai Vbi∗ = δab. (3.4) When M ≥ N ,

(14)

U†· U =1 (3.5)

M −1

X

i=0

U_ia∗ Uib = δab. (3.6)

If instead M < N , then singular values sk are zeros for M ≤ k ≤ N − 1, applying the

convention to sort all of them into descending order, and as a consequence equation 3.6 holds only for a, b < M − 1.

The decomposition 3.1 is always doable, no matter how singular the matrix is, and it is unique up to an equivalent reordering of the columns of U , the elements of s, and the rows of V†.

Note that equation 3.2 can be thought of as a sum of outer products of columns of

U and rows of V†, with the “weighting factors” being the singular values sk.

In the case that most of the singular values are very small, then we can well-approximate A, as shown Figure 3.4, using only a few terms in the sum

Aab ≈ D−1 X k=0 ˜ Uak s˜k V˜bk∗ (3.7)

where ∼ indicates that to U , V or s have been truncated the last (N − D) columns or elements. Therefore, we need to store only a few columns of U and V , decreasing computational cost. Indeed, the latter is exactly one of the main pillar of tensor network methods, which prevail over exact computations due to the extremely lower computational complexity required.

Figure 3.4: Singular value decomposition and subsequent compression, using tensor network diagrammatic. ˜U , ˜V† and ˜s are the "D-truncated" version of U , V† and s, where D is the chosen cut to apply in order to represent faithfully A.

(15)

In order to deal with multi-rank tensors, instead of just matrices, there exists a generalization of the SVD, namely the higher order singular value decomposition, or HOSVD [19]. To illustrate how it works, consider a n-rank tensor T . We can decompose it as Ti1i2...in = X a1,a2...an Sa1a2...an U (1) i1,a1 U (2) i2,a2 ... U (n) in,an, (3.8)

with {U(j)} column-unitary matrices while Sa1a2...an the core tensor, which is the

analogous of a singular values matrix. This tensor possesses the property of all-orthogonality

hS:..ab..:|S:..ab0..:i = δab,ab0 ∀ab, (3.9)

where the bracket indicates an inner product between the two rank-(n−1) tensors, and it is pseudo-diagonal

|S:..ab..:| ≥ |S:..ab0..:| if ab < ab0, (3.10)

with |S:..ab..:| a tensor-norm defined as the square root of all elements’ square sum.

These norms have a similar role as the singular values of a matrix.

3.3 Building the network

3.3.1 Imaginary time evolution

The first step, in representing the density matrix operator ˆρ as a tensor network, is to identify the inverse temperature β with the length of an imaginary time interval and, subsequently, discretize it in N smaller intervals of length τ . In other words, we perform an imaginary time evolution of the density matrix operator itself, starting from its infinite temperature limit, which correspond to the identity operator, to some target temperature. Thus

ˆ

(16)

with

β = τ N, (3.12)

where H is the Hamiltonian of the underlying quantum theory.

In the language of tensor network this is straightforward, as we can see from Figure 3.5.

Figure 3.5: Imaginary time evolution depicted in terms of tensor network language. In general we can represent ˆρ as a dH× dH matrix, where dH is the dimension of

Hilbert space.

We can make the parameter τ arbitrarily small, by increasing the number of steps N in the imaginary time evolution, and treat a single density matrix at small inverse temperature. The advantage of this first factorization is that now we can apply various techniques to further decompose each ˆρ(τ ) into a product of local density matrices, as it will be shown in the next section. Therefore, the latter will be our starting point for the next step.

(17)

3.3.2 Trotter-Suzuki decomposition

The first way with which we can construct the tensor network, and therefore represent the density matrix of the system, is through the use of the renowned Trotter-Suzuki decomposition [11].

For bounded operators A and B, holds

eA+ ˆˆ B = lim n→∞ eAnˆ e ˆ B n n (3.13)

The latter can be generalized for a set of non commuting operators { ˆAi}, such that

e ΣiAˆi _{= lim} n→∞ Πie ˆ Ai n n . (3.14)

We can apply this formula for systems with local Hamiltonian ˆH, and so such that

ˆ

H =X

i

ˆ

Hi, (3.15)

where { ˆHi} is the set of local Hamiltonians.

If now we identify n with the inverse of the negative imaginary time step τ , we have

e −τ ΣiHˆi _{= Π}

ie−τ ˆHi + o(τ ). (3.16)

Performing a singular value decomposition on each local density matrix ˆρl, see Figure

(18)

Figure 3.6: Tensor network representation of an SVD for the case of a local density matrix at inverse temperature τ .

The latter can be translated into tensor network language, as shown in Figure 3.7, in which a one-dimensional system is considered. Note that, in general this de-composition does not preserve internal symmetry of ˆH, and as a consequence the representation of ˆρ will not be hermitian. In some cases, by means of SVD, we can manually restore this symmetry in the network, indeed in Figure 3.7 and 3.8 we can distinguish two different schemes to create it.

Figure 3.7: (Scheme A) Trotter-Suzuki decomposition and subsequent singular value decomposition of the density matrix for a one-dimensional system, with local Hamiltonian, consisting of four sites. Since each T tensor has not equivalent bonds in the upper and lower part of the network, this scheme does not preserve the symmetry of ˆρ.

In the first case, scheme A, ˆρ is represented as a not hermitian object, while in the second case, scheme B, the symmetry is preserved. Note that, these two ways of decomposing the density matrix correspond to a first or second order Trotter-Suzuki decomposition respectively. Also, in Figure 3.9 is shown an example of Trotter-Suzuki decomposition for a two-dimensional local Hamiltonian, in which the her-mitian symmetry is broken.

(19)

Figure 3.8: (Scheme B). T-S decomposition and SVD of ˆρ for a one-dimensional system, with local Hamiltonian, consisting of four sites. The use of this scheme preserves the symmetry of the density matrix, because now T has equivalent bonds in both upper and lower part of the network.

Figure 3.9: An example of scheme to decompose a density matrix for a two-dimensional system made by four sites. Note that the hermitian symmetry of ˆρ is lost because of the non commuting different layers that form the network repres-entation.

(20)

Therefore, at first sight, it seems that Trotter-Suzuki decomposition might introduce some sort of systematic error due to the breaking of the hermitian symmetry, other than the expected error of order τ . Nevertheless, this kind of factorization is widely used giving, in general, accurate results.

At this point, using 3.11 and 3.16 we obtain,

ˆ

ρ(β) ≈Πie−τ ˆHi

N

, (3.17)

and so we can effectively represent ˆρ as a tensors network, up to an error of order τ , see Figure 3.10.

Figure 3.10: Tensor network representation of the density matrix operator for a one-dimensional system, with local Hamiltonian, where the horizontal dimension is physical space and the vertical is imaginary time. Note that in this framework is easy to identify the quantum-classical correspondence, indeed a one-dimensional quantum system is linked to a two-dimensional network, which can be thought as the corresponding classical statistical system.

3.3.3 Symmetric cluster expansion

As an alternative to Trotter-Suzuki decomposition, it is possible to use a recently developed method [12], also based on tensor network techniques, namely a symmetric cluster expansion. The advantages of it is that it preserves all internal symmetries of the underlying quantum theory while, at the same time, it is size extensive, and

(21)

therefore it scales well in the thermodynamic limit.

This symmetric cluster expansion of ˆρ can be seen as a sum of all possible clusters up to a fixed size, where for each cluster the expansion is exact up to infinite order. In order to illustrate the method, consider a spin chain with only nearest-neighbour interactions, such that the Hamiltonian of the system can be written as

ˆ

H =X

i

ˆ

hi,i+1, e−τ (ˆh1,2+ˆh2,3) (3.18)

with ˆhi,i+1 local Hamiltonian.

The idea is to construct an MPO, see Figure 3.11, which captures a systematically chosen number of terms from the Taylor expansion of e−τ ˆH_{, and in such a way that}

the perturbative expansion results exact in each cluster if we ignore terms outside of it. Thus, adding a certain component to the MPO will encode the action of a local density matrix in the respective cluster.

Figure 3.11: MPO representation of the density matrix for a one-dimensional spin chain with only nearest-neighbour interactions. Vertical indices are physical and horizontal are virtual.

The smallest cluster has size p = 1, and is encoded adding the first component to the MPO, which is just the identity operator acting on a single site, as shown in Figure 3.12. Note that, in this diagram, vertical legs are physical indices, while horizontal are virtual indices, and the label ”0” denotes their first entry.

(22)

For p = 2 we have the first non-trivial cluster, and to encode it we subtract the already incorporated unit to the local density matrix ˆρi,i+1. In order to add the

respective components to the MPO we introduce a new virtual level with label ”1” and degeneracy depending on the SVD of ˆρi,i+1, as shown in Figure 3.13. So, adding

O₀₁0 and O0₁₀ to the MPO O + O0 will encode the effect of a local density matrix acting on a two-site cluster.

Figure 3.13: Cluster of size p = 2.

Now, to obtain the component which encodes the action on a cluster of size p = 3, we consider a local density matrix acting on three sites instead, and subtract what was already contained in the two-cluster terms. Then we apply, to the right hand side of Figure 3.14, the inverse of O01 and O10, obtaining O110 without increasing

the bond dimension of the MPO, as shown in Figure 3.15. This procedure can be continued adding more and more components, and therefore, clusters of bigger sizes.

Figure 3.14: Cluster of size p = 3. The unlabeled indices denote a sum over all virtual levels.

The MPO constructed from this method preserves the symmetry of ˆρ and is correct up to order τp−1_{, with p the size of the largest cluster incorporated. Moreover, the}

expansion captures a large portion of the term of order τp too. Indeed, at this order, the terms that are lacking from the MPO are the ones coming from the cluster of size p + 1. To see how the error scales in this case, consider ˆρ(τ ) for a system with

(23)

Figure 3.15: Contraction of a three-sites cluster with the inverse of O01, to the left,

and O10, to the right, which return the new component O110 .

distinct terms and, among these, only p! are of the kind that the MPO misses. Thus, there is only a fraction p!/pp of terms at order τp which are missing. And so, for large p the error scales as

 ∼√pe−p.

The method can also be extended to higher dimensions and, as an example, will be briefly illustrated for a two-dimensional system with nearest-neighbour interactions in a square lattice.

Again, the goal is to construct a PEPO, which is the two-dimensional generalization of an MPO, that represent the density matrix of the system, see Figure 3.16. Note that, physical legs of the PEPO are now assumed to be normal to the sheet-plane, and therefore, they are not drawn.

Figure 3.16: PEPO representation of the density matrix for a two-dimensional spin system with only nearest-neighbour interactions in a square lattice. Note that phys-ical indices are normal to the sheet-plane and therefore they are not drawn for simplicity.

(24)

possible to encode some four- and five-sites clusters without the need to introduce a new virtual level, see Figure 3.18.

Figure 3.17: Clusters of one,two and three sites for a square lattice local Hamiltonian.

Figure 3.18: Clusters of four and five sites for a square lattice local Hamiltonian.

This technique can also be generalized for different kinds of Hamiltonians or lattice geometries, including the already introduced Shastry-Shuterland model. For the sake of clarity, the setup used for the Shastry-Sutherland model is shown in Figure 3.19. From this picture is possible to see that in this case there is not a translational invariant network, since A and B are not equivalent, and therefore we need to make

(25)

two expansions, one for each fundamental tensor. Indeed, in this setup the expansion was carried out using a dimer representation of the PEPO, in order to facilitate computation when the system was simulated in the dimer-product phase regime of parameters J and JD. Clearly this choice introduce a bias for the simulation that

must be taken into account during the computation.

Figure 3.19: Simulation setup used for the Shastry-Sutherland model in which next-nearest neighbour sites, corresponding to diagonal interactions, are blocked together. The physical open indices of each tensor are normal to the sheet-plane and are omitted for simplicity.

To conclude, the accuracy of the cluster expansion can be systematically controlled through the number of components encoded in the MPO, or PEPO, and thus the bond dimension, or degeneracy of the virtual levels. This, together with the fact that internal symmetries are preserved, is the advantage gained from the technique, that Trotter-Suzuki decomposition cannot provide.

3.4 Contracting the network

Once that the density matrix operator, by means of Trotter-Suzuki decomposition or symmetric cluster expansion, has been represented as a tensor network, we have to deal with the problem of contracting it, as shown in Figure 3.20, in order to obtain thermal expectation values of some local operator ˆO,

(26)

or simply the free energy

F = −1

βln Z. (3.20)

An exact contraction of the network is possible only for small system sizes, and so, in order to not exceed machine memory, we need a method to coarse-grain the network in such a way that computational cost can be kept low.

Figure 3.20: Example of network exact contraction, with an applied local two-sites operator ˆO, for a one-dimensional three-site system with periodic boundary con-dition. This contraction approximate the trace of the product (eβ ˆH · ˆO), up to a Trotter-Suzuki or symmetric cluster expansion error.

3.4.1 Tensor renormalization group method

To achieve the latter goal, a tensor renormalization group, TRG, method was first proposed by Levin and Nave in 2007 [20] for two-dimensional networks, based on the SVD. Later, in 2012, a new coarse-graining TRG scheme, based on the HOSVD, has been presented [13] to treat two- or three-dimensional networks. The advantage of this last method is that HOSVD takes into account more accurately the interplay

(27)

between different components of a tensor and thus, provides a better scheme to truncate a local tensor than the SVD.

To illustrate it, let’s first consider a one-dimensional, translational-invariant, quantum system. In this case we have to deal with a two-dimensional network, as already shown, where the horizontal dimension represents physical space, while vertical rep-resents the imaginary time direction, as shown in Figure 3.10. Contracting the full network we get as usual the partition function

Z = T r(ˆρ) =Y

i

Txix0iyiyi0. (3.21)

To coarse-grain in vertical direction we sequentially contract two tensors into a single one, reducing the vertical lattice size by a factor of two, as shown in Figure 3.21.

Figure 3.21: Contraction of a a two-dimensional tensor network along the vertical axis on a square lattice.

Therefore, at each step two local tensors are contracted defining

M_xx(n)0_yy0 = X i T_x(n) 1x0₁yiT (n) x2x0₂y0i, (3.22) with x ≡ x1⊗ x2, x0 ≡ x01⊗ x 0

2 and n denoting the nth iteration of the coarse-graining

process.

The dimension of the horizontal legs of M(n)_{, so called bond dimension, will be the}

(28)

M_xx(n)0_yy0 =

X

i,j,k,l

SijklUxiLUxR0_jU_ykUU_yD0_l, (3.23)

where the U ’s are unitary matrices and S is the core tensor. Note that, the same relationship is depicted in Figure 3.22.

Figure 3.22: Relation 3.23 depicted using tensor network diagrams. Note the great usefulness of this notation, capable to largely simplify standard notation of tensor contractions.

At this point, we need to renormalize the horizontal legs x and x0, but note that, due to the fact that we have a translational invariant network, the right leg of M(n)

is directly linked with the left one of an identical tensor in the neighbouring right site. Therefore, we only need to renormalize one of the two horizontal bonds to automatically truncate the opposite bond too. In order to choose which leg to treat, we compare the quantities

1 = X i>D |Si:::|2 2 = X j>D |S:j::|2, (3.24)

with D the chosen truncated bond dimension. If 1 < 2, the first i-th leg of S, or

UL_{, is truncated to D. Otherwise, the j-th leg of S, or U}R_{, is renormalized.}

To find, for instance, UL we first reshape M_xx(n)0_yy0 into a matrix M_x,x0 0_yy0,

(29)

with x as the row index and (x0, y, y0) as column index. If we perform an eigenvalue decomposition of M0M0†we can directly obtain UL, thanks to the fact that it is equal to the left unitary matrix of M0 under a singular value decomposition, see Figure 3.23, indeed

M0M0† = ULΛL(UL)†, (3.26)

with ΛL _{diagonal matrix. It can also be shown that}

|Si:::|2 = ΛLi . (3.27)

Figure 3.23: Cheapest procedure in terms of computational complexity to obtain, for instance, the left unitary matrix.

Thus, after the truncation of the U ’s, we can actually update M(n) _{into T}(n+1)_{, as}

shown in Figure 3.24, by contracting it with these isometries, in such a way that

T_xx(n+1)0_yy0 = X i,j (U_ixL)∗ M_ijyy(n)0 U_jxL0 ₁ < ₂ (3.28) T_xx(n+1)0_yy0 = X i,j (U_ixR)∗ M_ijyy(n)0 U_jxR0 ₂ < ₁ (3.29)

After this last update, a vertical coarse-graining step is concluded and now each ho-rizontal layer of the coarse-grained network represents a density matrix at a doubled inverse temperature. This procedure can be iterated and, using a finite number of vertical steps N , we can target some selected inverse temperature β = τ N . Moreover, for large values of N we can obtain ground state expectation values.

(30)

Figure 3.24: Renormalization step of the vertical HOTRG procedure for a two-dimensional network. Applying the truncated isometries U and U† to a vertically contracted pair of T(n) _{we obtain T}(n+1)_.

In any case, once the network is completely coarse-grained in the vertical direction, we are left with a one-dimensional tensor network, representing ˆρ(β), that can be exactly contracted completing the operation of trace of ˆρ, or can be first contracted with a local operators obtaining thermal expectation values from the trace.

Note that, due to the truncation done to the network during the HOTRG, there is no guarantee that it will represent an hermitian object after the coarse-graining procedure. Therefore, a variant of this algorithm, which was used in this project, is to interrupt the vertical HOTRG one step before the chosen target inverse temperature β and to perform exactly a contraction between the upper-horizontal layer of the network and its mirrored-conjugate image, as shown in Figure 3.25. In this way the final matrix obtained is hermitian by construction.

Figure 3.25: Final vertical step in which two layers, each representing e−β2Hˆ, are

contracted in such a way that the trace is performed on a matrix which is hermitian by construction.

(31)

The method is easily extendable to two-dimensional quantum systems, and therefore three-dimensional tensor networks, see Figure 3.26. For completeness, in Figure 3.27 is depicted the truncating scheme of the vertical coarse-graining procedure.

Figure 3.26: Contraction of the a three-dimensional tensor network along the vertical axis on a square lattice.

Figure 3.27: Renormalization step of the vertical HOTRG procedure for a three-dimensional network. In this case we need to apply two pairs of isometries V (y-axis) and U (x-axis).

Note that, the illustrated method provides only few data points, due to the fact that at each coarse-graining step the network vertical size is divided by a factor of two, and therefore imaginary time is doubled. Thus, an alternative way of coarse-graining the network is to interrupt the above explained scheme at some low inverse temperature β and continue to coarse-grain absorbing layer by layer, as depicted in

(32)

Figure 3.28: In order to obtain data points for more temperatures we can interrupt the standard coarse-graining procedure after k steps, at some small inverse temper-ature β = 2e kτ , and continue to contract vertically layer by layer. In this way we

would obtain thermodynamic observables at inverse temperature τ , 2τ , ..., β, 2e β,e

3β,..., β.e

Till now, the coarse-graining procedure has been illustrated only for the vertical dimension of the network, but in principle it can be used also for the horizontal one. However, in the case of quantum systems this direction represents physical space and, as it will be shown in next chapter, at high temperature the singular values of M decay extremely slow, making the coarse-graining procedure really expensive in that direction.

Moreover, our initial assumption of translational invariance is not fulfilled in general and does not apply to the Shastry-Sutherland model, therefore the latter scheme must be slightly modified in order to be applied. The generalization of it, to non-translational invariant networks, is partially inspired to the already tested generalized corner transfer matrix method, which will be treated in the next section.

3.4.2 Corner transfer matrix method

In the last section a method to perform an approximate vertical contraction of the network has been presented, but when we deal with quantum systems the latter scheme cannot be also applied efficiently to contract in the horizontal dimension. For one-dimensional quantum systems it is possible to ignore this problem because after the vertical coarse-graining we are left with a one-dimensional network, which can be exactly contracted. However, this is not the case for two-dimensional quantum systems, for which a two-dimensional network is left after the vertical HOTRG, or at least it is only for very small system sizes.

(33)

the horizontal directions, it is possible to employ the so called corner transfer matrix method, CTM [14]. The idea behind this technique is to recursively build an en-vironment of tensors which systematically approximates an infinitely large network around a bulk-tensor, see Figure 3.29.

Figure 3.29: Corner transfer matrix method, where {Ci} and {Ti} are

environment-tensors, with bond dimension χ, and a is the bulk-tensor of an infinite network. Picture taken from P. Corboz (Dresden 2018)

To illustrate the method in details, consider for simplicity an infinite two-dimensional rotational- and translational-invariant network which has been obtained by perform-ing a vertical HOTRG to a three-dimensional network. In this case we can build the environment using only a symmetric corner-matrix C together with an edge-tensor T , and to begin we initialize them by simply projecting the needed indices of the bulk-tensor a into, for instance, their zero-th components, see Figure 3.30. Then, we let the system grow in all directions and contract C and T with bulk-tensors as showed in Figure 3.31. Note that D denote the truncated bond dimension used during the vertical HOTRG, while χ indicate the bond dimension of C and T .

(a) (b)

(34)

Figure 3.31: System grows in all direction (rotational- and translational-invariant network). Picture taken from P. Corboz (Dresden 2018).

From a subsequent eigenvalues, or SVD, decomposition of corner-matrix C we can obtain the isometries needed to project the environment-tensors into a basis where it is possible to, as usual, efficiently truncate their bond dimension from χD to a χ, see Figure 3.32 and 3.33, by discarding smaller singular values. A complete CTM step is illustrated in Figure 3.34.

Figure 3.32: Eigenvalue decomposition of a corner, with U unitary matrix. Picture taken from P. Corboz (Dresden 2018).

Figure 3.33: Updating step for C and T , where ˜U is the version of U with truncated columns, which reduces the bond dimension to χ. Picture taken from P. Corboz (Dresden 2018).

(35)

Figure 3.34: CTM step for a rotational- and translational-invariant network. Picture taken from P. Corboz (Dresden 2018).

This procedure can be iterated until we reach convergence, which is dictated by the condition

X

k>χ

|sk− s0k| < tol, (3.30)

where {sk} and {s0k} are the normalized singular values, from two successive steps,

which have been truncated, while "tol" is a chosen small tolerance value.

Once the method reached convergence we can use these environment-tensors to com-pute thermal expectation values of local operators acting on a single network site, as shown in Figure 3.35, and similarly for operators acting on more sites.

(a) (b)

Figure 3.35: Thermal expectation value of ˆO for a 2D-quantum system (a). Tensor a is obtained by tracing all physical indices of the 3D-network fundamental tensor TH,

which has already been vertically coarse-grained through HOTRG. Instead tensor

b is obtained by contracting TH with a one-site operator ˆO. Picture taken from P.

(36)

networks [21] [22], which is the case of Shastry-Sutherland model. The procedure is basically unchanged apart from the fact that now we need to build an environment which takes into account the different features of the tensor network, as depicted in Figure 3.36.

Figure 3.36: CTM environment of Shastry-Sutherland model in dimer-product rep-resentation. Each corner-matrix, or edge-tensor, is distinct from the others since the network is not isotropic in this case.

(37)

Chapter 4 Results and Discussion

In this chapter the results from various simulations will be shown and discussed, ranging from a one-dimensional Heisenberg chain to the already presented Shastry-Sutherland model. In the first part the effect of the broken hermitian symmetry of ˆρ will also be shown in detail, in the case of networks constructed from a Trotter-Suzuki decomposition instead of a symmetric cluster expansion, together with the reason why coarse-graining in the horizontal direction of the network is really expensive in terms of computational cost and memory.

4.1 One dimensional Heisenberg model

In order to test the method with less challenging models than the Shastry-Sutherland, an antiferromagnetic Heisenberg chain with periodic boundary condition was first simulated, with the usual Hamiltonian

H = J L X i ˆ Si· ˆSi+1, (4.1)

with ¯Si the spin 1/2 operator and L the number of sites. The system was simulated

starting from either a Trotter-Suzuki decomposition or a symmetric cluster expansion using a value of J = 1 and, in the T-S case, the two different network construction schemes were tested.

(38)

the fact that the last vertical step of the network contraction is performed exactly and in such a way that the final object is symmetric under transposition and conjugation. Indeed, if the network is made by tensors from scheme "A", and so not symmetric, an error is introduced in this final contraction, since Txx0_yy0 is not equivalent to T_xx0_y0_y,

as already shown in Figure 3.7 and 3.8.

0.0 0.1 0.2 0.3 0.4 τ 0.000 0.005 0.010 0.015 Ab so lut e r ela tiv e e rro r

Free energy per site (absolute) relative errors

Vertical bond dimension 8 (scheme A) Vertical bond dimension 8 (scheme B)

(a) 0.0 0.1 0.2 0.3 0.4 τ 0.000 0.005 0.010 0.015 0.020 Ab so lut e r ela tiv e e rro r

Energy per site (absolute) relative errors

Vertical bond dimension 8 (scheme A) Vertical bond dimension 8 (scheme B)

(b)

Figure 4.1: Comparison between scheme "A" and scheme "B" (see 3.3.2) Trotter-Suzuki decomposition. Free energy (a) and energy (b) relative error with respect to the exact values, for a 4 sites (L) chain with periodic boundary condition, and network vertical size (N) of 4. The horizontal HOSVD bond dimension used is 8.

Also, a notable matter is the coarse-graining process in both vertical, or imaginary time, and horizontal directions. It has been already mentioned that the proced-ure was carried out only for the vertical dimension of the network, while an exact contraction was performed after the end of the imaginary time evolution. We can directly see from Figure 4.2 that singular values decay very differently in the two directions.

(39)

1 4 8 12 16 0.00 0.25 0.50 0.75 1.00 Si ng ula r v alu es

Horizontal singular values spectrum

τ=0.001 τ=0.01 τ=0.1 τ=1 (a) 1 4 8 12 16 0.00 0.25 0.50 0.75 1.00 Si ng ula r v alu es

Vertical singular values spectrum

τ=0.001 τ=0.01 τ=0.1 τ=1

(b)

Figure 4.2: Vertical (a) and horizontal (b) singular values, regarding the first re-spectively horizontal and vertical coarse-graining step for a 4 sites chain. In the first case, for τ small enough, the spectrum decays quickly and the bond dimension can be kept low, making the truncation scheme efficient. In the second case instead, for any practical imaginary time steps, singular values do not decay fast enough, behav-ing in the opposite way. The network was built from a second order Trotter-Suzuki decomposition (scheme "B")

To understand this, consider that in the HOSVD horizontal spectrum 3.27 we can link singular value magnitudes to entanglement between two neighbouring sites or blocks, while in the case of vertical spectrum we are actually compressing physical rather than virtual levels. As a consequence, and considering that for high temperature

ˆ

ρ = e−τ ˆH ∼ 1, (4.2)

which gives singular values of same magnitude, discarding any of them, and so trun-cating any states, would inevitable lead to a concrete loss of information. Therefore, physical space compression results to be too much expensive in terms of compu-tational cost, forcing us to keep a great portion of singular values at each coarse-graining step. Instead, as we can see in Figure 4.3, the procedure works perfectly for imaginary time direction, where only a few singular values need to be kept in order

(40)

0.2

0.4

0.6

0.8

1.0 Inverse truncated bond dimension 1/

D

0.0

0.2

0.4

0.6 En

erg

y r

ela

tiv

e e

rro

r

Energy per site relative error for different D

τ

=0.001

τ

=0.01

τ

=0.1

Figure 4.3: Relative error of energy density as function of truncated horizontal HOSVD bond dimensions D, for a 4 sites (L) chain with periodic boundary condition, and network vertical size (N) of 4. A second order Trotter-Suzuki decomposition was employed to construct the network.

Moreover, a comparison between Trotter-Suzuki decomposition and symmetric cluster expansion was performed. In order to exclude errors from the coarse graining trunca-tion scheme, an exact contractrunca-tion of a two-by-four network built from a T-S decom-position, and therefore a four sites chain at inverse temperature 4τ , was performed together with an exact contraction of a four-by-four network constructed from a symmetric cluster expansion with p = 3. In Figure 4.4 we can observe the usefulness of the latter method, which, in addition to providing a smaller error, preserves all internal symmetries of ˆρ without the use of artificial schemes as in the T-S case.

(41)

0.0 0.1 0.2 0.3 0.4 τ 0.000 0.001 0.002 0.003 0.004

Ab

so

lut

e r

ela

tiv

e e

rro

r

Free energy per site relative errors

Trotter-Suzuki

Symmetric cluster expansion

(a) 0.0 0.1 0.2 0.3 0.4 τ 0.000 0.002 0.004 0.006 0.008

Ab

so

lut

e r

ela

tiv

e e

rro

r

Energy per site relative errors

Trotter-Suzuki

Symmetric cluster expansion

(b)

Figure 4.4: Free energy (a) and energy (b) relative error with respect to the exact values, for a 4 sites (L) chain with periodic boundary condition, and network vertical

(42)

To conclude with the Heisenberg chain, a large system, with an external magnetic field h applied, was simulated and its results are shown in Figure 4.5, 4.6 and 4.7, where data from a different method, namely transfer-matrix renormalization-group method [23], are used as reference. In this case more simulations were run for different τ ’s, in order to obtain more data points.

0.0

0.5

1.0

1.5

2.0 Temperature

−1.5

−1.0

−0.5

Fr

ee

en

erg

p

er

sit

e

Free energ per site for different external fields

SCE h=0 T-S h = 0 SCE h = 1 T-S h = 1 SCE h = 2 T-S h = 2 SCE h = 2.5 T-S h = 2.5 TMRG h = 0 TMRG h = 1 TMRG h = 2 TMRG h = 2.5

Figure 4.5: Heisenberg spin-1/2 chain. Comparison between values obtained from Trotter-Suzuki decomposition or symmetric cluster expansion , τ ∈ {0.04, 0.05, 0.06}, with HOTRG bond dimension D = 12, and values obtained from transfer-matrix renormalization-group method [23].

(43)

0.0 0.5 1.0 1.5 2.0 Temperature −1.0 −0.5 0.0 E ne rg y pe r si te

Energy per site for different external fields

SCE h = 0 T-S h = 0 SCE h = 1 T-S h = 1 SCE h = 2 T-S h = 2 SCE h = 2.5 T-S h = 2.5 TMRG h = 0 TMRG h = 1 TMRG h = 2 TMRG h = 2.5

Figure 4.6: Heisenberg spin-1/2 chain. Comparison between values obtained from Trotter-Suzuki decomposition or symmetric cluster expansion with HOTRG bond dimension D = 24, τ ∈ {0.04, 0.05, 0.06}, and values obtained from TMRG method [23].

0.0 0.5 1.0 1.5 2.0 Temperature 0.1 0.2 Mz /h

Ratio beteen Mz(T) and the external field h

SCE h = 0.05 T-S h = 0.05 SCE h = 1 T-S h = 1 SCE h = 2 T-S h = 2 SCE h = 2.5 T-S h = 2.5 TMRG h = 0.05 TMRG h = 1 TMRG h = 2 TMRG h = 2.5

Figure 4.7: Heisenberg spin-1/2 chain. Comparison between values obtained from Trotter-Suzuki decomposition or symmetric cluster expansion with HOTRG bond

(44)

4.2 Two dimensional systems

4.2.1 Transverse Ising and Heisenberg model

Other models used as benchmarks for the method are the two-dimensional transverse Ising and Heisenberg, with familiar Hamiltonians

HIsing = − X hi,ji σi_z σj_z− hX i σ_xi (4.3) Hheis = J X hi,ji ˆ Si· ˆSj (4.4)

where, for Ising, h is the strength of an external magnetic field and, for Heisenberg, J is the magnitude of antiferromagnetic interactions, which was fixed again to 1, while hi, ji denotes a sum over nearest neighbour sites.

For these two systems a network was built from a symmetric cluster expansion, with p = 4, and contracted using HOTRG and CTM. Then, results were compared with values from pure HOTRG [13] and Monte Carlo [24], and shown in Figure 4.8, 4.9 and 4.10, in which it is possible to see that relatively low bond dimensions D and χ are needed to already obtain data with a discrete level of accuracy.

0.2 0.6 1.0 Temperature −0.65 −0.50 −0.35 E ne rg y pe r sp in

Energy per spin Heisenberg 2D

HOTRG D = 14 & CTM χ = 16 Monte Carlo

Figure 4.8: Comparison between values obtained from a symmetric cluster expansion, τ = 0.0015, with vertical HOTRG and CTM, and values from Monte Carlo [24]. Note that MC error bars are small compared to the symbol sizes.

(45)

0 1 2 3 4 5 6 T 0.4 0.5 0.6 0.7 0.8 0.9 1.0 Mx

Spin transversal component of 2D transverse Ising model

HOTRG D=24 h=2 HOTRG D=24 h=3 HOTRG D=24 h=4 HOTRG D=16 & CTM χ = 16 h = 2 HOTRG D = 16 & CTM χ = 16 h = 3 HOTRG D = 16 & CTM χ = 16 h = 4

Figure 4.9: Comparison between values obtained from a symmetric cluster expansion, imaginary time step τ = 0.0015, with vertical HOTRG and CTM, and values from a Trotter-Suzuki decomposition with vertical and horizontal HOTRG [13].

0 1 2 3 4 5 6 7 8 T −4.0 −3.5 −3.0 −2.5 −2.0 −1.5 −1.0 E

Energy density f 2D transverse Ising m del

HOTRG D=24 h=2 HOTRG D=24 h=3 HOTRG D=24 h=4 HOTRG D=16 & CTM χ = 16 h = 2 HOTRG D = 16 & CTM χ = 16 h = 3 HOTRG D = 16 & CTM χ = 16 h = 4

(46)

4.2.2 Shastry-Sutherland model

Finally, the Shastry-Sutherland model was simulated in the thermodynamic limit obtaining values of energy E(T ) and z-component of the spin Sz(T ) per dimer, from

which specific heat C = ∂E_∂T and magnetic susceptibility χ = ∂Mz

∂h

_h=0, where h is a

small external field in the z-direction, have been computed using numerical deriv-atives. The coupling constant JD from 2.1 was fixed to one and results had been

obtained for increasing values of J , but always keeping the system in the dimer-product phase. Indeed, as expected, the computational complexity needed to repro-duce faithfully the data increases as the system moves towards its phase transition between dimer and plaquette phases. Moreover, as the latter regime is approached, errors from the finite machine precision become more relevant, making more challen-ging to obtain accurate results. These rounding errors can be reduced by introducing a lower cutoff in the singular values kept at each vertical HOTRG or CTM step, but they can also be reduced by just decreasing the number of iterations in the vertical coarse-graining procedure, even though in this way a fewer number of data points would be obtained.

For low values of J data were compared with results obtained from Monte Carlo simulations [5], which carry a small error in this regime. In Figure 4.11, 4.12 and 4.13 these are shown. Note that, when J is small, even for small values of horizontal HOSVD bond dimension D results are in agreement and the energy converges to the right ground state value of −3₄, however this is also a consequence of the initial setup used for either Trotter-Suzuki decomposition or symmetric cluster expansion, as we can see in Figure 3.19, which introduces a strong bias in the network. Indeed, in the dimer setup the ground state of the system needs only a horizontal bond dimension χ = 1 to be exactly reproduced. For the same reason the main error comes from the vertical coarse-graining scheme rather than the horizontal CTM, which need just a few steps to converge, and a small χ is actually enough to capture the relevant states even at finite temperature.

(47)

0.0

0.2

0.4

0.6

0.8

1.0 Temperature

−0.75

−0.60

−0.45

−0.30

E

(

T

)

Energy per dimer of the Sha try-Sutherland model for

J

=0.2

D

=4 (SCE)

D

=4 (T-S)

D

=8 (SCE)

D

=8 (T-S)

(a)

0.0

0.2

0.4

0.6

0.8

1.0 Temperature

−0.75

−0.60

−0.45

−0.30

E

(

T

)

Energy per dimer of the Sha try-Sutherland model for

J

=0.3

D

=4 (SCE)

D

=4 (T-S)

D

=8 (SCE)

D

=8 (T-S)

D

=12 SCE

D

=12 T-S

(b)

(48)

Trotter-0.0

0.2

0.4

0.6

0.8

1.0 Temperature

0.2

0.4

0.6

0.8

1.0

C

(

T

)

Specific heat per dimer of the Shastry-Sutherland model for

J

=0.2

D=4 (SCE) D=4 (T-S) D=8 (SCE) D=8 (T-S) MC (a)

0.0

0.2

0.4

0.6

0.8

1.0 Temperature

0.2

0.4

0.6

0.8

1.0

C

(

T

)

Specific heat per dimer of the Shastry-Sutherland model for

J

=0.3

D=4 (SCE) D=4 (T-S) D=8 (SCE) D=8 (T-S) D=12 SCE D=12 T-S MC (b)

Figure 4.12: Shastry-Sutherland specific heat density obtained from either Trotter-Suzuki decomposition or symmetric cluster expansion, for different values of HOTRG bond dimension D, and compared with results from Monte Carlo simulations [5]. Also, a CTM bond dimension χ = 8 was used. Note that MC error bars are small

(49)

0.0

0.5

1.0

1.5

2.0 Temperature

0.05

0.15

0.25

0.35

χ

(

T

)

Magnetic susceptibility per dimer of the Shastry-Sutherland model for J=0.2

D=4 (SCE) D=4 (T-S) D=8 (SCE) D=8 (T-S) MC (a)

0.0

0.5

1.0

1.5

2.0 Temperature

0.05

0.15

0.25

0.35

χ

(

T

)

D=4 (SCE) D=4 (T-S) D=8 (SCE) D=8 (T-S) MC (b)

Figure 4.13: Shastry-Sutherland magnetic susceptibility obtained from either Trotter-Suzuki decomposition or symmetric cluster expansion, for different values of HOTRG bond dimension D, and compared with results from Monte Carlo simu-lation [5]. Also, a CTM bond dimension χ = 8 and a small external magnetic field

(50)

Instead, for values of J near its critical value, data obtained from iPEPS, which is another method based on tensor network techniques [6], had been used as reference and are shown in Figure 4.14, 4.15 and 4.16. When J/JD approaches its critical

value, the data shows a concrete loss in accuracy due to the larger bond dimension needed to include all relevant states. Moreover, in this regime results can also be concretely affected by rounding errors, which increase for larger bond dimension D, and as a consequence only few vertical HOTRG steps can be performed.

0.0 0.2 0.4 0.6 0.8 1.0 Temperature −0.75 −0.60 −0.45 −0.30 E ( T )

Energy per dimer of the Sha try-Sutherland model for J=0.5

D=8 (SCE) D=8 (T-S) D=12 (SCE) D=12 (T-S) D=16 SCE D=16 T-S (a) 0.0 0.2 0.4 0.6 0.8 1.0 Temperature −0.75 −0.60 −0.45 −0.30 E ( T )

Energy per dimer of the Sha try-Sutherland model for J=0.6

D=8 (SCE) D=8 (T-S) D=12 (SCE) D=12 (T-S) D=16 SCE D=16 T-S (b)

Figure 4.14: Shastry-Sutherland energy density obtained from either Trotter-Suzuki decomposition or symmetric cluster expansion, for different values of HOTRG bond dimension D. Also, a CTM bond dimension χ = 16 was used.

(51)

0.0

0.2

0.4

0.6

0.8

1.0 Temperature

0.2

0.4

0.6

0.8

1.0

C

(

T

)

Specific heat per dimer of the Shastry-Sutherland model for

J

=0.5

D

=8 (SCE)

D

=8 (T-S)

D

=12 (SCE)

D

=12 (T-S)

D

=16 SCE

D

=16 T-S

iPEPS

(a)

0.0

0.2

0.4

0.6

0.8

1.0 Temperature

0.2

0.4

0.6

0.8

1.0

C

(

T

)

Specific heat per dimer of the Shastry-Sutherland model for

J

=0.6

D

=8 (SCE)

D

=8 (T-S)

D

=12 (SCE)

D

=12 (T-S)

D

=16 SCE

D

=16 T-S

iPEPS

(b)

(52)

0.0

0.5

1.0

1.5

2.0 Temperature

0.05

0.15

0.25

0.35

χ

(

T

)

D=8 (SCE) D=8 (T-S) D=12 (SCE) D=12 (T-S) D=16 (SCE) D=16 (T-S) iPEPS (a)

0.00

0.25

0.50

0.75

1.00 Temperature

0.05

0.15

0.25

0.35

χ

(

T

)

D=8 (SCE) D=8 (T-S) D=12 (SCE) D=12 (T-S) D=16 (SCE) D=16 (T-S) iPEPS (b)

Figure 4.16: Shastry-Sutherland magnetic susceptibility obtained from either Trotter-Suzuki decomposition or symmetric cluster expansion, for different values of HOTRG bond dimension D, and compared with results from iPEPS simulation [6]. Also, a CTM bond dimension χ = 16 and a small external magnetic field of strength h = 0.01 were used for these simulations.

(53)

4.3 Symmetric cluster expansion vs Trotter-Suzuki

decomposition

Comparisons between results obtained from an initial Trotter-Suzuki decomposition (T-S) and results from a symmetric cluster expansion (SCE) have been made in almost all the simulated models. For the one-dimensional Heisenberg model it is possible to see that SCE yields a smaller error than T-S, even though both methods provide results in close agreement. Instead, for the Shastry-Sutherland model it-self, Trotter-Suzuki produces results manifestly less accurate than SCE in almost all regimes of coupling constants. Therefore, it seems that the symmetric cluster expan-sion outperforms the Trotter-Suzuki decomposition for both one and two-dimenexpan-sional quantum systems, providing a more systematic way to represent operators, even though a SCE requires an initial bond dimension D which is in general higher than the one needed from a Trotter-Suzuki decomposition.

(54)

Chapter 5 Conclusions

The main goal of this project was to simulate the Shastry-Sutherland model at finite temperature T . In doing so, new tensor network methods have been tested and used to mainly perform two tasks: represent the density matrix operator e−β ˆH _{of the}

system as a tensor network, and realize an approximate contraction of this network in order to obtain thermal expectation values.

For the initial task two schemes were presented, a Trotter-Suzuki decomposition and a symmetric cluster expansion. The first method does not preserve, in general, all internal symmetries of the theory and produces an MPO, or PEPO, which is not hermitian. Instead, the second method provides objects that preserve all internal symmetries of the original Hamiltonian, as well as a systematically higher order error in τ = _Nβ.

Furthermore, for the second task, a higher order tensor renormalization group method was used to coarse-grain the network in the imaginary time direction, and it has been shown that HOTRG should not be used to coarse-grain also horizontal dimensions at low values of the inverse temperature β. Therefore, for two-dimensional quantum systems, a subsequent corner transfer matrix method was employed to contract a remaining two-dimensional network.

Results obtained from simulations of benchmark models are all in agreement with previous results, for both Trotter-Suzuki decomposition and symmetric cluster ex-pansion, and actual data of the Shastry-Sutherland model are almost in agreement with reference data, which are systematically approached with increasing values of

(55)

the truncated bond dimensions D and χ used in the simulations. Also, a neat im-provement in accuracy is guaranteed by the use of a symmetric cluster expansion, instead of a Trotter-Suzuki decomposition, to represent a D-dimensional density matrix operator as a (D+1)-dimensional tensor network.

Therefore, the symmetric cluster expansion has proven to be a novel remarkable tool, combined with higher order tensor renormalization group and corner transfer matrix methods, to simulate non-trivial quantum spin systems in the thermodynamic limit, keeping computational complexity reasonably low.

Its applications might help to shed light over highly non-trivial models and exotic quantum phases of matter, thanks to the high flexibility of this method.

(56)

References

[1] A. W. Sandvik, A. Avella and F. Mancini, ‘Computational studies of quantum spin systems’, 2010 (cit. on p. 1).

[2] Shastry and Sutherland, ‘Exact ground state of a quantum mechanical antifer-romagnet’, Physica B+C, vol. 108, no. 1, pp. 1069–1070, 1981 (cit. on pp. 1, 4, 5).

[3] R. Smith and D. Keszler, ‘Synthesis, structure, and properties of the orthob-orate SrCu2(BO3)2’, Journal of Solid State Chemistry, vol. 93, pp. 430–435,

Aug. 1991 (cit. on pp. 1, 5).

[4] M. Ueda, ‘Exact dimer ground state of the two dimensional heisenberg spin system SrCu2(BO3)2’, Phys. Rev. Lett., vol. 82, no. 18, 1999 (cit. on pp. 1, 5).

[5] S. Wessel, I. Niesen, J. Stapmanns, B. Normand, F. Mila, P. Corboz and A. Honecker, ‘Thermodynamic properties of the shastry-sutherland model from quantum monte carlo simulations’, Physical Review B, vol. 98, no. 17, Nov. 2018 (cit. on pp. 1, 42, 44, 45).

[6] A. Wietek, P. Corboz, S. Wessel, B. Normand, F. Mila and A. Honecker, ‘Ther-modynamic properties of the shastry-sutherland model throughout the dimer-product phase’, Physical Review Research, vol. 1, no. 3, Oct. 2019 (cit. on pp. 1, 46–48).

[7] J. Biamonte and V. Bergholm, Tensor networks in a nutshell, 2017 (cit. on p. 2).

[8] J. C. Bridgeman and C. T. Chubb, ‘Hand-waving and interpretive dance: An introductory course on tensor networks’, Journal of Physics A: Mathematical and Theoretical, vol. 50, no. 22, p. 223 001, May 2017 (cit. on p. 2).

[9] R. Penrose, ‘Applications of negative-dimensional tensors’, 1971 (cit. on p. 2). [10] U. Schollwöck, ‘The density-matrix renormalization group in the age of matrix product states’, Annals of Physics, vol. 326, no. 1, pp. 96–192, Jan. 2011 (cit. on pp. 2, 8).

[11] M. Suzuki, ‘Generalized Trotter’s formula and systematic approximants of exponential operators and inner derivations with applications to many-body problems’, Communications in Mathematical Physics, vol. 51, no. 2, pp. 183– 190, Jun. 1976 (cit. on pp. 3, 13).

(57)

[12] B. Vanhecke, L. Vanderstraeten and F. Verstraete, ‘Symmetric cluster expan-sions with tensor networks’, 2019 (cit. on pp. 3, 16).

[13] Z. Y. Xie, J. Chen, M. P. Qin, J. W. Zhu, L. P. Yang and T. Xiang, ‘Coarse-graining renormalization by higher-order singular value decomposition’, Phys-ical Review B, vol. 86, no. 4, Jul. 2012 (cit. on pp. 3, 22, 40, 41).

[14] T. Nishino and K. Okunishi, ‘Corner transfer matrix renormalization group method’, Journal of the Physical Society of Japan, vol. 65, no. 4, pp. 891–894, Apr. 1996 (cit. on pp. 3, 29).

[15] Introduction to Frustrated Magnetism. Lacroix, Mendels, Mila, 2011 (cit. on

p. 5).

[16] P. Corboz and F. Mila, ‘Tensor network study of the shastry-sutherland model in zero magnetic field’, Physical Review B, vol. 87, no. 11, Mar. 2013 (cit. on p. 6).

[17] F. Verstraete, V. Murg and J. Cirac, ‘Matrix product states, projected en-tangled pair states, and variational renormalization group methods for quantum spin systems’, Advances in Physics, vol. 57, no. 2, pp. 143–224, Mar. 2008 (cit. on p. 8).

[18] Numerical recipes. Press, Teukolsky, Vetterling, Flannery, 2007 (cit. on p. 9).

[19] L. Lathauwer and B. De Moor, ‘A multi-linear singular value decomposition’, Society for Industrial and Applied Mathematics, vol. 21, pp. 1253–1278, Mar. 2000 (cit. on p. 11).

[20] M. Levin and C. P. Nave, ‘Tensor renormalization group approach to two-dimensional classical lattice models’, Physical Review Letters, vol. 99, no. 12, Sep. 2007 (cit. on p. 22).

[21] P. Corboz, S. R. White, G. Vidal and M. Troyer, ‘Stripes in the two-dimensional t-j model with infinite projected entangled-pair states’, Physical Review B, vol. 84, no. 4, Jul. 2011 (cit. on p. 32).

[22] P. Corboz, T. M. Rice and M. Troyer, ‘Competing states in thet-jmodel: Uniformd-wave state versus stripe state’, Physical Review Letters, vol. 113, no. 4, Jul. 2014 (cit. on p. 32).

[23] T. Xiang, ‘Thermodynamics of quantum heisenberg spin chains’, Physical Re-view B, vol. 58, no. 14, pp. 9142–9149, Oct. 1998 (cit. on pp. 38, 39).

[24] A. Cuccoli, V. Tognetti, R. Vaia and P. Verrucchi, ‘Two-dimensional quantum heisenberg antiferromagnet: Effective-hamiltonian approach to the thermody-namics’, Physical Review B, vol. 56, no. 22, pp. 14 456–14 468, Dec. 1997 (cit. on p. 40).

Thermodynamics of the Shastry-Sutherland model using tensor network methods