Simulation of the two-dimensional Ising model using the Corner
Transfer Matrix Renormalisation Group method on the square,
honeycomb and triangular lattice
Fenna van der Ploeg July 3, 2018
Student-id 10799907
Study program Bachelor Physics and Astronomy
Course Bachelor Project Physics and Astronomy, 15 EC
Conducted between 03-04-2018 and 03-06-2018
Supervisor dr. P. R. Corboz
Second Assessor dr. E. Lerner
Abstract
To simulate the two-dimensional Ising model, its partition function is represented as a tensor network. This network is contracted using the Corner Transfer Matrix Renormalization Group (CTMRG) method, in order to calculate the magnetization m, free energy f and correlation length ξ as a function of temperature. This method is first applied to the square lattice and then a generalization of the CTMRG method to the honeycomb and triangular lattices is presented.
Using finite size scaling, the critical temperatures Tc of the different lattice types are
determined in three different ways. All values for Tc found are in agreement with the exact
values. Also values for the critical exponents β and ν are found using finite size scaling. ν is determined to be ν = 0.97 ± 0.04, ν = 1.01 ± 0.07 and ν = 0.99 ± 0.03, for the square, honeycomb and triangular lattice, respectively. These values are in agreement with the exact value of ν = 1. For the ratio of the critical exponents β/ν, the values obtained are β/ν = 0.1249 ± 0.0003 for the square lattice, β/ν = 0.1249 ± 0.0004 for the honeycomb lattice and β/ν = 0.1254 ± 0.0005 for the triangular lattice. Also these values agree with the exact value of β/ν = 1/8.
As all results obtained agree with their exact values, it can be concluded that the CTMRG is a valid method to simulate the partition function not only for the square lattice, but also for the honeycomb and triangular lattices.
Populair wetenschappelijke samenvatting
Faseovergangen kennen we allemaal: je ziet bijvoorbeeld water verdampen tijdens het koken van je pasta en als het ’s winters koud genoeg is, bevriezen de sloten. De temperatuur waarbij zo’n overgang tussen twee verschillende fases plaatsvindt noemen we de kritische temperatuur. Als je de dichtheid van het water bij verschillende temperaturen zou meten en hier een grafiek van zou maken, zou je zien dat de dichtheid een sprong maakt bij de kritische temperatuur.
Een van de meest eenvoudige modellen die gebruikt wordt in de statistische fysica om faseovergangen te bestuderen en te begrijpen is het Ising model. Dit theoretische model gaat uit van een versimpelde beschrijving van een materiaal: het bestaat uit een rooster met op ieder roosterpunt een klein magneetje. Dit magneetje noemen we spin en kan ofwel omhoog, ofwel omlaag wijzen (en dus niet ergens daartussen in!).
Bij hoge temperaturen zorgen thermische fluctuaties ervoor dat de spins een random kant op wijzen; de helft van de spins wijst omhoog, de helft omlaag. Bij lage temperaturen ordenen de spins zich. De meeste spins wijzen dan in dezelfde richting, omdat dit de laagst mogelijke energietoestand geeft. Net zoals het bevriezen van water, is de transitie van de spins in ongeordende toestand naar geordende toestand een faseovergang. In plaats van de dichtheid verandert bij de kritische temperatuur nu de magnetisatie (deze keer niet met een sprong maar gelijdelijker; dit proces is een zogenaamde tweede-orde-faseovergang).
Er zijn verschillende roosters waarop de spins geplaatst kunnen worden; in figuur 1 zijn drie van die typen roosters ge¨ıllustreerd. Voor het Ising model op een vierkant rooster is een analytische oplossing bekend: in 1944 vond Lars Onsager een formule die de magnetisatie van het systeem als functie van de temperatuur beschrijft. Voor andere roosters, zoals een honingraatrooster en een driehoekig rooster, bestaat er echter geen exacte oplossing. Om deze systemen toch te bestuderen, kun je gebruik maken van computersimulaties.
Figuur 1: De drie verschillende roostertypen die bestudeerd zijn in dit onderzoek: het vierkante rooster is weergegeven in figuur (a), het honingraatrooster is weergegeven in figuur (b) en het driehoekig rooster is weergegeven in figuur (c).
Zo’n computersimulatie kun je op verschillende manieren doen. De methode die ik heb gebruikt heet de Corner Transfer Matrix Renormalization Group method (afgekort CTMRG). In deze methode wordt het systeem beschreven als een netwerk van tensoren. Een tensor is een wiskundig object waarin je getallen kunt opslaan. In feite is het niets meer is dan een generalisatie van een matrix; een matrix heeft twee dimensies, een tensor kan meer dimensies hebben.
Met de CTMRG simulatie worden eindige tensoren gevonden die het oneindige systeem het beste benaderen. Deze kunnen vervolgens gebruikt worden om thermodynamische groot-heden, zoals de magnetisatie, te berekenen. Ik heb dit gedaan voor het Ising model op een vierkant rooster en de resultaten vergeleken met de exacte oplossing. Deze kwamen goed overeen, zoals ook al was aangetoond in eerder onderzoek. Vervolgens heb ik de CTMRG gegeneraliseerd voor het Ising model op een honingraat- en driehoekig rooster. Ook voor deze roostertypen gaf de CTMRG methode goede resultaten. De generalisatie die ik heb ge¨ıntroduceerd is dus valide.
Contents
1 Introduction 6
2 Theory 7
2.1 The Ising model . . . 7
2.2 Tensor networks . . . 8
2.3 Tensor network representation of the Ising model . . . 10
2.4 Critical exponents . . . 12
2.5 Finite size scaling . . . 13
3 Methods 15 3.1 The Corner Transfer Matrix Renormalization Group method . . . 15
3.2 Modifications for other lattice types . . . 18
3.2.1 Honeycomb lattice . . . 18
3.2.2 Triangular lattice . . . 22
3.3 Calculating thermodynamic quantities . . . 25
3.3.1 Magnetization . . . 26 3.3.2 Free energy . . . 26 3.3.3 Correlation length . . . 27 4 Results 28 4.1 Thermodynamic quantities . . . 28 4.1.1 Magnetization . . . 28 4.1.2 Free energy . . . 31 4.1.3 Correlation length . . . 33
4.2 The critical temperature and critical exponents . . . 35
1
Introduction
In theoretical condensed matter physics, theoretical models are studied to describe and un-derstand the properties of condensed matter. A famous example of such a model is the Ising model, first introduced in 1920 (Lenz, 1920). Its two-dimensional case is one of the simplest statistical models that provides insight into phase transitions.
In the Ising model, the magnetic dipole moments of atomic spins are represented by a discrete value of +1 or -1. The one- and two-dimensional cases of the Ising model (on a line and on a square lattice, respectively) can be solved analytically (Onsager, 1944). In three or more dimensions, however, the model has not been solved exactly, nor was the two-dimensional Ising model on a honeycomb or triangular lattice (Taroni, 2015).
To study these analytically unsolvable models, different numerical methods can be utilized, like Monte Carlo simulations (Boer, 2011) and the Corner Transfer Matrix Renormalization Group method. In this thesis, the accuracy and validity of the latter method will be studied. The Corner Transfer Matrix Renormalization Group (CTMRG) method for the square lattice was proposed by Nishino & Okunishi (1995), unifying the Corner Transfer Matrix method (Baxter, 1978) and the Density Matrix Renormalization Group method (White, 1992). The method is based on the idea that the partition function of a system can be represented by a so called tensor network. The CTMRG provides an efficient way to contract this network, i.e. to calculate the partition function.
In this thesis, I will confirm the original CTMRG method for the square lattice, finding results in agreement with previous studies. Moreover, I will generalize the CTMRG and apply it to the honeycomb and triangular lattices, demonstrating the validity of this method for these lattice types. For the three different lattices studied, different thermodynamic quanti-ties will be calculated and compared to their analytical solutions. Furthermore, a value for the critical exponents β and ν and critical temperatures of these lattices will be determined.
This thesis is organized as follows. In section 2, I will provide theoretical background on the Ising model and tensor networks. The Corner Transfer Matrix Renormalization Group method will be explained in section 3, as well as the methods used to calculate the different observables. In section 4, the results are presented, which will be discussed in section 5.
2
Theory
2.1 The Ising model
In 1920 Wilhelm Lenz introduced the Ising model. In this infinite lattice model for ferromag-netism, each lattice point has a discrete value of +1 or -1. This value represents the magnetic dipole moment of an atomic spin that is either up or down (Lenz, 1920). The model was named after Lenz’ student Ernst Ising, who solved the one-dimensional case of this model and proved it has no phase transition (Ising, 1925). The two-dimensional Ising model, on the other hand, does have a phase transition; an analytic description of this was provided by Onsager (1944).
In the Ising model, each spin interacts with its direct neighbors; two adjacent sites i, j have an interaction energy Jij. Moreover, each site interacts with its external magnetic field
hi. Thus, the Hamiltonian of the system can be written as
H(s) = −X hi,ji Jijsisj− µ X i hisi, (1)
where hi, ji denotes that sites i and j are nearest neighbors; every pair is counted once. µ is the magnetic moment, which I will set to 1.
If for every pair hi, ji the interaction energy Ji,j is positive, the interaction is called
ferromagnetic. Conversely, if the interaction energy Ji,j is negative for every pair hi, ji, the
interaction is called antiferromagnetic.
Assuming that the interaction strength is the same for every pair of nearest neighbors and that the external magnetic field is the same for every site, the Hamiltonian becomes
H(s) = −JX hi,ji sisj− h X i si. (2) Using β = k1
BT, this Hamiltonian results in a canonical partition function of
Z(β) =X s e−βH(s)=X s eβJ P hi,jisisj+βhPisi =X s Y hi,ji eβJ sisjY i eβhsi, (3)
where T is the temperature of the system and kB is Boltzmann’s constant. Computing the
partition function by evaluating this sum directly scales exponentially with system size. The partition function enables us to find the magnetization m of the system:
m = hsri = 1 Z X s sr Y hi,ji eβJ sisjY i eβhsi. (4)
Moreover, with the partition function we can determine the free energy f of the system, which is defined as
f = −kBT ln(Z)
N , (5)
where N is the total number of lattice sites.
Lars Onsager (1944) provided an analytic description for the two-dimensional Ising model on a square lattice with zero external magnetic field, h = 0. He proposed a solution for the free energy and the magnetization, but did not provide prove of the latter. A complete derivation was later given by Yang (1952).
The exact solutions for the magnetization and the free energy are −βf = ln(2) 2 + 1 2π Z π 0 ln cosh2(2βJ ) +1 k p 1 + k2− 2kcos(2θ), (6)
with k = sinh−2(2βJ ), and
m(T ) = [1 − sinh−4(2βJ )]18 T < Tc 0 T > Tc. (7)
Tc is called the critical temperature; it is the temperature at which a phase change of the
system occurs. For temperatures lower than Tc, spontaneous magnetization appears. In this
phase most of the spins are aligned (either in the +1 or -1 state), as this generates the lowest energy states. For temperatures higher than Tc, thermal fluctuations cause the system to
become disordered. In this phase, approximately half of the spins is in the +1 state and half of the spins is in the -1 state, resulting in no magnetization.
The exact value of Tc for the two-dimensional Ising model was derived exactly (Kramers
& Wannier, 1941): Tc= 1 kBβc = 2 pln(1) + 2 ≈ 2.26918531421 J/kB (8)
for a system with J = 1 and h = 0.
2.2 Tensor networks
In the field of numerical simulation algorithms, there has been increasing interest in tensor network methods to simulate strongly correlated systems (Verstraete et al., 2008). A tensor is a multidimensional array of (complex) numbers, which can be represented with a diagram-matic notation, as is illustrated in figure 2. In this notation, a tensor is represented by a closed shape and its indices are represented as lines emanating from it. As depicted in figure 2, scalars, vectors and matrices can also be represented using this notation. They can be interpreted as tensors with zero, one and two indices, respectively.
Figure 2: Tensor network diagrams representing a scalar, a vector, a matrix and a tensor.
If a line connects two tensors, the corresponding index can be contracted. This is done by summing over all the possible values of the common index. For example, the summation of Pαβ and Oβγ over β results in new matrix Rαγ, as shown in figure 3 in both diagrammatic
notation and its equivalent equation.
Figure 3: Diagrammatic notation of the contraction of an index, here representing a matrix multipli-cation.
A tensor network consists of multiple tensors with connected indices. In figure 4, an example of such a network N with four open indices, α, β, δ and γ is given. Contracting the network (that is, summing over all connected indices), results in a new tensor T that has four open indices α, β, δ and γ as well.
Figure 4: Diagrammatic notation of a tensor network N with four open indices, α, β, δ and γ. After contracting the network, the resulting tensor T has four open indices as well.
2.3 Tensor network representation of the Ising model
To represent the partition function of the Ising model as a tensor network, I will first rewrite the partition function. For a two-dimensional square lattice, the product over lattice sites can also be written as a product over nearest neighbors:
Y
i
eβhsi = Y
hi,ji
eβh(si+sj)/4, (9)
where the division by 4 in the exponent compensates for the fact that on a square lattice, every site has 4 nearest neighbors. Now, the partition function can be rewritten as
Z(β) =X s Y hi,ji eβJ sisjY i eβhsi =X s Y hi,ji eβJ sisj+βh(si+sj)/4≡X s Y hi,ji Qsi,sj, (10)
where Qsi,sj is the transfer matrix, which can be written as a 2 × 2 matrix. This becomes
clear when writing out the different possible values for Q explicitly. Since the spins in the Ising model can be either up or down, siand sj can take the discrete value of ±1. This means
there are four different combinations of si and sj, which result in the following four values for
Qsi,sj: Q+1,+1= eβ(J +h/2) Q−1,+1= e−βJ Q+1,−1= e−βJ Q−1,−1= eβ(J −h/2). (11)
Next, writing Q as a matrix, we get: Q = " Q+1,+1 Q+1,−1 Q−1,+1 Q−1,−1 # = " eβ(J +h/2) e−βJ e−βJ eβ(J −h/2) # . (12)
With this expression for matrix Q, the partition function can be represented as an infinite network. On each bond between two lattice sites a tensor Q is placed, and on each lattice site a four dimensional Kronecker delta function δijklis placed, where δijkl= 1 if i = j = k = l and
δijkl = 0 otherwise. In figure 5a, a part of the resulting lattice is shown. The delta function
represents the possible states the lattice site can be in, and Q contains the local Boltzmann weights of the four possible local spin configurations. Contracting the tensor network, that is, summing over all indices, gives all the possible configurations of the model, with the proper weight.
Figure 5: The diagrammatic representation of the partition function of the two dimensional square lattice Ising model as a tensor network. In figure (a), a δ-function is placed on each lattice site. Between all the lattice sites, the transfer matrix Q is placed, which contains the local Boltzmann weights. Then all Q’s are split into two identical matrices √Q, shown in (b). Each δ-tensor is then contracted with four√Q-tensors, as indicated by the blue circles in (c). The resulting lattice, as shown in (d), consist solely of identical tensors a.
Next, the representation is simplified into a network with only one type of tensor. To do this, Q is written as the contraction of two identical tensors√Q:
Qsi,sj = X k p Qs i,k· p Qk,s j. (13)
In figure 5b, this replacement of Q in the network is given in diagrammatic notation. Next, each Kronecker delta tensor δijkl is contracted with four tensors
√ Q: aijkl= X s p Q is p Q js p Q ks p Q ls, (14)
which results in a network consisting of only one type of tensor aijkl, shown in figure 5d.
Like the partition function, equation 4 for the magnetization can be rewritten using equation 9, obtaining m(β) = 1 Z X s sr Y hi,ji Qsi,sj. (15)
The numerator of this equation, P
ssr
Q
hi,jiQsi,sj, can be displayed as a tensor network
similar to the network for Z, except that at lattice site r, the tensor aijkl is replaced by a
tensor bijkl that measures the spin on that site:
bijkl = X sr sr p Q isr p Q jsr p Q ksr p Q lsr . (16)
Thus, as shown in diagrammatic notation in figure 6, the magnetization m can be displayed as a division of two tensor networks.
Figure 6: The expression for the magnetization in diagrammatic notation. It is a division of two infinite networks that are identical except for one tensor: in the network in the numerator, one a-tensor is replaced by tensor b, as defined in equation 16.
2.4 Critical exponents
The Ising model shows a second order (i.e. continuous) phase transition; at a certain temper-ature, called the critical temperature Tc, the system undergoes a phase transition between an
exhibit a power law behavior characterized by critical exponents. To express the behavior of thermodynamic quantities in terms of a power law around the critical point, we define the reduced temperature t as
t = T − Tc Tc
, (17)
where T is the temperature of the system and Tc is the critical temperature. At the phase
transition, T = Tc and thus the reduced temperature is zero.
The magnetization m and the correlation length ξ obey a power law in the reduced temperature:
m ∼ tβ and (18)
ξ ∼ |t|−ν. (19)
Here, β and ν are the critical exponents. Note that the former one should not be confused with the inverse temperature β = 1/kBT .
Critical exponents are universal; they only depend on the dimension and symmetry of the system, not on its microscopic details. Thus, the values of β and ν should not depend on the type of lattice. For the two-dimensional Ising model, they are known exactly: β = 1/8 and ν = 1 (Yeomans, 1992).
2.5 Finite size scaling
When using a finite system of size L to simulate an infinite system, systematic errors called finite size effects can occur.
Away from the critical point, the correlation length ξ is finite. Here, finite size effects are negligible for systems with L ξ. As T approaches Tc, however, the correlation length
diverges. In practical simulations, usually a finite system size L is used. This imposes a cutoff on ξ; the correlation length can not grow larger than L. This results in an order parameter that is smeared out around Tc, rather than showing a sharp transition.
Finite size effects can be exploited by incorporating L as a parameter in the scaling functions for thermodynamic quantities. This approach is called finite size scaling. The cutoff L imposes on ξ affects the system as a finite distance from the critical point would, implying L ∼ ξ ∝ |t|−ν. This relation can be rewritten as
T∗(L) = a · L−1/ν+ Tc, (20)
where T∗(L) is the critical temperature at system size L, Tc is the critical temperature at
Finite size scaling can also be applied to other quantities of interest. For an observable X that diverges as t−χ for L ξ, the following ansatz is made (Fisher & Barber, 1972):
X(t, L) = |t|−χF1(L/ξ), (21)
where F1 is an unknown function. Using |t| ∝ ξ−1/ν and L/ξ ∝ tνL, this can be rewritten as
X(t, L) = |t|−χF2(Lχ/ν · L1/ν). (22)
For the magnetization, this becomes:
m(t, L) = L−β/νF (t · L1/ν), (23)
Where F is again an unknown function. Taking the natural logarithm of this equation at t = 0 gives
ln[m(t = 0, χ)] = −β
νln(L) + C, (24)
where C is an arbitrary constant. If Tcis known, this equation makes it possible to obtain the
ratio β/ν by finding m(t = 0, χ) for different system sizes and performing a fit to the data. Moreover, equation 23 can be used to perform a data collapse. Rewriting it gives
m(t, L) · Lβ/ν = F (t · L1/ν) or (25)
y = F (x), (26)
where I have defined y = m(t, L) · Lβ/ν and x = t · L1/ν. For the right critical exponents, all
points for different system sizes L should fall on one curve. Thus, the critical exponents and critical temperature can be obtained by plotting y = m(t, L) · Lβ/ν against x = t · L1/ν for different t and L and finding β, ν and Tc such that the data fall on a single curve.
3
Methods
3.1 The Corner Transfer Matrix Renormalization Group method
To approximate the infinite tensor network of the partition function with a finite tensor network, I run a simulation based on the Corner Transfer Matrix Renormalization Group (CTMRG) method. In this simulation, I consider the Ising model without an external mag-netic field (h = 0) and set the interaction energy and Boltzmann’s constant to 1: J = 1, kB= 1.
Figure 7: The ansatz of the Corner Transfer Matrix Renormalization Group Method: the infinite tensor network that represents the partition function (shown on the left) can be approximated by a finite tensor network, where each corner tensor represents one quadrant of the infinite lattice and each edge tensor replaces an infinite half-row or half-column of the lattice (shown on the right).
Considering the square lattice, the partition function can be expressed as an infinite network of tensors a, as shown in figure 7 on the left. The CTMRG method is based on the ansatz that the environment of one bulk tensor a can be approximated by four corner tensors C and four edge tensors T. As shown in figure 7, each corner tensor represents one quadrant of the infinite lattice and each edge tensor replaces an infinite half-row or half-column of the lattice. The higher the bond dimension of these edge and corner tensors, the more accurate this approximation is; for infinite bond dimensions, it is exact.
Since the square lattice is rotationally and translationally symmetric, all four corner ten-sors are the same and so are all four edge tenten-sors. Thus, only one corner tensor C and one edge tensor T need to be found in the simulation. Moreover, both C and T are symmetric under exchange of legs, implying Cij = Cji and Tijk = Tjik, where index i and j have the
same bond dimension.
The tensors C and T that approximate the partition function Z at a given temperature are obtained trough iteration. To start, the algorithm initializes C and T with random values
and bond dimensions (χ, χ) and (χ, χ, d), respectively. In each iteration, the system is grown from size n×n to size (n+2)×(n+2) by inserting eight bulk tensors a and eight edge tensors T, step 1 in figure 8.
Figure 8: The CTMRGM. Each iteration starts with a tensor network as shown in the upper left corner. Then, in step 1, the system is grown by inserting eight bulk tensors a and eight edge tensors T. In step 2, the renormalization tensors are inserted. The network is divided into nine parts (step 3), which are contracted (step 4) to get a tensor network with improved corner and edge tensors. These steps are repeated until convergence is reached.
In principle, the new corner tensor C0 can now be found by contracting the new corner (consisting of tensors C, T, T and a) and the new edge tensor T0 can be found by contracting the edge (consisting of tensors T and a). However, this would make the bond dimensions of C and T grow with a factor d in each iteration, d being the number of possible states for each spin. This growth would cause computational time to exponentially increase with the number of iterations. Therefore, when the bond dimensions exceed a set maximum bond dimension χ, a renormalization on C and T is done. To do this, renormalization tensors are inserted in the network before the edges and corners are contracted, step 2 in figure 8. This limits the bond dimension to χ.
The process preceding the renormalization, to find the right renormalization tensors, is illustrated in figure 9. First, the tensors C, T, T and a are contracted into a tensor M of size χd × χd. Next, we would like to find the tensor ˘M with size χ × χ that best approximates M . To find ˘M , a singular value decomposition on M is done. This is a factorization of M of the form U · s · VT. Here, U and VT are unitary matrices and s is a diagonal matrix with
as its diagonal entries the singular values of M in descending order. To get ˘M , only the χ largest singular values are kept; tensor U is truncated from size χd × χ × d to size χ × χ × d, giving ˜U . ˜U† is then used as the normalization tensor shown in figure 8.
Figure 9: The construction of the normalization tensor. First, the tensors C, T, T and a are contracted into a tensor M . Then, a singular value decomposition on M is done, which factorizes M into U ·s·VT.
The tensor U is then truncated from size χd × χ × d to size χ × χ × d, giving ˜U . ˜U† is then used as the normalization tensor.
In the field of quantum many-body physics, a similar procedure called the density ma-trix renormalization group (DMRG) is used (White, 1992). This method is useful because in quantum many body-systems, Hilbert space grows exponentially with system size. For example, a spin-12 chain of length K has 2K degrees of freedom. DMRG reduces the degrees of freedom to the most important ones.
The new boundary tensors C0 and T0, with sizes χ × χ and χ × χ × d respectively, are found by contracting the corners and edges with ˜U†, as shown in figure 10. To keep numbers bounded, C0 and T0 are divided by their largest element. The obtained tensors are then used in the next iteration.
The leading computational cost of this algorithm is O(χ3), its memory cost is O(χ2).
Figure 10: The production of the new corner tensor C0and edge tensor T0. The renormalization tensor ˜
U† ensures C0 and T0 are of size χ × χ and χ × χ × d, respectively.
To check for convergence, the normalized singular value spectrum sk of corner tensor C0 is
is below a given tolerance τ :
X
k
|sk− s0k| < τ. (27)
Once convergence is reached, C and T can be used to calculate the thermodynamic quantities of interest.
3.2 Modifications for other lattice types
For the honeycomb and triangular lattices, not all four corner and edge tensors are the same. This requires some modifications to the Corner Transfer Matrix Renormalization Group method for the square lattice as outlined in section 3.1. The infinite tensor network of the partition function is still approximated by a tensor network consisting of nine tensors, but in this case I give all corners and edges a different label, as shown in figure 11.
Figure 11: The finite network used in the CTMRG method, this time with four different corner and edge tensors.
3.2.1 Honeycomb lattice
In figures 12a and 12b, the construction of the tensor network for a honeycomb lattice is shown. Here, bulk tensor a is similar to the one used for the square lattice, except for its size. Now, instead of the contraction of four√Q-tensors and one four-dimensional Kronecker delta tensor δijkl, it is a contraction of three
√
Q-tensors and one three-dimensional Kronecker delta tensor δabc:
aijk= X a,b,c δabc( p Q)ai( p Q)bj( p Q)ck . (28)
Figure 12: Mapping the honeycomb lattice onto a square lattice. To represent the honeycomb lattice as a tensor network, a δ-function is placed on each lattice site. Between all the lattice sites, the transfer matrix Q is placed, which contains the local Boltzmann weights. Then all Q’s are split into two identical matrices√Q. Each δ-tensor is then contracted with four√Q-tensors, as indicated by the blue circles in (a). In the resulting lattice (b), each pair of tensors is contracted to a four dimensional tensor ˜a. This results in a square lattice, shown in (d).
To apply the Corner Transfer Matrix Renormalization Group method on the honeycomb lattice, it is mapped onto a square lattice. In figures 12c and 12d it is illustrated how this is done. The a-tensors are contracted in pairs to tensors ˜a. The result is a square lattice network consisting of four-dimensional tensors ˜a.
Bulk tensor ˜a is not fully symmetric, as the bulk tensor for the square lattice was. There-fore, not all four corner tensors are equal. There are, however, still some symmetries we can exploit. In figure 13, ˜a is depicted as its original two a-tensors. In this way, the symmetries of the corners can be well observed; it is clear that for this lattice type C1 = C3 and C2= C4.
As can be deduced from figure 13b, all four edge tensors are the same: T1 = T2 = T3 = T4
still applies. As is the case for the square lattice, the corner tensors are symmetric under exchange of legs: Cµij = Cµji, where µ = 1, 2, 3 or 4. The legs of the edge tensors, however,
(a) (b)
Figure 13: The corner and edge tensors for the honeycomb lattice (before renormalization is done). To illustrate the symmetries of the system, ˜a is depicted as its original two a-tensors. From (a), it can be inferred that C1= C3 and C2= C4. From (b) it can be deduced that T1 = T2 = T3 = T4 still
applies.
Because the system shows less symmetries, instead of ˜U†, two different renormalization tensors are used, called Pu and Pd. In figure 14, it is shown how these tensors are inserted in
the network.
Figure 14: An iteration of the CTMRG for the honeycomb lattice. Contrary to the square lattice, two different renormalization tensors are used, indicated by two different colors.
Figure 15: The construction of the two renormalization tensors Pu and Pd. In (a), C1 and C4 are
contracted and a singular value composition is done, giving U , V and s. Since V†, U† and s−1 would be the result of a singular value decomposition on the contraction of C−11 and C−14 (shown in (b)), a contraction of C4, V†, s−
1
2,s−
1
2, U† and C1 approximates 1, as shown in (c). This network is then
split, as shown in (d), to get the renormalization tensors Puand Pd.
In figure 15, the construction of the renormalization tensors is depicted. First, C1 and C4
are contracted and a singular value composition is done. The resulting tensors U , V and s are then used to construct Pu and Pd. As shown in figure 15d, the former is a contraction of
C4, V†and s−
1
2. The latter is a contraction of C1, U† and s− 1 2.
Since C4 = C2, a similar execution with C1 and C2 (instead of C1 and C4) results in
the same renormalization tensors Pu and Pd. Thus, Pu is used to renormalize C1 on both
sides and Pd is used for the renormalization of C4, as shown in figure 16a and 16b. The
Figure 16: The production of the new corner tensors C01 = C03 and C04 = C02 and edge tensors T0
1= T02= T03= T04for the honeycomb lattice.
3.2.2 Triangular lattice
In figure 17a, the tensor network representation of the triangular lattice is illustrated. It consists again of two-dimensional Q-tensors and, this time six-dimensional, Kronecker delta functions δijklmn. To apply the Corner Transfer Matrix Renormalization Group method on
this type of lattice, it is first mapped to a honeycomb lattice. To do this, δijklmn is written
as a tensor network consisting of four tensors δijk: δijklmn=
X
a,b,c
δaijδbklδabcδcnm, (29)
and this network is used to replace δijklmn, as illustrated in figure 17b in diagrammatic
notation.
Next, a contraction of 3 δ-tensors and 3 Q-tensors is done, resulting in a new tensor e, shown in figure 17e:
epqr =
X
i,j,k,l,m,n
QijQklQmnδpinδjqkδlmr (30)
Contracting all Q-tensors in this way results in a honeycomb lattice with, as can be seen in figure 17f, two different tensors: e and δabc.
Figure 17: The tensor network representation of the partition function for the triangular lattice. In figure (a), a δ-function is placed on each lattice site. Between all lattice sites, the transfer matrix Q (in yellow) is placed, which represent the interaction between the sites. Then, as shown in (c), each six dimensional δ-function is replaced by a network of four three dimensional δ-functions as illustrated in (b). In the new network (d), a group of three Q-tensors and three δ-tensors is contracted into a new tensor e, as in (e). This results in a honeycomb lattice with two different tensors, as displayed in (f).
Figure 18: The bulk tensor ˜a as used for the triangular lattice. It is the result of a contraction of e and δabc.
The resulting honeycomb lattice can be mapped onto a square lattice as before. To do this, the bulk tensors ˜a are constructed as demonstrated in figure 18.
This definition of ˜a shows less symmetries than the bulk tensor used for the honeycomb lattice. Looking at figure 19a, we see that C1 = C3 still applies, but C4 6= C2. Moreover,
from figure 19b we see that T1 = T2 and T3 = T4, but T16= T3. Thus, we need to find
C1, C2, C4, T1 and T3.
For this lattice type, Cµij = Cµji, still applies for µ = 2 or 4. However, the other two
corner tensors are no longer symmetric: Cνij 6= Cνji for ν = 1 or 3. Neither are the edge
tensors: Tijkω 6= Tjikω for ω = 1, 2, 3 or 4.
(a) (b)
Figure 19: The corner and edge tensors for the triangular lattice (before renormalization is done). To illustrate the symmetries of the system, ˜a is depicted as its original two tensors e and δ. From (a), it can be inferred that C1= C3 still applies but C26= C4. From (b) we see that T1= T2 and T3= T4,
but T16= T3.
For the renormalization, Pu and Pd are again constructed as described before. However,
C2 are no longer equal to Pu and Pd. Using the same procedures as outlined in figure 15, the
renormalization tensors Pl and Pr are found, as can be seen in figure 20. Pu, Pd, Pl and Pr
are used to renormalize the corners and edges as depicted in figure 21.
Figure 20: The construction of renormalization tensors Pland Pr, similar to the construction of Pu,
Pd as detailed in figure 15.
Figure 21: The production of the new corner tensors C01= C03(a), C02(b) and C04(c) and edge tensors
T04= T03(d) and T01= T02(e) for the triangular lattice.
3.3 Calculating thermodynamic quantities
With the corner and edge tensors obtained from the simulation using the CTMRG method, different thermodynamic quantities can be determined. Here, the magnetization m, free
energy f and correlation length ξ will be calculated.
3.3.1 Magnetization
The tensor network representation of the magnetization, as shown in figure 6, can be approx-imated by the division of two tensor finite networks, as shown in figure 22. As in figure 6, the tensors a and b correspond to equations 14 and 16 respectively. The corner and edge tensors are obtained through simulation based on the CTMRG method as outlined in sections 3.1 and 3.2.
Figure 22: The magnetization m approximated by the division of two tensor networks. As in figure 6, the tensors a and b correspond to equations 14 and 16 respectively. The corner and edge tensors C and T are obtained through simulation based on the CTMRG method as outlined in sections 3.1 and 3.2.
3.3.2 Free energy
To calculate the free energy for a given temperature, as given by equation 5, first the partition function per site κ is determined:
κ = lim
N →∞Z
1/N, (31)
which can be approximated by a multiplication and division of three different tensor networks, as shown in figure 23 (Chan, 2013). Subsequently, κ is used to calculate the free energy:
f = −T · ln(κ), (32)
Figure 23: The free energy f approximated by a multiplication and division of three different tensor networks.
3.3.3 Correlation length
Using the edge tensors, another thermodynamic quantity can be calculated, namely the cor-relation length ξ. This is a measure of the characteristic distance over which the behavior of one spin is correlated to the behavior of another spin.
The correlation length ξ at a given temperature is found by first contracting two edge tensors T to a matrix TM, as illustrated in figure 24. Then, the two largest eigenvalues λ1
and λ2 of TM are determined and used to find the correlation length, according to
ξ = 1
ln(λ1/λ2)
, (33)
where λ1 > λ2 (Yeomans, 1992).
Figure 24: The contraction of two edge tensors T to a matrix TM. The two largest eigenvalues λ1and
4
Results
In this section, I present the simulation results of the different thermodynamic quantities, followed by the finite size scaling results.
4.1 Thermodynamic quantities
4.1.1 Magnetization
In figure 25, the simulation results for the magnetization of the square lattice as a function of temperature are presented. For temperatures ranging from T = 2.23 to T = 2.30, with increments of 0.0001, the simulation was run to find the corner and edge tensors. The mag-netization was then calculated as detailed in section 3.3.1. This was done for different bond dimensions: χ = 4, 8, 12 and 24, with a tolerance for the convergence criterion of τ = 10−7. The results are shown in figure 25, along with the exact solution provided by Onsager (1944). The higher the bond dimension χ, the more closely the exact solution is approached. Thus, for higher bond dimensions, the simulation is more accurate. This can be explained by the fact that a finite bond dimension χ imposes a cutoff on the correlation length. Since this effect plays a bigger role in simulations using a smaller bond dimension χ, those simulations produce less accurate results.
2.23 2.24 2.25 2.26 2.27 2.28 2.29 2.30 T 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 m χ=4 χ=8 χ=12 χ=24 exact
Figure 25: The magnetization m as a function of temperature for the square lattice. The magnetization is found using corner and edge tensors with bond dimension χ = 4, 8, 12 and 24. The exact solution is indicated in black.
The same was done for the honeycomb lattice, at bond dimensions χ = 4, 8, 12, 16, 20 and 24. For this lattice, the temperature was ranged from T = 1.49 to T = 1.54, with
increments of 0.0001. The tolerance used is τ = 10−7. In figure 26, the simulation results are shown. The exact value of the critical temperature T 7c = 2
ln(2+√3) ≈ 1.5186514 is indicated
by a dashed line (Yeomans, 1992). This shows that the phase transition occurs closer to the critical temperature as the bond dimension χ is increased. Thus, also for this lattice type the CTMRG method produces more accurate results for higher χ, as expected.
1.49 1.50 1.51 1.52 1.53 1.54 T 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 m χ=4 χ=8 χ=12 χ=16 χ=24 Tc
Figure 26: The magnetization m as a function of temperature for the honeycomb lattice. The magne-tization is found using corner and edge tensors with bond dimension χ = 4, 8, 12, 16, 20 and 24. The critical temperature Tc is indicated by a dashed line.
3.56 3.58 3.60 3.62 3.64 3.66 3.68 3.70 T 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 m χ=4 χ=12 χ=16 χ=20 χ=24 Tc
Figure 27: The magnetization m as a function of temperature for the triangular lattice. The mag-netization is found using corner and edge tensors with bond dimension χ = 4, 12, 16, 20 and 24. The critical temperature Tc is indicated by a dashed line.
Also for the triangular lattice the magnetization as a function of temperature was cal-culated for bond dimensions χ = 4, 12, 16, 20 and 24, with a tolerance of τ = 10−7. The temperature was ranged from T = 3.55 to T = 3.70 and was incremented in steps of 0.0001. The results are shown in figure 27, where the dashed line indicates the exact value of the criti-cal temperature Tc4 = ln(3)4 ≈ 3.6409569 (Yeomans, 1992). Also for this simulation, the phase
transition occurs closer to the critical temperature as the bond dimension χ is increased. The CTMRG method modified for this lattice type produces more accurate results for higher χ, as expected.
For the square lattice, the magnetization was also calculated for different tolerances for con-vergence, as used in equation 27. This was done for bond dimension χ = 4. The tolerance τ was varied from 10−7 to 10−2, increasing with factors of 10. These results are shown in figure 28, along with the exact solution (Onsager, 1944). The smaller the tolerance, the smaller the ‘tail’ of the curve and the better the exact solution is approximated. This agrees with expec-tation based on the convergence criterion as expressed in equation 27, as a smaller tolerance τ implies the corner and edge tensors take longer to converge and thus provide a more accurate approximation of the infinite tensor network.
For all other results, the tolerance is set to τ = 10−7.
2.23 2.24 2.25 2.26 2.27 2.28 2.29 2.30 T 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 m τ=10−7 τ=10−6 τ=10−5 τ=10−4 τ=10−3 τ=10−2 exact
Figure 28: The magnetization m as a function of temperature for the square lattice with bond dimen-sion χ = 4. The tolerance τ for convergence is varied between 10−7 and 10−2, increasing with powers of 10. The exact solution, provided by Onsager (1944), is indicated in black.
4.1.2 Free energy
In figure 29, the simulation results for the free energy of the square lattice as a function of temperature are presented. For temperatures ranging from T = 0.5 to T = 3.5, with increments of 0.01, the simulation was run to find the corner and edge tensors. The free energy was then calculated as detailed in section 3.3.2. This was done for bond dimensions χ = 6, 8, 12, 16 and 24. The resulting curves are shown in figure 29a, along with Onsager’s exact solution, as given in equation 5. As all results seem similar to the exact solution, I have also calculated the relative error δf for each bond dimension χ:
δf =
|f − fexact|
fexact
, (34)
where f is the free energy as obtained from the simulation and fexact is the free energy
according to Onsager’s solution. In figure 29b, the relative error as a function of temperature is shown for the different bond dimensions.
Contrary to the results for the magnetization, the results for the free energy already are quite accurate for low bond dimensions: the curve with χ = 4 has a maximum relative error of the order of 10−5 compared to the exact result. This can be explained by the fact that the free energy is a global variable, whereas the magnetization is a local variable.
For all bond dimensions χ, the relative error shows a maximum close to the critical temperature. This indicates the simulation is least accurate at the phase transition.
0.5 1.0 1.5 2.0 2.5 3.0 3.5 T 2.8 2.6 2.4 2.2 2.0 f χ=4 χ=8 χ=12 χ=16 χ=24 exact
(a) The free energy f
2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 T 10-16 10-15 10-14 10-13 10-12 10-11 10-10 10-9 10-8 10-7 10-6 10-5 10-4 δf χ=4 χ=8 χ=12 χ=16 χ=24 Tc
(b) The relative error δf
Figure 29: The free energy f and its relative error δf as a function of temperature for the square
lattice. In (a), the dashed line indicates Onsager’s exact solution as given in equation 5. In figure (b), the relative error compared to this analytical solution is shown for different bond dimensions.
For the honeycomb lattice, the free energy was calculated for temperatures ranging from T = 0.5 to T = 2.0, incrementing in steps of size 0.01. This was done for χ = 4, 8, 12, 16, 20
and 24, resulting in the curve shown in figure 30a. To illustrate how the free energy changes as a function of bond dimension, the deviation with respect to the free energy for χ = 24 is shown in figure 30b. As one would expect, this deviation gets smaller for increasing bond dimensions. The biggest deviation, for the curve with χ = 4, is still relatively small: it is of the order of 10−5. The critical temperature is indicated by a dashed line; showing the biggest deviations occur at the phase transition.
0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 T 1.80 1.75 1.70 1.65 1.60 1.55 1.50 1.45 1.40 f χ=4 χ=8 χ=12 χ=16 χ=20 χ=24
(a) The free energy f
1.49 1.50 1.51 1.52 1.53 1.54 T 10-15 10-14 10-13 10-12 10-11 10-10 10-9 10-8 10-7 10-6 10-5 10-4 δf χ=4 χ=8 χ=12 χ=16 χ=20 Tc (b) The deviation δf
Figure 30: The free energy f as a function of temperature for the honeycomb lattice. In figure (b), the deviation of the free energy for bond dimensions χ = 4, 8, 12, 16, 20 compared to the free energy for bond dimension χ = 24 is shown. The critical temperature Tc is indicated by a dashed line.
0.5 1.0 1.5 2.0 2.5 3.0 T 1.62 1.60 1.58 1.56 1.54 1.52 1.50 f χ=4 χ=8 χ=12 χ=16 χ=20 χ=24
(a) The free energy f
3.4 3.5 3.6 3.7 3.8 3.9 T 10-14 10-13 10-12 10-11 10-10 10-9 10-8 10-7 10-6 10-5 10-4 10-3 10-2 δf χ=4 χ=8 χ=12 χ=16 χ=20 Tc (b) The deviation δf
Figure 31: The free energy f as a function of temperature for the triangular lattice. In figure (b), the deviation of the free energy for bond dimensions χ = 4, 8, 12, 16, 20 compared to the free energy for bond dimension χ = 24 is shown.The critical temperature Tc is indicated by a dashed line.
For the triangular lattice, the free energy was calculated for temperatures ranging from T = 0.5 to T = 3.0, incrementing again in steps of size 0.01. This was done for χ = 6, 8, 12, 16, 20 and 24. The results are shown in figure 31a. In figure 31b, the deviation with respect to the free energy for χ = 24 is shown. Also for this lattice type, the deviation gets smaller for increasing bond dimensions, as expected. The biggest deviation, for the curve with χ = 4 is of the order of 10−4. The critical temperature is indicated by a dashed line. This shows the deviations are biggest at the phase transition.
4.1.3 Correlation length
In figure 32, the simulation results for the correlation length of the square lattice as a function of temperature are presented. For temperatures ranging from T = 2.23 to T = 2.30, with increments of size 0.0001, the simulation was run to find the corner and edge tensors. The correlation length was then calculated as detailed in section 3.3.3. This was done for χ = 4, 8, 12, 16, 20 and 24, resulting in the curves in figure 32. The exact critical temperature Tc = 2
ln(1+√2) ≈ 2.2691853 is indicated by a dashed line (Yeomans, 1992).
2.23 2.24 2.25 2.26 2.27 2.28 2.29 2.30 T 0 200 400 600 800 1000 ξ χ=4 χ=8 χ=12 χ=16 χ=20 χ=24 Tcrit
Figure 32: The correlation length ξ for the square lattice as a function of temperature for bond dimensions χ = 4, 8, 12, 16 and 24. The exact critical temperature Tc is indicated by the dashed line.
For an infinite system, the correlation length should converge at the critical temperature. This is not the case for the obtained results, as can be seen in the graph. Instead, a maximum can be observed that increases in height for increasing bond dimensions χ. This behavior can be explained by the fact that a finite bond dimension χ imposes a cutoff on the correlation
length χ.
Moreover, figure 32 shows the phase transition (i.e. the maximum of the curve) occurs closer to the critical temperature for increasing bond dimensions χ. This agrees with the results found for the magnetization.
The same was done for the honeycomb lattice, at bond dimensions χ = 4, 8, 12, 16, 20 and 30. For this lattice, the temperature was ranged from T = 1.515 to T = 1.525 and was incre-mented in steps of 0.0001. In figure 33, the simulation results are shown. The exact value of the critical temperature T 7c = 2
ln(2+√3) ≈ 1.5186514 is indicated by a dashed line (Yeomans,
1992). This shows that the phase transition occurs closer to the critical temperature as the bond dimension χ is increased.
The curve for χ = 30 is not completely smooth. This behavior is unexpected and further investigation is needed to explain this.
1.516 1.518 1.520 1.522 1.524 T 0 500 1000 1500 2000 ξ χ=4 χ=8 χ=12 χ=16 χ=20 χ=30 Tc
Figure 33: The correlation length ξ for the honeycomb lattice as a function of temperature for bond dimensions χ = 4, 8, 12, 16 and 30. The critical temperature Tc is indicated by a dashed line.
For the triangular lattice, the correlation length as a function of temperature was cal-culated for bond dimensions χ = 4, 8, 12, 16, 20 and 24. The temperature was ranged from T = 3.60 to T = 3.75, incrementing in steps of 0.0001. The results are shown in figure 34. The dashed line indicates the exact value of the critical temperature Tc4 = ln(3)4 ≈ 3.6409569
(Yeomans, 1992). Also for these simulations, the phase transition occurs closer to the critical temperature as the bond dimension χ is increased, as expected.
3.60 3.62 3.64 3.66 3.68 3.70 3.72 3.74 T 0 200 400 600 800 1000 1200 ξ χ=4 χ=8 χ=12 χ=16 χ=20 χ=24 Tc
Figure 34: The correlation length ξ for the triangular lattice as a function of temperature for bond dimensions χ = 4, 8, 12, 16, 20 and 24. The critical temperature Tc is indicated by a dashed line.
4.2 The critical temperature and critical exponents
Even though the system considered in these simulations is infinite, we can still observe finite size effects. These arise because the bond dimension χ is finite, which results in an effective correlation length ξef fχ at the critical temperature
ξχef f = ξ(t = 0, χ) ∼ χb ∼ ξef f ∗(T = T∗, χ), (35)
playing a similar role as a finite system size L. Exploiting these finite size effects enables us to obtain values for the critical temperature Tc and critical exponents β and ν. In order to
do this, equation 20 can be rewritten as
T∗(χ) = c · χ−b+ Tc or (36)
T∗(χ) = d · ξef f ∗(T = T∗, χ)−1/ν+ Tc, (37)
where c and d are arbitrary constants and T∗(χ) is the temperature at which the phase transition occurs if the system is simulated using a maximum bond dimension χ.
Now, for different bond dimensions χ = 8, 12, 16, 20, 24, 30, 36, 40, 44, 48 and 52, a simula-tion is performed to find T∗(χ). The temperature for which the slope of the magnetization curve is the steepest is used as Tm∗(χ). The error on Tm∗(χ) is determined by the step size used for T , which is 0.0001. The data is then plotted against χ and fitted to equation 36 to find Tc. The results are shown in figure 35, 36 and 37 for the square, honeycomb and triangular
0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 1/χ 2.2690 2.2692 2.2694 2.2696 2.2698 2.2700 2.2702 T ∗ m
Figure 35: Tm∗ as a function of the inverse bond dimension 1/χ for the square lattice. A fit to y = c · x−b+ T
c yields c = 0.93 ± 0.11, b = 2.5 ± 0.5 and Tc = 2.26914 ± 7 · 10−5. The dashed line
indicates the exact value for the critical temperature T c = 2 ln(1+√2) ≈ 2.2691853 (Yeomans, 1992). 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 1/χ 1.5180 1.5185 1.5190 1.5195 1.5200 1.5205 1.5210 1.5215 1.5220 T ∗ m
Figure 36: Tm∗ as a function of the inverse bond dimension 1/χ for the honeycomb lattice. A fit to y = c · x−b+ Tcyields c = 2.434 ± 0.023, b = 3.40 ± 0.04 and Tc= 1.51869 ± 6 · 10−5. The dashed line
Figure 37: Tm∗ as a function of the inverse bond dimension 1/χ for the triangular lattice. The data are fitted to y = c · x−b+ T
c, yielding c = 1.99 ± 0.09, b = 2.57 ± 0.02 and Tc= 3.64096 ± 7 · 10−5. The
dashed line indicates the exact value for the critical temperature Tc4= ln(3)4 ≈ 3.6409569 (Yeomans, 1992).
For the square lattice, a fit to y = c · x−b+ Tc yields c = 0.93 ± 0.11, b = 2.5 ± 0.5 and
Tc= 2.26914 ± 7 · 10−5. For the honeycomb lattice, the results of the fit are c = 2.434 ± 0.023,
b = 3.40 ± 0.04 and Tc = 1.51869 ± 6 · 10−5. For the triangular lattice, the resulting values
are c = 1.99 ± 0.09, b = 2.57 ± 0.02 and Tc= 3.64096 ± 7 · 10−5.
For each bond dimension, T∗(χ) is also determined using the correlation length: the tem-perature at which the curve of ξ shows a maximum is used as Tξ∗(χ). The error on Tξ∗(χ) is determined by the step size used for T , which is 0.0001. Then, these data are also fitted to equation 36 to find Tc. The results are shown in figure 38, 39 and 40 for the square,
honeycomb and triangular lattice, respectively.
For the square lattice, a fit to y = c · x−b+ Tc yields c = 1.01 ± 0.18, b = 2.57 ± 0.04 and
Tc= 2.26917 ± 4 · 10−5. For the honeycomb lattice, the results of the fit are c = 2.78 ± 0.15,
b = 3.450 ± 0.026 and Tc = 1.51868 ± 4 · 10−5. For the triangular lattice, the fit yields
0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 1/χ 2.2690 2.2692 2.2694 2.2696 2.2698 2.2700 2.2702 2.2704 T ∗ ξ
Figure 38: Tξ∗(χ) as a function of the inverse bond dimension 1/χ for the square lattice. The data are fitted to y = c · x−b+ T
c, yielding c = 1.01 ± 0.18, b = 2.57 ± 0.04 and Tc = 2.26917 ± 4 · 10−5.
The dashed line indicates the exact value for the critical temperature T c = 2 ln(1+√2) ≈ 2.2691853 (Yeomans, 1992). 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 1/χ 1.5180 1.5185 1.5190 1.5195 1.5200 1.5205 1.5210 1.5215 1.5220 T ∗ ξ
Figure 39: Tξ∗(χ) as a function of the inverse bond dimension 1/χ for the honeycomb lattice. The data are fitted to y = c · x−b+ Tc, yielding c = 2.78 ± 0.15, b = 3.450 ± 0.026 and Tc= 1.51868 ± 4 · 10−5.
The dashed line indicates the exact value for the critical temperature T 7c = 2
ln(2+√3) ≈ 1.5186514
0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 1/χ 3.6405 3.6410 3.6415 3.6420 3.6425 3.6430 T ∗ ξ
Figure 40: Tξ∗(χ) as a function of the inverse bond dimension 1/χ for the triangular lattice. The data are fitted to y = c · x−b+ T
c, yielding c = 1.48 ± 0.08, b = 2.44 ± 0.03 and Tc= 3.64095 ± 4 · 10−5. The
dashed line indicates the exact value for the critical temperature Tc4= ln(3)4 ≈ 3.6409569 (Yeomans, 1992).
Furthermore, the obtained values for T∗(χ) are plotted against ξef f ∗(T = T∗, χ) and fitted
to equation 37, once again to find Tc. The results for the square lattice are shown in figure
41; a fit to y = d · x−1/ν+ Tcyields d = 0.373 ± 0.018, ν = 0.97 ± 0.04, Tc= 2.26915 ± 4 · 10−5.
For the honeycomb lattice, the fit yields d = 0.1248 ± 0.0013, ν = 1.01 ± 0.07, Tc =
1.51866 ± 3 · 10−5, as shown in figure 42. The results for the triangular lattice, as shown in figure 43, are d = 0.846 ± 0.026, ν = 0.99 ± 0.03, Tc= 3.64089 ± 7 · 10−5.
All three values found for ν agree with the exact value of the critical exponent ν = 1. The fact that ν is the same for each of the three lattice types confirms its universality.
0.000 0.001 0.002 0.003 0.004 0.005 0.006 1/ξeff 2.2690 2.2695 2.2700 2.2705 2.2710 2.2715 T ∗
Figure 41: T∗ as a function of the inverse effective correlation length 1/ξef f for the square lattice. A fit to y = d · x−1/ν+ Tc is done, yielding d = 0.373 ± 0.018, ν = 0.97 ± 0.04, Tc= 2.26915 ± 4 · 10−5.
The dashed line indicates the exact value for the critical temperature T
c = ln(1+2√2) ≈ 2.2691853 (Yeomans, 1992). 0.000 0.001 0.002 0.003 0.004 0.005 0.006 1/ξeff 1.5186 1.5188 1.5190 1.5192 1.5194 T ∗
Figure 42: T∗as a function of the inverse effective correlation length 1/ξef f for the honeycomb lattice.
A fit to y = d · x−1/ν+ Tcis done, yielding d = 0.1248 ± 0.0013, ν = 1.01 ± 0.07, Tc= 1.51866 ± 3 · 10−5.
The dashed line indicates the exact value for the critical temperature T 7c = ln(2+2√3) ≈ 1.5186514
0.0000 0.0005 0.0010 0.0015 0.0020 0.0025 0.0030 0.0035 0.0040 0.0045 1/ξeff 3.641 3.642 3.643 3.644 3.645 T ∗
Figure 43: T∗as a function of the inverse effective correlation length 1/ξef f for the triangular lattice. A fit to y = d·x−1/ν+Tcis done, yielding d = 0.846±0.026, ν = 0.99±0.03, Tc = 3.64089±7·10−5. The
dashed line indicates the exact value for the critical temperature T4
c = ln(3)4 ≈ 3.6409569 (Yeomans,
1992).
The results for the critical temperatures for the different lattice types are summarized in table 1. All results agree with the exact critical temperatures of Tc= ln(1+2√
2) ≈ 2.2691853,
T 7c = 2
ln(2+√3) ≈ 1.5186514 and T 4
c = ln(3)4 ≈ 3.6409569 for the square, honeycomb and
triangular lattice, respectively (Yeomans, 1992).
Table 1: An overview of the critical temperatures Tcfor the square, honeycomb and triangular lattice,
found by fits based on finite size scaling.
fit on Tc T 7c Tc4
Tm∗(χ) 2.26914 ± 7 · 10−5 1.51869 ± 6 · 10−5 3.64096 ± 7 · 10−5 Tξ∗(χ) 2.26917 ± 4 · 10−5 1.51868 ± 4 · 10−5 3.64095 ± 4 · 10−5 T∗(ξef f) 2.26915 ± 4 · 10−5 1.51866 ± 3 · 10−5 3.64089 ± 7 · 10−5
With relation 35, equation 24 can be rewritten, giving: ln[m(t = 0, χ)] = −β
νln(ξ
ef f
χ ) + C. (38)
For different bond dimensions χ = 8, 12, 16, 20, 24, 30, 36, 40, 44, 48, 52, 60 and 80, a simulation is performed to find the effective correlation length ξχef f and the magnetization at the critical
temperature m(t = 0, χ). Then, the logarithm of these values is calculated. The data are fitted to equation 38 in order to obtain the ratio β/ν.
The results for the square lattice are shown in figure 44. A fit of the data to y = −βνx + C yields β/ν = 0.1249 ± 0.0003 and C = −0.2544 ± 0.0021. The results for the honeycomb lattice, shown in figure 45, are β/ν = 0.1249 ± 0.0004 and C = −0.2308 ± 0.0028. For the triangular lattice, the results are β/ν = 0.1254 ± 0.0005, C = −0.1778 ± 0.0023, as shown in figure 46.
For each of the three lattice types, the results agree with the exact value of the ratio of critical exponents β/ν = 1/8 = 0.125. The results confirm the universality of critical exponents, since β/ν is the same for the different lattice types.
4 5 6 7 8 ln[ξeff(χ)] 1.4 1.3 1.2 1.1 1.0 0.9 0.8 0.7 0.6 ln [m ( χ )]
Figure 44: Finite size scaling of the square lattice. For different bond dimensions, ln[m(t = 0, χ)] is plotted against ln(ξef f
χ ). The data are then fitted to y = − β
νx + C, yielding β/ν = 0.1249 ± 0.0003,
4 5 6 7 8 ln[ξeff(χ)] 1.4 1.3 1.2 1.1 1.0 0.9 0.8 0.7 0.6 ln [m ( χ )]
Figure 45: Finite size scaling of the honeycomb lattice. For different bond dimensions, ln[m(t = 0, χ)] is plotted against ln(ξef f
χ ). The data are then fitted to y = − β νx + C, yielding β/ν = 0.1249 ± 0.0004, C = −0.2308 ± 0.0028. 4 5 6 7 8 ln[ξeff(χ)] 1.4 1.3 1.2 1.1 1.0 0.9 0.8 0.7 0.6 ln [m ( χ )]
Figure 46: Finite size scaling of the triangular lattice. For different bond dimensions, ln[m(t = 0, χ)] is plotted against ln(ξef f
χ ). The data are then fitted to y = − β
νx + C. This yields β/ν = 0.1254 ±
0.0005, C = −0.1778 ± 0.0023.
As a consistency check, equation 26 is used to perform a data collapse. Using the exact values β = 1/8 and ν = 1, y = m(t, χ) · (ξχef f)β/ν is plotted against x = t · (ξχef f)1/ν for
of the data collapses for the square, honeycomb and triangular lattice, respectively. For each lattice type, all points fall on a single curve F , as expected. This validates the results of the simulation. 3 2 1 0 1 2 3 t·L1/ν 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 m ( t, L) · L β/ ν χ=12 χ=16 χ=20 χ=24
Figure 47: Data collapse of the square lattice. Using the exact values β = 1/8 and ν = 1, all points fall on a single curve F , as expected.
1.5 1.0 0.5 0.0 0.5 1.0 1.5 t·L1/ν 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 m ( t, L) · L β/ ν χ=12 χ=16 χ=20 χ=24
Figure 48: Data collapse of the honeycomb lattice. Using the exact values β = 1/8 and ν = 1, all points fall on a single curve F , as expected.
3 2 1 0 1 2 3 t·L1/ν 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 m ( t, L) · L β/ ν χ=12 χ=16 χ=20 χ=24
Figure 49: Data collapse of the triangular lattice. Using the exact values β = 1/8 and ν = 1, the points around t = 0 fall on a single curve F , as expected.
5
Discussion and conclusion
The original Corner Transfer Matrix Renormalization Group method, as introduced by Nishino & Okunishi (1995), provides a useful approach to simulate a system on a square lattice. To confirm the validity of this original method, which was already shown in previous studies, I have used the CTMRG to simulate the partition function of the two-dimensional Ising model on a regular square lattice. Then, I have made generalizations of the CTMRG method for the honeycomb and triangular lattices, which was not done before.
To study the validity and accuracy of the generalized CTMRG method for the different lattice types, three thermodynamic quantities have been calculated: the magnetization, free energy and correlation length. All results show the expected behavior. Furthermore, values for the critical exponents ν and the ratio β/ν have been determined, as well as the critical temperature Tc of each of these lattices.
Table 2: An overview of the results for the critical exponents and their exact values.
ν β/ν
0.97 ± 0.04 0.1249 ± 0.0003 7 1.01 ± 0.07 0.1249 ± 0.0004 4 0.99 ± 0.03 0.1254 ± 0.0005
exact 1 0.125
As is shown in table 2, all obtained values for ν and β/ν agree with their exact values. In table 3, it can be seen that the results for Tc obtained through a fit on Tξ∗(χ) agree with the
exact values. Also the values for Tc obtained through fits on Tm∗(χ) and T∗(ξeff) agree with
the exact values. Therefore, I can conclude that the generalizations of the CTMRG method to the honeycomb and triangular lattice are valid. The generalization does not change the leading computational cost of the CTMRG algorithm; this is O(χ3) for all three lattices. Table 3: An overview of the results for the critical temperatures found by a fit on Tξ∗(χ) and the exact results.
Tc T 7c Tc4
fit on Tξ∗(χ) 2.26917 ± 4 · 10−5 1.51868 ± 4 · 10−5 3.64095 ± 4 · 10−5
exact result 2.2691853 1.5286514 3.6409569
are of the same order of accuracy. In general, the accuracy of the results can be improved by doing simulations for bigger maximum bond dimensions χ and using smaller step sizes when varying the temperature. However, this requires more running time, more computer power or a more efficient code. To avoid this, instead of computing the magnetization over a fine grid in an extended region, a function which locates a maximum of the slope could be used. This would require fewer evaluations than computing m on many points.
Also, a different algorithm, like the Tensor Network Renormalization approach (Evenbly & Vidal, 2015), might help to increase efficiency. Comparing the efficiencies and running times of different algorithms might be an interesting topic to study in further research. Also, one could simulate the partition function of other models, such as the Potts model, on a honeycomb or triangular lattice. Moreover, since the d-dimensional quantum Ising model can be mapped to the d + 1-dimensional classical Ising model (Barma & Shastry, 1978), the one dimensional quantum Ising model might be an interesting topic to study next.
Acknowledgements
I would like to thank dr. Philippe Corboz and Schelto Crone for their supervision and feedback during this project. I would also like to thank Edan Lerner for acting as my second assessor.
References
Barma, M., & Shastry, B.M. 1978, Phys. Rev. B., 18, 3351 Baxter, R.J. 1978, J. Stat. Phys. Appl., 19, 461
Boer, A. 2011, J. Physica A Stat. Mech. Appl., 390, 4203 Chan, Y. 2013, J. Phys. A: Math. Theor., 46, 125009
Evenbly, G., & Vidal, G. 2015, Phys. Rev. Lett., 115, 180405 Fisher, M.E., & Barber, M.N. 1972, Phys. Rev. Lett., 28, 1517 Ising, E. 1925, Z. Phys., 31, 253
Kramers, H.A., & Wannier, G.H. 1941, Phys. Rev., 60(3), 252 Lenz, W. 1920, Physikalische Zeitschrift, 21, 613
Levin, M., & Nave, C.P. 2007, Phys. Rev. Lett., 28, 1517 Nishino, T., & Okunishi, K. 1995, J. Phys. Soc. Jpn., 65, 891 Onsager, L. 1944, Phys. Rev., Series II, 65 (3–4), 117
Or´us, R. 2012, Phys. Rev. B85, 205117
Potts, R.B. 1952, Math. Proc. Camb. Philos. Soc., 48, 106 Stanley, H.E. 1999, Rev. Mod. Phys., 71, S358
Taroni, A. 2015, Nat. Phys., 11, 997
Verstraete, F., Cirac, J.I. & Murg, V. 2008, Adv. Phys. 57, 143 White, S.R. 1992, Phys. Rev. Lett. 69(19), 2863
Yang, C.N. 1952, Phys. Rev., 85 (5), 808
Yeomans, J.M. 1992, Statistical Mechanics of Phase Transitions. (Oxford Science Publica-tions, Oxford), 117
