MSc Physics
Theoretical Physics
Master Thesis
Bethe Ansatz Equations for the
Square-Triangle Random Tiling Model
by
Kevin A.M.B. Smit
6037844
August 8, 2014
60 ECTS
November 2013 - August 2014
Supervisor:
Prof. dr. B. Nienhuis
Examiner:
Dr. V. Gritsev
Contents
1 Introduction 5
1.1 Quasicrystals . . . 6
1.1.1 Discovery of quasicrystals . . . 6
1.1.2 Properties of quasicrystals . . . 7
1.1.3 Long range-order and quasiperiodicity . . . 7
1.1.4 Forbidden crystalline rotational symmetries . . . 8
1.1.5 Quasiperiodicity as described by a restricted periodic func-tion in hyperspace . . . 8
1.1.6 Stabilization of quasicrystals: energy versus entropy . . . 9
1.2 Random tiling models . . . 10
1.2.1 Introduction . . . 10
1.2.2 Random tiling models, as described in hyperspace . . . . 10
1.2.3 Random tiling hypotheses . . . 12
1.3 Partition function for random tiling models . . . 12
1.4 Transfer matrix method . . . 13
1.5 Coordinate Bethe Ansatz . . . 14
2 Construction of the statistical lattice model 16 2.1 Square-triangle random tiling model . . . 16
2.2 Transformation of the square-triangle model into a triangular lat-tice model . . . 18
2.3 Construction of the transfer matrix for the lattice model . . . . 19
3 Coordinate Bethe Ansatz 22 3.1 Preliminaries to the Coordinate Bethe Ansatz . . . 22
3.2 Derivation of the amplitude relations . . . 22
3.3 Derviation of the Coordinate Bethe Ansatz equations . . . 27
3.4 Logarithmic Coordinate Bethe Ansatz equations . . . 29
4 Numerical analyses 31 4.1 Application of the Newton-Raphson method to the Bethe Ansatz equations and numerical analyses of the transfer matrix . . . 31
4.2 Parameter space of the system . . . 32
4.3 Results of the numerical analyses on the Bethe Ansatz equations and the transfer matrix . . . 33
5 Analytic solution to the Bethe Ansatz equations 36 5.1 Bethe Ansatz equations . . . 36 5.2 Continuum limit of the system . . . 37 5.2.1 Density of roots on the solution curves . . . 37 5.3 Application of perturbation theory on the Bethe Ansatz equations 38 5.3.1 Preliminaries of perturbation theory . . . 38 5.3.2 Zeroth order in perturbation theory . . . 39 5.3.3 First order in perturbation theory . . . 41 5.3.4 First order in perturbation theory in expansion parameterw 44
6 Discussion 46
7 Conclusion 47
A Long-range order 48
B Geometrical constraint 49
C Application of the transfer matrix to the diamond tiling 50
D Perron-Frobenius theorem 53
D.1 Theorem for regular matrices . . . 53 D.2 Theorem for nonnegative matrices . . . 53 D.3 List of definitions . . . 54
E Newton-Raphson method 55
E.1 Single variable method . . . 55 E.2 Multivariate method . . . 56 F Mathematica 8 program for the transfer matrix 57 G Mathematica 8 program for the Newton-Raphson method 72
Abstract
The hexagonal phase of the square-triangle random tiling models is inves-tigated. We transform this tiling to a triangular lattice model, for which we construct a transfer matrix. The Coordinate Bethe Ansatz is used to diagonal-ize the transfer matrix, thereby obtaining the Bethe Ansatz equations. We are concerned with the solutions to the Bethe Ansatz equations, corresponding to the largest eigenvalue of the transfer matrix. Numerical solutions to the Bethe Ansatz equations are obtained and studied. We investigate the continuum limit of the system with perturbation theory, thereby attempting to obtain a pertur-bation series expansion for the begin- and endpoints of the 2 solution curves. We are able to obtain the zeroth order terms for the begin- and endpoints in the thermodynamic limit. However we are not able to obtain the higher order terms in neither the thermodynamic or the continuum limit. Our results con-tradict the results of B. Nienhuis, who was able to obtain a perturbation series expansion for the begin- and endpoints, in his model for the hexagonal phase of the square-triangle random tiling model.
Acknowledgments
I want to thank Bernard Nienhuis for his outstanding job as my supervisor. He always took the time to explain the subjects of this thesis to me and to guide me through this project. He even took the time to give me feedback on my thesis while he was on his holiday, for which I want to thank him greatly.
I also want to thank Vladimir Gritsev for his interest in my work.
Last but not least I want to thank my family, especially my parents, for their ongoing love and support. Which really helped me through the tough last months of my work on my thesis. Without there support I would not be able to finish my master thesis in time.
Chapter 1
Introduction
In this thesis we will examine the square-triangle random tiling model. For this random tiling model M. Widom and P.A. Kalugin have made important contri-butions. In 1992 Widom [1] discovered the Coordinate Bethe Ansatz equations and a year later, Kalugin [2] used a different lattice representation to obtain the same Bethe Ansatz equations. However, Kalugin was able to derive the analytic expression for the entropy density, in a part of the phase diagram and in the thermodynamic limit.
In this first chapter we present a solid introduction in quasicrystals and ran-dom tiling models, which are used to describe some (quasi-)crystals. Further-more we will introduce some important techniques such as the transfer matrix method and the Coordinate Bethe Ansatz. Although a good understanding of quasicrystals is not necessary for our study of random tiling models, we include this information as a useful overview for future students, who want to do re-search in these fields. After having introduced the basic concepts and methods, we will continue, in chapter 2, to transform our square-triangle random tiling model into a lattice model and apply the transfer matrix method to this model. Chapter 3 concerns the diagonalization of the transfer matrix by means of the Coordinate Bethe Ansatz. In this chapter we will obtain the important Bethe Ansatz equations. In chapter 4 we will study the Bethe Ansatz equation nu-merically, by means of the Newton-Raphson method, and we will discover some useful properties of the numerical solutions to the Bethe Ansatz equations. Then in chapter 5 we will study the Bethe Ansatz equations with perturbation theory in order to obtain the begin- and endpoints of the solution curves of the roots to the Bethe Ansatz equations. After obtaining the begin- and endpoints of the solution curves, one can calculate the largest eigenvalue of the system. From the largest eigenvalue of the system, one can calculate the entropy, the free energy and the partition function of the system.
The objective of this thesis is to obtain the perturbation series expansions for the begin- and endpoints of the solution curves and to investigate whether these are summable. My supervisor B. Nienhuis, already obtained such expansions from a lattice model and requested my to confirm his result by studying a slightly different lattice model, for the square-triangle random tiling.
1.1
Quasicrystals
1.1.1
Discovery of quasicrystals
Quasicrystals were discovered in 1984 by Shechtman and his co-workers, in their research on the metallurgical properties of Al-Fe and Al-Mn alloys [6]. They discovered diffraction patterns of icosahedral rotational symmetry, consisting out of sharp Bragg peaks, in the metastable phase of an Al-Mn alloy. The research group obtained this Al-Mn alloy in its metastable icosahedral phase by rapidly quenching the sample from the melt [14]. Up until that point in time, the scientific community did not believe that 5-fold, or any other non-crystalline, rotational symmetry in two dimensions could exist in materials with long-range order and sharp Bragg peaks [4]. They were convinced that only crystalline rotational symmetries could exist in these types of materials, i.e. 1-,2-,3-,4- and 6-fold rotational symmetries [3].
By now many compounds with quasicrystalline phases have been observed. Even compounds exhibiting other non-crystalline rotational symmetries, such as 8-, 10- and 12-fold, have been discovered [5].
In the first couple of years after the discovery of quasicrystals, there were some speculations that quasicrystals might be inherently disordered and un-stable. One of the reasons, for these speculations, was that the known qua-sicrystalline materials were only weakly metastable against a transformation to a lower-energy conventional crystalline state. Fortunately, at least a dozen thermodynamically stable compounds have been discovered since then. Thereby demonstrating this argument for speculations to be false [5] and reinforcing the understanding that quasicrystals are a new class of solids.
An important and frequently used method to investigate the structure of condensed matter systems and to discover their associated properties, is the use of diffraction experiments. These experiments produce an intensity distribution related to the Fourier transform of the underlying structure, for the material un-der investigation. From diffraction experiments, one can deduce properties such as rotational symmetry and long-range order [4]. An example of a decagonal diffraction pattern is given in figure 1.1.
Figure 1.1: Electron diffraction pattern of AlCuCoSi alloy exhibiting a decago-nal rotatiodecago-nal symmetry, figure adopted from [4].
1.1.2
Properties of quasicrystals
Quasicrystals are neither crystals, which are periodically ordered and have a discrete rotational symmetry [5], nor amorphous solids, which are disordered and do not posses any discrete rotational symmetry whatsoever [6]. They are a new class of solids on their own. Quasicrystals are characterized by exhibiting a non-crystalline rotational symmetry, quasiperiodicity, long-range order and a minimum separation distance between the atoms [4].
Their non-crystalline rotational symmetry, is a discrete point group symme-try which is incompatible with periodic translational order. As a consequence, quasicrystals (or materials in a quasicrystalline phase) do not posses crystalline rotational symmetry. Furthermore the exhibited quasiperiodicity is a different kind of translational order [5]. Since the diffraction patterns of quasicrystals consist out of a finite number of Fourier components with sharp diffraction peaks, quasicrystals posses long-range order [4].
1.1.3
Long range-order and quasiperiodicity
As an example of long range-order and quasiperiodicity, we examine the Fi-bonacci chain, which is a sequence consisting out of short (S) and large (L) segments. We start the sequence with a large segment L and apply the iterative rule S→L and L→LS, to obtain successive sequences with increasing lengths. The Fibonacci chain has the property that the occurrence of a short or large seg-ment, at any point in the chain, is uniquely determined by the starting sequence of the chain. Moreover, the sequence has no repetition length (and hence is not periodic), if the ratioL/S is an irrational number. For the canonical Fibonacci chain this ratio is equal to the golden meanτ , which is an irrational number.
The demonstration that the Fibonacci chain exhibits perfect long-range or-der is displayed by the fact that one can calculate the distance D from the origin of the first segment of the chain to the end of thenth segment through
D = S n + E τ n + 1 τ (1.1) In whichE[n+1τ ] is the integer part in (n+1)/τ [4]. For a more elaborate example on why systems can exhibit long-range order without possessing translational periodicity, we refer to appendix A.
Due to the long-range order, the diffraction pattern of the Fibonacci chain consist out of a set of sharp Bragg peaks which fill the reciprocal space in a dense manner. The first step in obtaining the diffraction pattern for the Fibonacci chain is to renormalize S → 1 and L → τ in (1.1), resulting in
xn=n + E τ n + 1 τ (1.2) Where xn is the nth atomic position. This equation consist out of two terms
which correspond to periodic spacings but with incommensurate periods, as in (A.3). Now by performing a Fourier transformation on (1.2), one concludes that the non-zero Fourier components are obtained at the vectors
Qhh0= 2πτ 2 τ2+ 1 h +h 0 τ (1.3)
Whereh and h0 are integers. And the Fourier transform is given by
F (Q) = Fh,h0
X
h,h0
δ(Q − Qhh0) (1.4)
WhereFh,h0 is some complicated amplitude, depending onh and h0. Note that
the one-dimensional structure of diffraction peaks is indexed by two integers, this feature will become more evident in 1.1.5. Furthermore, the sharp peaks of the diffraction pattern fill the reciprocal space densely.
All long-range ordered structures, which are not periodic and posses such a similar singularly continuous Fourier spectrum, posses quasiperiodic transla-tional order [4].
1.1.4
Forbidden crystalline rotational symmetries
Crystalline rotational symmetries (in a two-dimensional space) are restricted to 1-,2-,3-,4- and 6-fold rotations. These rotational symmetries are implied by the existence of a shortest basis vector as demanded by translational symmetry. This requirement restricts the number of possible crystalline rotational symmetries, to 5. The proof for this statement is constructed as follows: denote the shortest basis vector, of a periodic set of points in two dimensions, by r. Let every rotation R(θ), over an angle of θ = 2πk/n with k ∈ {1, ..., n − 1}, keep the periodic set of points invariant, if the lattice structure (formed by these points) has an-fold rotational symmetry. As a consequence r0 =R(θ)r and every linear
combination of the vectors r and r0 represents a point on this lattice structure. These considerations lead to the following set of inequalities
| r0− r |≥| r | and | r0+ r |≥| r | Implying
2 | cos(θ) |≤ 1 (1.5) For all the values of k ∈ {1, ..., n − 1}. From equation (1.5), we conclude that the only allowed values forn, hence the rotational symmetries compatible with translation periodicity, are 1,2,3,4 and 6 [8].
Another more intuitive proof is obtained by means of tiling a plane with regular polygons, in such a way that there are no gaps or overlaps between the tiles. In order to tile the plane without gaps or overlaps, the regular polygons should have vertex angles equal to integral fractions of 2π. The only regular polygons obeying this restriction are the triangle, square, rectangle and hexagon. Therefore, according to this proof, the allowed n-fold rotations are 3,4 and 6. The other two n-fold rotations are obviously compatible with translational periodicity [4].
1.1.5
Quasiperiodicity as described by a restricted
peri-odic function in hyperspace
A very useful method for describing quasicrystals is by making use of the hy-perspace concept. This concept means that any quasiperiodic function can be derived from a periodic function in a higher D-dimensional space (the hyper-space). This hyperspace-dimension is the number of integer independent basis
vectors necessary to span the reciprocal lattice. The quasiperiodic function in the realdk-dimensional space, is then obtained by slicing the higher dimensional
periodic function with ad⊥ = (D − dk)-dimensional subspace (the perp space).
This subspace is incommensurate with the D-dimensional lattice., i.e. it is not parallel to any of the D-dimensional reciprocal lattice planes [5].
As an example of describing quasiperiodicity with a restricted periodic func-tion in a higher dimensional space, we examine the Fourier transform of the density of the solid. Which for an ordinary crystal can be written as a Fourier series ρ(r) = 1 V X k ρ(k)eik·r (1.6)
Where V is the volume of the 3-dimensional crystal and the wave vectors k define a discrete reciprocal lattice. Each wave vector in the sum can be written as an integer linear combination of the three basis vectors vi, which span the
3-dimensional reciprocal lattice, as follows:
k =n1v1+n2v2+n3v3 (1.7)
The basis vectors are integer linearly independent and are three dimensional. Although, in a quasicrystal, the transform is also a Fourier series and the wave vectors also form a discrete reciprocal lattice, the number of integer linearly independent basis vectors required to span the reciprocal lattice exceeds the spatial dimension. As an example one might consider a 3-dimensional qua-sicrystal with an icosahedral symmetry, for which 6 basis vector are needed to span the reciprocal lattice [5].
A slightly different approach is obtained by defining the basis vectors in (1.6) to be higher dimensional independent basis vectors, which are independent over all real numbers. For the icosahedral case we would define 6-dimensional basis vectors pi = vi+ bi, such that bi are orthogonal to any position vector
r in the 3-dimensional subspace. Consequently, one can replace vi by pi in
(1.6), while restricting r to the lower dimensional subspace (three in this case). Therefore one can conclude that a quasiperiodic function is a restriction of a higher dimensional periodic function [14].
1.1.6
Stabilization of quasicrystals: energy versus entropy
There are two possible mechanism for the stabilization of quasicrystals: ener-getic and entropic stabilization. The two main models for the stabilization of quasicrystals are the perfect tiling models, in which energetic effects dominate, and the random tiling models, in which the entropy is dominant. First a model, based on energetic stabilization was proposed by Levine and Steinhardt. Ac-cording to their model “ a structure which possess perfect quasiperiodic trans-lational order, constitutes the minimum-energy configuration of the system.” Another model based on entropic stabilization was proposed by Widom and co-workers. They thought that, in their random tiling model, entropic effects could be sufficient to stabilize the quasicrystalline state, even without favorable energetics [5]. Consider the case in which the quasicrystal, at sufficiently high temperatures, is in its thermal equilibrium state. Then, the entropic contri-bution to the free energy would overcome the energy contricontri-bution. And, as a
consequence, one could conclude that configurational entropy is intrinsic to the quasicrystalline state, in the random tiling model [5].
Both types of models have their disadvantages. An obstacle for the perfect tiling models, such as the Penrose tilings, is that the matching rules can only be enforced by unnaturally fine-tuned atomic interactions [5]. Matching rules impose constraints on the global structure of tilings and may be regarded as schematically representing the energetics of some atomic decoration of tiles, such that all tilings which satisfy the rules are degenerate minimum-energy states [5]. This scenario seems very unlikely, since quasicrystals are common among a large set of metal alloys [5]. Moreover the matching rules are only likely to be applicable to real quasicrystals, if the energetics of the short-range interactions have a significant effect in determining the long-range structure [5]. The random tiling has its obstacles as well. The model supposes a near-degeneracy in energy over a vast and various set of tile clusters. As a conse-quence the system is free to fluctuate from one arrangement to another. There-fore the energy needs to be almost constant over a wide range of atomic clusters. Even though, the model assumes a strong interaction between the atoms of a cluster, in order to prevent the cluster from breaking up into smaller pieces. Besides these strong interactions between the atoms of a cluster, there need to be strong interactions to enforce the tiles tot join edge-to-edge without any holes [5].
Fortunately both the random tiling models as the matching rule models lead to the same reciprocal lattice of Bragg peaks, which is the most important factor for many important physical properties such as phonon properties [5].
1.2
Random tiling models
1.2.1
Introduction
A random tiling model is a model for the tiling of the plane by two or more different types of rigid geometrical units. These geometrical units are all ori-entated differently and these orientations are related by symmetry [5]. Packing rules, such as the constraint that faces should adjoin or Penrose matching rules, are used to determine which tilings are allowed. The packing rules are referred to as matching rules, if they force quasiperiodicity [5]. All the possible tilings are assumed to have an almost degenerate energy [5]. And the number of such tilings increases as exp(S0N ) for some S0, in which N is the total number of tiles
[5]. Furthermore the system under consideration, tries to maximize its entropy by creating the structure with the highest allowed (non-)crystallographic rota-tional symmetry. The fluctuations around the mean rotarota-tional symmetry are described by a gradient-squared elasticity, which sustains the long-range order [5].
1.2.2
Random tiling models, as described in hyperspace
Similar to the discussion of section 1.1.5, we can describe random tiling mod-els by making use of the hyperspace concept. This picture will enable us to introduce some useful concepts.
In our real space we can represent every vertex position xk, of the tiling, as a sum of basis vectors with integer coefficients, given by
xk =
D
X
α=1
nαekα (1.8)
In which xkand the basis vectors ekα are of dimensiondk. Although this
struc-ture can be quasiperiodic, it does not need to be so. To each vertex one can assign integer coordinates nα, which can be interpreted as the lattice points of
the hyperspace. Then (1.8) can be interpreted as a projection from the hyper-space to the real hyper-space, hence a projection from D → d dimensions.
We choose the basis vectors ekα such that (1.8) is unique and such that the
components nα are independent in hyperspace. In order to obtain all
well-defined discrete perp variables, we should choose the maximum value ofD (also known as the lifting dimension) for which uniqueness is possible. However the basis vectors ekαneed not be linear independent over the integers.
We construct a linear independent basis for the hyperspace, by defining a complementary basis e⊥α, such that eα= e⊥α⊕ e
k
α. As a result we define for each
vertex the perp space coordinate x⊥ by x⊥=
D
X
α=1
nαe⊥α (1.9)
These basis vectors span the perp space of dimensiond⊥ =D − dk.
An important notion in the hyperspace concept, is that of the representative surface. As a consequence of the continuous embedding of the physical space into the hyperspace, the tiling raises to a faceted hypersurface which is composed out of a set ofd-faces of the D-space unit cells. This set of d-faces is known as the representative surface. Note that the image of the representative surface, regarded as an undulating cut through a periodic structure (used in section 1.1.5), concurs only in the coarse-grained limit with the image of the faceted hypersurface.
The basis vectors eαdetermine the orientation of the physical space,
embed-ded in the hyperspace. Consequently, the tiling possesses a (non-)crystallographic rotational symmetry, if this orientation is (ir)rational. The phason space, which is a linear subspace of the perp space, is responsible for the irrational orientation [5].
In (1.8) we chose the maximal value forD for which this equation is unique. However, instead of choosingD maximal, we will demand the basis vectors ekαto
be linearly independent over the integers. As a consequence, the hyperspace has a smaller dimension, referred to as the indexing dimension,Dindex =dk+d⊥ph.
Here the perpendicular part of this restricted space is the phason space with dimensiond⊥
ph.
In this thesis we will be examining a random tiling model for which we de-mand the basis vectors ekαto be linearly independent over the integers.
There-fore we will be working with a phason space in stead of a perp space, hence each perp space coordinate x⊥ is replaced by a phason space coordinate x⊥ph. Furthermore, the random tiling model we will be examining has an indexing dimension of 4 and a phason space dimension of 2.
An important quantity in the hyperspace representation of the random tiling model, is the phason strain, which is defined by
E ≡ ∂x
⊥ ph
∂xk (1.10)
This is a (d⊥
ph× dk)-dimensional tensor which is used to measure the statistical
inequality between densities of by symmetry related tiles in different orienta-tions. The state of zero phason strain is the quasicrystalline state. In this state all the by symmetry related orientations of a given polygon are equal in number.
1.2.3
Random tiling hypotheses
Now that we defined the phason strain, we can elaborate on the important random tiling hypotheses.
First we define the entropy, as a function of phason strain, to be σ(E) ≡ lim
l→∞S(E; L)/L d
(1.11) In whichS(E; L) is the entropy of the sub-ensemble of tilings, contained in a box of volumeLd. And the phason strain, averaged over the whole system, in
this volume is denoted by E. Note that the entropy density cannot be dependent on the phason coordinate, since the displacement of the phason coordinate is a macroscopic symmetry in a quasicrystal.
Then, the first random tiling hypothesis states that the entropy densityσ(E) is maximized at E = 0. The second random tiling hypothesis, which does not apply to the tilings we will be considering, applies to maximally random tilings. In these tilings, the weight of each tile is set to equal one. Then the second random tiling hypothesis states that entropy density is quadratic in E near E = 0, hence
σ(E) − σ(0) = −1 2
X
KiljmEilEjm+O(| E |n) (1.12)
In which n > 2 and Kiljm are the elastic constants [5].
1.3
Partition function for random tiling models
The physics of our random tiling model is encapsulated in the partition function for the grand canonical ensemble. In this ensemble the set of our random tiling systems is allowed to exchange energy and tiles with a much larger reservoir. This reservoir is in thermal and chemical equilibrium with the set of random tiling systems [13]. For our model, the partition function is given by
ZA= X tilings exp m X k=1 nkµk ! (1.13) Where the summation is over all the possible tilings, obeying the constraint that the total area A is constant, and we have setkB=T ≡ 1. One obtains the total
area A =Pm
individual areas of the m tiles. A chemical potential µk is assigned to each tile
k.
We choose to examine random tilings for which the internal energy of the system equals zero. As a consequence, there is only an entropic contribution to the free energy. Our partition function is related to the grand potential (or grand free energy [13]), from which we can derive our desired thermodynamic quantities, through
Ω(µk) = − log(ZA) (1.14)
For example, for a given tilel we can calculate its average number by hnli = −
∂Ω(µk)
∂µl
Since we have set the internal energy of the system to zero, the grand po-tential is reduced to −Ω = S(h{nk}i) + m X k=1 hnkiµk =max {nk} S({nk}) + m X k=1 nkµk ! (1.15) From which we can obtain the entropyS and where the maximization is done under the constraint A =Pm
k=1nkak.
Therefore in order to calculate the entropy, and other thermodynamic quan-tities, of a random tiling, one has to calculate the partition function [8].
1.4
Transfer matrix method
In order to calculate the partition function we will be making use of the trans-fer matrix method. This method is one of the standard techniques of lattice statistical mechanics. Here we will give a short example of the transfer matrix method, for a thorough explanation see appendix C. Consider a square lattice, which is taken to be periodic in both the horizontal and the vertical direction. Denote the states of two consecutive rows of horizontal edges by |Ψi and |Ψ0i. Then the transfer matrix element TΨ,Ψ0, for the transition from the state |Ψi
to |Ψ0i, equals the product of the weights of the elementary lattice faces that can construct an intervening layer and equals zero if no such layer can be con-structed [14]. Due to the periodicity in both directions, the partition function reduces to the trace of the transfer matrix, raised to power of the number of vertical layersM , and is given by
Z = Tr TM (1.16)
In the thermodynamic limit of M → ∞ the partition function is equal to the largest eigenvalue of the transfer matrix, raised to M : limM →∞Z = ΛMmax.
We conclude that, in the thermodynamic limit, the calculation of the partition function is reduced to an eigenvalue problem. Hence by calculating the largest eigenvalue of the transfer matrix we can evaluate the partition function in the thermodynamic limit [7].
1.5
Coordinate Bethe Ansatz
In the process of transforming the random tiling model to a lattice model, it is common that one can define conserved quantities. These quantities are conserved in the sense that they are invariant under the action of the transfer matrix on the state they are in. The existence of conserved quantities is due to the construction of the transfer matrix and the periodic boundary conditions. By assigning the horizontal direction of the lattice to denote the position of these conserved quantities and the vertical direction as a time axis, one can interpret the transfer matrix as a time evolution operator and, consequently, the conserved quantities as “particles.” Similar to ordinary particles, these “particles” are or are not able to interact with each other, depending on the construction of the transfer matrix.
Due to the existence of these particles the transfer matrix splits up into a series of diagonal blocks, which relate states with an equal number of particles to each other. Therefore the transfer matrix adopts a block diagonal form and as a consequence the diagonalization of the transfer matrix reduces to the diagonalization of each block, separately.
Consider a block withn particles and let |x > denote the general state of a row of the lattice. For the diamond tiling the row would consist out of L horizontal edges, n would denote the number of occupied edges and the vector x = (x1, ..., xn) would denote the position of each occupied edge. Now in order
to diagonalize the transfer matrix we consider, for each block separately, the following eigenvalue equation:
T X
x0
ψ(x0)|x0>= ΛX
x0
ψ(x0)|x0> (1.17) In which the transfer matrix T only acts on the orthonormal base |x0> and not on the coefficientsψ(x0). Note that the eigenvector of the eigenvalue Λ is a
linear combination of the orthonormal base.
In order to select a specific coefficientψ(x), we perform a scalar product of (1.17) with the state |x >. We then obtain the useful relationship, which will be referred to as the eigenvalue equation from now on
X
x0
< x|T |x0> ψ(x0) = Λψ(x) (1.18) The left hand side of this equation consists out of the sum of the old co-efficients ψ(x0) weighted by the transfer matrix element < x|T |x0 >, for the
transition from state |x0> to |x >. And the right hand side consists out of the new coefficient and the corresponding eigenvalue. Therefore this equation can be interpreted as an evolution equation for the coefficients [7].
At this point in the calculation one assumes that the coefficients can be expressed as linear combinations of products of plane waves. This is known as the Coordinate Bethe Ansatz. Application of this Ansatz to the eigenvalue equation (1.18) results in consistency equations on the momenta, which are known as the Bethe Ansatz equations [14].
The Coordinate Bethe Ansatz method is a very useful analytic technique to diagonalize a transfer matrix [14]. For a general model, the transfer matrix is of orderkL, in whichL is the system size and k is some number which is related to
the local degrees of freedom. The large advantage of the Bethe Ansatz method is that it reduces the eigenvalue problem, for a possibly very large matrix, to a problem for solving a system of L coupled non-linear equations. A model is called solvable if such a reduction is possible [7].
Chapter 2
Construction of the
statistical lattice model
2.1
Square-triangle random tiling model
s1 s2 s3 t1 t2
Figure 2.1: The elementary tiles of the square-triangle random tiling model, as adopted from fig. 2.1 of [8].
In this thesis we will be examining the square-triangle random tiling model. The elementary tiles of our model are squares and equilateral triangles, which we take to have equal edge lengths. The two-dimensional plane is tiled is such a way that there are no gaps or overlaps between the tiles [8] and such that the edges adjoin [7]. The 3 elementary squares and 4 elementary triangles, which have the respective weights ofsi and tj, are shown in figure 2.1. Note that we
assign the same weight to a triangle of a given orientation and its π-rotated versions. Since, due to the imposed boundary conditions [14], their respective amounts are equal in number [8].
The square-triangle random tiling ensemble has three distinct symmetric phases: the quasicrystalline, the hexagonal and the square phase. Typical tilings of the plane for the three symmetric phases are shown in figure 2.2.
The quasicrystalline phase is characterized by a 12-fold rotational symmetry [8]. In this phase both the squares and the triangles, each cover half of the area of the plane [7]. Moreover, all the orientations of the squares (triangles) occur in the same number [8]. The other two phases are incommensurate.
The characteristic rotational symmetry for the hexagonal phase is 6-fold. In this phase all the orientations of the squares occur in the same number and one
Figure 2.2: Typical tiling configurations of the plane for the three symmetric phases of the square-triangle random tiling model: (A) hexagonal, (B) qua-sicrystalline and (C) square phase. Figure adopted from [1].
orientation of the triangles is dominant over the other. As can be seen from A in figure 2.2, the squares form domain walls between large hexagonal patches, consisting of the dominant triangle. The domain walls are connected through the triangles of the other orientation. Since either orientation of the triangles can dominate over the other, this phase is 2-fold degenerate.
The square phase is characterized by a 4-fold rotational symmetry. One orientation of the squares is dominant over the others and the orientations of the triangles occur in equal numbers. The domain walls between large rect-angular patches, consisting out of the dominant orientation of the squares, are constructed by the triangles. The other two orientations of the squares connect the domain walls. A typical configuration of the tiling for this phase is shown in C in figure 2.2. Since each orientation of the triangles can dominate over the other two, this phase is 3-fold degenerate.
However, in both incommensurate phases the large patches, of either the dominant triangle or square, do not have to be of the same size and shape [8].
The parameter space of this model is seven dimensional, since there are 7 different elementary tiles. Each of the seven parameters corresponds to a tile density, for which each (tile density) has a conjugate chemical potential. Fortunatelly the thermodynamic degrees of freedom can be reduced to three parameters. Two degrees of freedom can be removed, since for each orientation of the triangles theirπ-rotated versions occur in equal numbers. Another degree of freedom can be removed, since the total area of the tiling is constant. Fur-thermore the last degree of freedom can be removed by a geometrical constraint [8], which is given by
nt1nt2 = 4(ns1ns2+ns2ns3+ns1ns3) (2.1)
In which nx is the number of tiles of type x. This geometrical constraint is
2.2
Transformation of the square-triangle model
into a triangular lattice model
In this thesis we will examine the hexagonal phase of the square-triangle tiling. We will investigate this phase by applying the transfer matrix approach on a triangular lattice of unit edge length. In order to apply this transfer matrix approach, we will transform the tiles of the square-triangle model in such a way that they fit onto a triangular lattice. We set this transformation up, precisely, as follows: set the edges of the 7 elementary tiles to be of unit length and denote the weights of the tiles before the transformation with lower case letters and after the transformation with upper case letters. Rotate the edges of the tiles over an angle of π
12 to the left or to the right, in such a way that
the edges end up with an angle with the horizontal equal to the their closed multiple of π3. This procedure transforms the squares to rhombi, which occupy two adjacent triangles of the lattice, and rotates the triangles [14]. Since t1
and t2 are transformed to the same set of triangles, we decorate the triangles
of T1 or T2 with a Y-vertex, depending on whether the hexagonal patches are
constructed fromT2 or T1. We decorateS2 and S3 with a horizontal line and
decorateS1with a tilted line. The elementary tiles after the transformation are
depicted in figure 2.3. Furthermore, the decorated lines have to be connected to each other.
S1 S2 S3 T1 T1 T2 T2
Figure 2.3: Elementary tiles of the square-triangle random tiling model after the transformation to the lattice model.
Through the above described transformation for the square-triangle tiling, we have obtained a lattice model. In our triangular lattice model the decorated lines, which represent the domain walls, form a hexagonal net on top of a tri-angular lattice. We let the transfer matrix run along one of the directions of the decorated lines and define our vertical axis of the lattice, parallel to this direction. Then by interpreting the vertical axis of the lattice as a time axis, the transfer matrix can be interpreted as a time evolution operator and the decorated lines as world lines of particles. In our model we can identify left moving particles, right moving particles and bound states, which consist out of a left and a right moving particle. Bound states can split into a right and left moving particle and a right and a left moving particle can merge into a bound state [14].
Note that the partition function is unchanged, since each original tiling is mapped, in an unique way, into the deformed tiling [14]. The transformation to a triangular lattice model is possible, since the transformed tiles fit perfectly onto the triangular lattice. Moreover, the transformation is unique, since each elementary tile of the square-triangle random tiling model is transformed in an
unique way into the elementary tiles for out lattice model.
2.3
Construction of the transfer matrix for the
lattice model
Before constructing the transfer matrix, we will represent our lattice model in a graphically simpler way. Instead of drawing the transformed tiles on top of a triangular lattice, we will slightly transform the obtained decorated tiling and only draw the vertices of the new triangular lattice, occasionally the center of each horizontal edge (the horizontal distance between two vertices) and always the world lines of the particles. We denote each vertex as a dot, the center of each horizontal edge (mid edge) as a cross and the world lines link dots and/or mid edges together. An example of this new representation is given in figure 2.4. H H H H H H H H H H H H H H H H H H H HH H H H H H H H H H s s s s s s s s s s s s s s s s s s s s s s s s s s s s × s × s × s × s
Figure 2.4: Simplified representation of the transformed tiles on top of a trian-gular lattice.
The rows of the triangular lattice consist out of L edges. Hence our rows consist out of 2L points, of which there are L vertices and L mid edges. Since we impose periodic boundary conditions along the horizontal direction of the lattice, the first vertex is identified with the last vertex of each row. Therefore for each row starting from the bottom row, we label these points from 0 to 2L and identify point 2L with point 0.
We choose the flow of time to be in the upward direction and we divide the time axis in time steps. We define a time step as the distance in the vertical direction between a vertex and a mid edge. The bottom row is defined to be in the zeroth time step.
We assign an integer number from 0 to 5 to each point, this notation is used in the transfer matrix program given in appendix F. We will postpone the use of the transfer matrix program until chapter 4. Behind each integer number we denoted the labeling used in chapter 3 in between brackets, i.e. (j) for a just bound state. Then 0 denotes an empty point, 1 (o) denotes a bound state (at an odd number of time steps), 2 (r) denotes a point which is crossed by a right mover, 3 (l) denotes a point which is crossed by a left mover, 4 (j) denotes a bound states (which has just formed) and 5 (e) denotes a bound state (at an even number of time steps). Now each row of points can be denoted by a string of numbers, which we shall refer to as a configuration. Each configuration
consists out of 2L elements, of which each element can take the value 0 to 5. Now that we have a procedure for labeling the rows of the lattice, we can construct the transfer matrix. The construction goes as follows: we assign to each transition between two configurations a weightw for each bound state that continues one time step in the vertical direction, a weights for each bound state that splits into a right and a left moving particle and a weight of 1 for each of the following transitions: for each merger of a right and left moving particle into a just formed bound state, for each right (left) moving particle which propagates to the right (left) in one time step and for the vertical propagation of an empty point. For each transition between two configurations we multiply the weights, of the individual processes, that we have obtained through application of these rules. For example we assign the weights ∗ w for the case of 2 bound states, for which one state continues and one splits. Furthermore we impose an additional constraint that the bound state can fall apart only after an even number of time steps. This implies that a bound state, which is created at a vertex, can only split at a vertex. Since the weights w and s are positive and some elements of the transfer matrix will be zero, we obtain a nonnegative transfer matrix. According to the Perron-Frobenius theorem, the largest eigenvalue of this transfer matrix is real and nonnegative. The Perron-Frobenius theorem is summarized in appendix D.
We impose further constraints on the possible configurations. Each point in the lattice cannot be occupied by more than one right and one left moving par-ticle, or equivalently a bound state. Therefore configurations in which a point is occupied by bound state and a right and/or left moving particle and ones in which a left and right moving particle form a bound state which immediately disintegrates into a left and right moving particle are forbidden, or previous configurations which inevitably lead to these forbidden configurations. Further-more we forbid two bound states to be situated at an adjacent vertex and mid edge, since in the hexagonal phase the domain walls have a width of one square. We also restrict our starting configuration, which represents the row at the ze-roth time step, to consist out of just formed bound states, which are all situated at vertices and none at the mid edges (or visa versa). We make this choice in order to prevent the existence of forbidden states such as {2, 3} (the allowed state would be {2, 0, 3}, for example). These rules restrict the possible integer numbers which can be assigned to the vertices and the mid edges considerably. In light of these rules the integers 0 to 5, without 1, can be assigned to each vertex, however only the integers 0 and 1 can be assigned to the mid edges.
The dimension of our transfer matrix is 42L× 42L, since there are 2L points
in a configuration and each point can be empty, occupied by a bound state, occupied by a right moving particle or occupied by a left moving particle. How-ever, this dimension is much larger than necessary to include all the allowed configurations. This is the consequence of the imposed restrictions on all the possible configurations, given by simple combinatorics. In particular the mid edges can only take two values, as opposed to the vertices. Therefore large parts of the transfer matrix will be zero: large numbers of blocks will be completely empty and within nonzero blocks there will be empty rows and columns.
From section 2.1 we know that the parameter space of the square-triangle random tiling model is three dimensional. We choose the number of left moving particlesnl, the number of right moving particlesnr, the weightss and w to be
than necessary to our parameter space, since we do not know wheterw, s or a combination of the two weights is the most convenient third parameter. As we will conclude in chapter 5, the weight w will be sufficient and we do not need the weights.
Chapter 3
Coordinate Bethe Ansatz
3.1
Preliminaries to the Coordinate Bethe Ansatz
The number of right and the number of left moving particles are conserved under the action of the transfer matrix. Note that we defined a bound state to consist out of one right and one left moving particle. Therefore our transfer matrix adopts a block diagonal form. We label each block, which relates states with an equal number of right and an equal number of left moving particles, by its corresponding numbers of right and left moving particles. We will set up the eigenvalue equation for each block separately. And we adopt a similar scheme as in [14], to obtain the Bethe Ansatz equations.
Due to our division of each row in vertices and mid edges, we choose to measure our position variablesxiin units of one half edge (the distance between
a vertex and a mid edge), wherei runs from 1 to the number of given states the particles can be in. For example, we could denote the position of the only right moving particle in a configuration by x, or we could denote the right moving particle of a bound state at an even number of time steps byx.
We label the coefficients ψ(x) with a subscript, as follows: a left (right) moving particle is denoted by l [left] (r [right]), a bound state on a vertex by e [even], a bound state at a mid edge by o [odd] and a just formed bound state by j [join].
3.2
Derivation of the amplitude relations
The first block consist out of zero particles and has a dimension of 1 × 1. Since there are no particles, there is only one transition from the zeroth to the first configuration. The corresponding transfer matrix element is 1, therefore the eigenvalue for this sector isλ = 1.
For the block withnr= 1, we consider a right moving particle which
prop-agates one step to the right in one time step from position x − 1 to x. The corresponding eigenvalue equation reads
λψr(x) = ψr(x − 1) (3.1)
At this point in the calculation we apply the coordinate Bethe Ansatz. We ex-pect that this eigenvalue equation has a plane-wave solutionψr(x) = exp(ikx) =
ux, where k is the momentum of the right moving particle. The motivation for
considering a plane-wave solution is translational invariance of the system in the horizontal direction. Inserting the plane-wave solution into (3.1) results in the eigenvalue
λ = u−1 (3.2)
For the block withnl = 1, we consider a left moving particle which
prop-agates one step to the left in one time step from position x + 1 to x. The corresponding eigenvalue equation reads
λψl(x) = ψl(x + 1) (3.3)
Again we expect a plane-wave solutionψl(x) = exp(iqx) = vx, where q is the
momentum of the left moving particle, for this eigenvalue equation. Inserting the plane-wave solution into (3.3) results in the eigenvalue
λ = v (3.4)
For the block withnr=nl= 1, we have to consider two situations. Firstly
the situation in which the particles do not interact with each other, since they are at a large distance from each other. Secondly the situation in which the particles do interact with each other, resulting in a just formed bound state.
In the first situation the particles behave as if they were isolated from each other. Therefore the resulting eigenvalue equations are
λψlr(y, x) = ψlr(y + 1, x − 1) with y < x − 2 (3.5)
λψrl(x, y) = ψrl(x − 1, y + 1) with x < y (3.6)
Since the particles behave as if they were isolated from each other, we consider the piecewise plane-wave solution
ψlr(y, x) = Alruxvy fory < x − 2 (3.7)
ψrl(x, y) = Arluxvy forx < y (3.8)
Inserting the piecewise plane-wave solution into the eigenvalue equations (3.5) and (3.6) results, for each of the eigenvalue equations, in the eigenvalue
λ = v
u (3.9)
For the second situation the eigenvalue equation reads
λψj(x) = ψrl(x − 1, x + 1) (3.10)
We consider the plane-wave solutionψrl(x − 1, x + 1) = Arlux−1vx+1 and insert
the eigenvalue (3.9) for this block into this eigenvalue equation, resulting in the relation between a just formed bound state and a right and left moving particle ψj(x) = ψrl(x, x) (3.11)
To complete the investigation of this block we have to consider the eigenvalue equations for the following transitions: the propagation of a bound state at a mid edge to a vertex, the separation of a bound state into a right and left moving
particle and the propagation of a bound state at a vertex and the propagation of a just formed bound state to a bound state at a mid edge. The eigenvalue equations are
λψe(x) = wψo(x) (3.12)
λψlr(x − 1, x + 1) = sψe(x) (3.13)
λψo(x) = wψe(x) + wψrl(x, x) (3.14)
We consider the plane-wave solutionψsub(x, y) = Asubuxvy, in whichsub ∈
{e, o, lr, rl}. Inserting these solutions into equations (3.12),(3.13) and (3.14) results in the following relations between the amplitudesAsubof the plane-wave
solutions v uAe=wAo (3.15) Alr =sAe (3.16) v uAo=wAe+wArl (3.17) As will be explained in section 3.3 we need the amplitude relations between the interchange of a right and a left moving particle and the interchange of two like particles, in order to derive the Bethe Ansatz equations. We can derive the amplitude relation for the interchange of a right and a left moving particle by exploiting equations (3.15), (3.16) and (3.17). Inserting (3.15) and (3.16) into (3.17) and eliminatingAeandAo, results in the desired amplitude relation
Arl= v2 u2w2s− 1 s Alr (3.18)
In order to deduce the amplitude relations for the exchange of two like parti-cles, we need to consider the block with two like particles and one unlike particle. We are obliged to investigate this sector since two like particles cannot interact with each other without one unlike particle assisting in the process.
We start with the calculation of the eigenvalues for the sector with 3 particles, of which two are like particles. To calculate the eigenvalues, we consider the case in which all particles are spatially separated and hence do not interact with one another.
For the block withnr= 1 and nl= 2 we have the eigenvalue equation
λψ12r(x, x + 2, x + 6) = ψ12r(x + 1, x + 3, x + 5) (3.19)
Where 1 and 2 label the left moving particles, each with a different momentum. We consider the the plane-wave solution ψ12r(y1, y2, x) = A12rvy11v
y2
2 ux, which
inserted into the above eigenvalue equation results in the eigenvalue λ = v1v2
u (3.20)
Note that we obtain exactly the same eigenvalue (3.20) if we had interchanged labels 1 and 2 in (3.19). For the block with nr = 2 and nl = 1 we have the
eigenvalue equation
Where, this time, 1 and 2 label the right moving particles, each with a dif-ferent momentum. We consider the the plane-wave solution ψl12(y, x1, x2) =
Al12vyux11u x2
2 , which inserted into the above eigenvalue equation results in the
eigenvalue
λ = v u1u2
(3.22) Here we would also obtain exactly the same eigenvalue (3.22) if we had inter-changed labels 1 and 2 in (3.21).
We are now able to investigate the sector in which the 3 particles interact with each other. Once more we start with examining the sector nr = 1 and
nl = 2. For this sector we expect a wave function of the piecewise plane-wave
type, wherey1,y2, x are in increasing order, which are given by
ψllr(y1, y2, x) = A12ullr v y1 1 v y2 2 u x+A21u llr v y1 2 v y2 1 u x (3.23) ψlrl(y1, x, y2) =A1u2lrl v y1 1 uxv y2 2 +A2u1lrl v y1 2 uxv y2 1 (3.24) ψrll(x, y1, y2) =Au12rll v y1 1 v y2 2 u x+Au21 rll v y1 2 v y2 1 u x (3.25)
For the situation in which the 3 particles interact with each other we need to use these 3 structures for the wave function. However there are two situations in which we can treat the two terms in the above wave functions separately. The first situation is for the calculation of the eigenvalues, since all the particles are spatially separated. The second situation is for the generalization of (3.18), since two unlike particles interact while the other particle (the spectator) is spatially separated from the interaction.
First, consider the calculation of the eigenvalues. By examining the struc-ture of the eigenvalues, for the case we have considered, we observe that each eigenvalue factorizes in factors for each momentum. Therefore we expect that we can simply generalize the eigenvalue for an arbitrary number of left and number of right movers to obtain
λ = nl Y j=1 vj nr Y k=1 u−1k ! (3.26) Secondly, consider the case of two unlike particles interacting with each other and a spatially separated spectator particle. This spectator particle does not influence the interaction process and can therefore be treated independently. This is reflected by the observation that the contribution of the interacting particles and that of the spectator particle factorize. Therefore we can gen-eralize equations (3.15)-(3.18) by considering the situation of two interacting unlike particles and a spectator particle. This result in the generalization of the amplitude relation (3.18) to A...p1p2... ...rl... = (v p2) 2 (up1)2w2s −1 s A...p2p1... ...lr... (3.27)
Wherep1andp2denote the momenta of the particles.
We are now able to deduce the amplitude relations for the exchange of two like particles. We examine the situation in which all three particles interact with each other and we start by examining the block withnr = 1 andnl = 2.
H H H H H H H H H H H H H H H H H H H HH H H H s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s s
Figure 3.1: Configuration for the derivation of the amplitude relation for the exchange of two like particles, modified picture of B. Nienhuis.
The final configuration consist out of a left moving particle at positionx − 4, a left moving particle at positionx and a right moving particle at x + 2. The previous history of this configurations is unique for up to six time steps. As our starting configuration we choose the configuration four time before the final configuration. In the starting configuration there is a bound state at x and a left mover atx + 2. The eigenvalue for this transition is given by
λ4ψllr(x − 4, x, x + 2) = s2w2ψel(x, x + 2) (3.28)
Now by inserting, (3.20), (3.23) and a similar equation forψel into (3.28) we
obtain v1v2 u 4 A1,2,rvx−41 v2xux+2+A2,1,rv2x−4v1xux+2 = s2w2 Ae1,2u xvx 1v2x+2+Ae2,1u xvx+2 1 v2x Simplification of this equation and inserting the generalization of (3.16) in
the form of Ae1,2 = s
−1A
1r2 and Ae2,1 = s
−1A
2r1 (in order to eliminate the
label e), results in v1v2 u 4 A1,2,rv−41 u 2+A 2,1,rv2−4u 2 = sw2 A 1,r,2v22+A2,r,1v12 (3.29) Now we make use of (3.27) in the formA1,r,2=
v2 2
w2u2s−1s
A1,2,r, resulting
in the amplitude relation for the interchange of two left moving particles A21r = −A12r
v2
v1
2
(3.30) Note that, although the right moving particle is essential in order for the process of interchanging two left moving particles to take place, it does not enter this amplitude relation.
We can do a similar computation for the block with nr = 2 and nl = 1,
in order to find the amplitude relation for the interchange of two right moving particles. Starting from the eigenvalue equation
λ4ψ
We obtain, after similar computations, the amplitude relation for the inter-change of two right moving particles
Al21=Al12
u2
u1
2
(3.32) From which we note that the left moving particle does not enter this amplitude relation, although it is necessary for the process of interchanging two right moving particles.
Since we can consider the situation in which three particles interact, spatially isolated from other particles, we can generalize the amplitude relations (3.30) and (3.32) for the interchange of two like particles to
A...p2p1... ...ll... = −A ...p1p2... ...ll... vp2 vp1 2 (3.33) A...p2p1... ...rr... = −A ...p1p2... ...rr... up2 up1 2 (3.34) The structure of the considered wave functions, inspires us to construct an Ansatz for the sector with arbitrary particle numbers. This Ansatz is a linear combination of products of plane waves and is given by
ψc(x)= X p Ac p Y j (wcj pj) xj (3.35)
In which x = (x1, x2, ...) represents the sequence of positions in increasing
order (xj< xj+1). The sequence of particle types are encoded in c = (c1, c2, ...),
in which nr elements are right movers and nl elements are left movers. The
summation is over the permutations p = (p1, p2, ...) of the indexes, only mixing
particles of the same type. pi denotes the momentum of a particle. And the
exponential wave numbers are denoted as wr
j ≡ uj and wjl ≡ vj. Note that we
interchanged the momentum and particle types in the notation.
We have to remark that the momenta of particles of the same type must all be different in order to have a non vanishing wave function. This can be derived from the generalized amplitude relations. For example one can observe from (3.33) and (3.34) that when the momenta of two like particles are the same, the amplitudes and hence the wave function vanishes.
However, as a final note, we have not proven that the the eigenvalue equa-tions, the Ansatz and the amplitude relations are valid for arbitrary particle numbers and arbitrary configurations. In order to proof that the Ansatz is valid in general, one has to show that each group of interacting particles can be treated independently. And one has to consider arbitrary (horizontal) lengths of such groups of interacting particles, i.e. consisting out of an arbitrary number of particles. Also one has to consider a arbitrary number of these groups of interacting particles.
3.3
Derviation of the Coordinate Bethe Ansatz
equations
Since we have imposed periodic boundary conditions on our system, we want the Ansatz to be consistent with these periodic boundary conditions. Therefore
the Ansatz has to satisfy
ψ...,c(..., 2L) = ψc,...(0, ...) (3.36)
Which requires that the amplitudes have to satisfy the relation A...c ...p(w c p) 2L=Ac... p... (3.37)
We can use equations (3.27), (3.33) and (3.34) to relate the right and left hand side of this equation, by successive interchanges of the momenta and parti-cle types. This will enable us to eliminate the amplitudes. We will then obtain consistency equations on the momenta, which are known as the Coordinate Bethe Ansatz equations (BAE). Since we have two types of particles, we will obtain two sets of BAE, one set for each type of particle.
To derive the first BAE we examine (3.37), for the right moving particles Alllrrr
123123(u3)2L=Arlllrr312312 (3.38)
By applying (3.27) three times to this equation, we obtain Alllrrr 123123(u3)2L= v2 1 u2 3w2s −1 s v2 2 u2 3w2s −1 s v2 3 u2 3w2s −1 s Alllrrr 123312 (3.39)
By applying (3.34) three times to this equation, we obtain Alllrrr123123(u3)2L= v2 1 u2 3w2s −1 s v2 2 u2 3w2s −1 s v2 3 u2 3w2s −1 s (−1) u3 u1 2 (−1) u3 u2 2 Alllrrr123123 (3.40) Which, after canceling the amplitudes, reduces to
(u3)2L= (−1)2 3 Y k=1 ( vk)2 (u3)2w2s −1 s 3 Y i=1 u3 ui 2 (3.41) Which can be generalized to the first BAE given by
(uj)2L= (−1)nr+1 nl Y k=1 (v k)2 (uj)2w2s −1 s nr Y i=1 uj ui 2 (3.42) To derive the second BAE we examine (3.37), for the left moving particles
Arrrlll
123123(v3)2L=Alrrrll312312 (3.43)
By applying (3.27) three times to this equation, we obtain Arrrlll123123(v3)2L =Arrrlll123312 (v 3)2 (u1)2w2s −1 s −1 (v 3)2 (u2)2w2s −1 s −1 (v 3)2 (u3)2w2s −1 s −1 (3.44) By applying (3.33) three times to this equation, we obtain
Arrrlll 123123(v3)2L=Arrrlll123123 3 Y i=1 (v 3)2 (ui)2w2s −1 s −1 (−1)2 v3 v1 2 v 3 v2 2 (3.45)
Which, after canceling the amplitudes, reduces to (v3)−2L = 3 Y i=1 (v 3)2 (ui)2w2s −1 s (−1)2 3 Y k=1 vk v3 2 (3.46) Which can be generalized to the second BAE given by
(vk)−2L= (−1)nl+1 nr Y j=1 ( vk)2 (uj)2w2s −1 s nl Y i=1 vi vk 2 (3.47) Now that we have obtained the Coordinate Bethe Ansatz equations, we can attempt to obtain the solutions to these equations by means of numerical analyses. We will do this in chapter 4. Each solution to the BAE provides an eigenvalue which is obtained by inserting the solutions into the equation for the eigenvalue (3.26).
3.4
Logarithmic Coordinate Bethe Ansatz
equa-tions
It is convenient to make the change of variablesu2=τ and v2= in equations
(3.42) and (3.47), this results in the following Bethe Ansatz equations: τL j = (−1) nr+1 nl Y k=1 k τjw2s −1 s nr Y i=1 τj τi (3.48) −Lk = (−1) nl+1 nr Y j=1 k τjw2s −1 s nl Y i=1 i k (3.49) And the eigenvalue equation (3.26) transforms to
λ = nr Y j=1 τj−1/2 nl Y k=1 1/2k ! (3.50) By applying the natural logarithm to both sides of the BAE (3.48) and (3.49), we obtain the logarithmic BAE, which are given by
LF+(τj) = (nr+ 1)πi mod 2πi (3.51)
LF−(k) = (nl+ 1)πi mod 2πi (3.52)
Where the F-functions are given by: F+(z) = (1 + nl− nr L ) log(z) − 1 L nl X k=1 log(k− zw2) + nl L log(w 2s) +1 L nr X i=1 log(τi) (3.53) F−(z) = (−1 + nl L) log(z) − 1 L nr X j=1 log(z − τjw2) − log(τjw2s) − 1 L nl X i=1 log(i) (3.54)
We can conclude that, from the form of the logarithmic BAE, it is obvious that the roots of the BAE are situated on curves in the complex plane, for which the real part of the F-functions F+ and F−, respectively, vanish: i.e.
Re[F+(τ )] = Re[F−()] = 0.
In order to obtain the largest eigenvalue, we expect the solutions to be as close to each other on the solution curves as possible. Therefore we expect that the imaginary part of the F-functions, between two successive roots, changes by
2π
L and not by a multiple amount, hence
F+(τj+1) −F+(τj) = 2πi L (3.55) F−(k+1) −F−(k) = 2πi L (3.56)
Chapter 4
Numerical analyses
4.1
Application of the Newton-Raphson method
to the Bethe Ansatz equations and
numeri-cal analyses of the transfer matrix
In this chapter we apply numerical analyses to the BAE (3.48) and (3.49) in order to find the numerical solutions to these equations. We employ the Newton-Raphson method, which is a root finding algorithm to find the approximate roots of a set of functions. A thorough explanation of this method is presented in appendix E.
In order to apply the Newton-Raphson method to the Bethe Ansatz equa-tions, we first rewrite these equations to
τL j − (−1) nr+1 nl Y k=1 k τjw2s −1 s nr Y i=1 τj τi = 0 (4.1) −Lk − (−1)nl+1 nr Y j=1 k τjw2s −1 s nl Y i=1 i k = 0 (4.2) The first BAE gives rise tonr equations, the second BAE gives rise to nl
equations. Both equations depend on the parameters L, w, s,nr,nl,nrdifferent
τ -variables and nl different-variables. In order to find the roots, τj andk, to
the BAE we fix the parameters L, w, s,nrandnl. In our numerical analyses we
set the number of right movers equal to the number of left movers, i.e. nr=nl.
This choice will make it easier to find the roots.
To perform the numerical analyses, we have written a program in Math-ematica 8, which is given in appendix G. This program applies the Newton-Raphson method to the Bethe Ansatz equations (4.1) and (4.2) for the given set of fixed parameters. Moreover this program calculates the (largest) eigenvalue of the given solution set of the roots and it calculates the difference between the change in the imaginary part of the F-functions for two successive roots. How-ever, we verified in another Mathematica 8 file that the roots (corresponding to the largest eigenvalue) satisfy the BAE.
Furthermore, we have written programs for calculating the transfer matrix and the corresponding largest eigenvalue, for a given starting configuration and
for the given set of fixed parameters. These programs are given in appendix F. We can calculate the transfer matrix and corresponding largest eigenvalue for starting configurations consisting out of 1 or 2 bound states. However we did not write programs for starting configurations consisting out of a higher number of bound states. Since we can calculate the largest eigenvalue of a system by solving the BAE or by using the transfer matrix, the programs for the transfer matrix give us an independent way to calculate the largest eigenvalue of the system. Therefore, by using the transfer matrix programs, we can verify that we have obtained the correct largest eigenvalue of the BAE (for a given choice of fixed parameters) and we can verify that our obtained BAE are correct. Moreover, if both methods give an identical largest eigenvalue for a particular system, we can conclude with large certainty that the BAE are correct and that we have indeed obtained the largest eigenvalue.
4.2
Parameter space of the system
Unfortunately, we expect that our Bethe Ansatz equations (3.48) and (3.49) are to complicated to solve for general values of the parameters. Therefore we are going to specify in which part of the parameter space we expect the BAE to be solvable. This restriction of the parameter space will simplify the BAE, in such a way that we expect to be able to solve these BAE in the thermodynamic limit.
First we transform our discrete lattice model to a continuous lattice model, by taking the system size L to infinity. Remember that we had set our lattice spacing to be of unit length. As a consequence we get rid of some of the microscopic physics and obtain a simpler model. However, we will not start making use of the continuum limit until section 5. In this chapter we will apply the numerical methods on the discrete lattice model.
Secondly, we put some restrictions onto the average length of the domain walls. We have seen in section 2.3 that the domain walls occur in three different directions. We define an average length of the domain walls, along each of the three different directions. If the hexagonal net is isotropic, the average lengths of the domain walls are all equal to each other. However, if the hexagonal net is not isotropic, the average lengths of the domain walls are not all equal to each other. We represent the average lengths of the domain walls graphically in a phase diagram, see figure 4.1. This phase diagram is used as follows: draw a point inside the triangular phase diagram, then draw a line from the point to each of the edges of the triangle. Each line should be perpendicular to an edge. Then the length of each line represents the average length of the domain walls along the direction of the line. Now we choose to examine the partition function along the bottom edge of the phase diagram. Along this edge, the average length of the domain wall in the vertical direction approaches zero. In order to achieve this situation in the lattice model, we will need to stimulate the disintegration of the bound states and oppress the propagation of the bound states. We achieve this by taking the limit of the weight s to infinity and of the weight w to zero, respectively. Furthermore one can move along the bottom edge of the phase diagram by adjusting the number of right and the number of left movers.
Figure 4.1: Phase diagram for the average size of the domain walls. The blue dashed lines represent the average length of the domain walls in the correspond-ing direction.
4.3
Results of the numerical analyses on the Bethe
Ansatz equations and the transfer matrix
Since we are examining the hexagonal phase of the square-triangle tiling, we choose the system size in our numerical calculations to be significantly larger than the number of particles. If we do not make this choice and let the system size be of similar magnitude to the number of particles, the system would be be very close to the quasicrystalline phase. Examining this section of the phase di-agram makes it very hard to obtain the solutions to the BAE, which correspond to the largest eigenvalue of the system. By choosing the number of particles to be significantly smaller than the system size, the hexagonal patches will be relatively large. And as a consequence the roots of the BAE can be found more easily.
However, in stead of taking our system size to be significantly larger than the number of bound states, we start by setting the system length L equal to 6. We make this choice for the system size in order to keep the computation time for the transfer matrix programs to a minimum. The computation time, as well as the size of the transfer matrix itself, increases substantially by increasing the system size. By keeping the computation time to a minimum, we can quickly verify that the solutions to the BAE (obtained with the Newton-Raphson method), corresponding to the largest eigenvalue, are correct. Fortunately, we can always acquire the solutions to the BAE, for a larger system size, from the solutions to the BAE for the smaller system size. These new solutions can be obtained by gradually increasing L, using the old solutions as the initial guess and then applying the Newton-Raphson method on the BAE, obtaining the new solutions. This process can be continued until the solutions for the desired system size are obtained. This procedure is valid since, although only the solutions with integer L correspond to physical systems, the BAE are valid mathematical equations for any L, even for non integer values. Note that we can also apply this scheme to the statistical weights of the system.
We will refer to the solutions to the BAE corresponding to the largest eigen-value, simply, as the solutions to the BAE. We start our numerical analyses of the BAE with the following choice for the parameters: L = 6, w = 2.5 and
