9512.net

# Numerical study of the 6-vertex model with domain wall boundary conditions_图文

arXiv:cond-mat/0502314v1 [cond-mat.stat-mech] 13 Feb 2005

Numerical study of the 6-vertex model with domain wall boundary conditions
David Allison and Nicolai Reshetikhin February 2, 2008
Abstract A Markov process is constructed to numerically study the phase separation in the 6-vertex model with domain wall boundary conditions. It is a random walk on the graph where vertices are states and edges are elementary moves. It converges to the Gibbs measure of the 6-vertex model. Our results show clearly that a droplet of c vertices is created when Boltzamnn weights are in the antisegnetoelectric region. The droplet is a diamondlike shaped curve with four cusps.

Introduction
It is well known that the 6-vertex model is exactly solvable and has phase transitions. The history and the classi?cation of phases in the 6-vertex model as well as many interesting facts about the structure of the partition function of the 6-vertex model with periodic boundary conditions can be found in [10], [3]. There is an important function of Boltzmann weights of the model which is usually denoted by ? (see [3],[10] and section 3). The 6-vertex model with periodic boundary conditions has 3 phases in the thermodynamical limit, depending on the value of ?. One is the totally ordered (frozen) phase, with ? > 1, the second is the disordered (critical) phase, with |?| < 1, and the third is the partially ordered (antisegnetoelectic) phase with ? < ?1. The 6-vertex model with domain wall boundary conditions on a square N ×N grid perhaps was ?rst considered in [?] in the process of computation of norms of Bethe vectors. The partiction function of this system can be written as the determinant of a certain N × N matrix [5]. Its asymptotics in the thermodynamic limit N → ∞ were analyzed in [7]. It is related to matrix models, which was pointed out and exploited in [17]. The 6-vertex model with ? = 1/2 is also known as the ice-model. This model with domain wall (DW) boundary condition is closely related to the enumeration of alternate sign matrices [9]. It also has other interesting combinatorial features (see for example [18], [15]). When ? = 0 the 6-vertex model is equivalent to the

1

problem of counting of weighted tilings of the Aztec diamond (see for example [7], [6] and references therein). The spatial coexistence theory of di?erent phases and the interfaces separating phases is an important part of statistical mechanics. Growth of crystals is one of the well known phenomena of this type. This is also closely related to the limit shape e?ect in statistics of Young diagrams [16] and plane partitions. In dimer models related to enumeration of plane partitions and domino tilings, the interface between the disordered and totally ordered phases is also known as an arctic circle phenomenon [1]. Dimer models on bipartite planar graphs with periodic weights are exactly solvable models where this phenomenon has been studied in [11]. In dimer models the limit shapes or interfaces (curves separating phases), under broad conditions, are real algebraic curves [11]. Since at ? = 0 the 6-vertex model is equivalent to a dimer model these results imply that such a phenomenon exists in the 6-vertex model for ? = 0. The natural question is whether the spatial coexistence of phases happens only at the free fermionic point or if it occurs for all values of ?. Numerical evidence suggesting the existence of a limit shape in the 6-vertex model with domain wall boundary conditions for all weights was obtained in [14]. Here we report results of numerical study of the 6-vertex model with DW boundary conditions in all phases of the model. Our method is di?erent from [14]. To generate a random con?guration in the 6-vertex model we construct the Markov process which is equivalent to a random weighted walk on the graph where the vertices are states of the model and edges are local moves which transform states into other states. This process satis?es the detailed balance condition and therefore converges to the Gibbs state of the 6-vertex model. It is also known as Monte-Carlo with local update and as a heat-bath algorithm. In statistical mechanics such processes are known as Kawasaki, or Glauber dynamics. For a more e?ective version of this algorithm known as the “coupling from the past” algorithm, see [12]. For periodical boundary conditions the system can be either in the ordered (segnetoelectric) phase, disordered (critical) phase, or antisegnetoelectric (noncritical) phase, depending on the values of Boltzmann weights. Our results con?rm the conclusion from [4] and [14] that there is a coexistence of ordered and disordered phases in the 6-vertex model. They also clearly indicate that for ? < ?1 there is a coexistence of all three phases. The outer layer is an ordered phase. It follows by the ring of disordered phase. Finally, there is an inner droplet of the antisegnetoectric phase. This phenomenon was ?rst conjectured in [14] using a di?erent numerical method. The shape of the inner droplet has four cusps and is reminiscent of one of the limit shapes for dimers on a square-octagon grid [11] equivalent to the diablo tiling . Acknowledgments. We are grateful to R. Kenyon and A. Okounkov for many illuminating discussions, to K. Palamarchuk for valuable comments, and to T. Yates for help with the implementation of the algorithm in the C language.

2

1
1.1

Weights and local moves
States

States of the 6-vertex model on a square lattice are con?gurations of arrows assigned to each edge. The 6-vertex rule is that the total number of arrows coming into any vertex should be equal to the total number of arrows going out of this vertex. Each con?guration of arrows can be equivalently regarded as a con?guration of empty edges (arrows oriented South-North and East-West) and occupied edges, or thick edges (arrows pointing in the opposite directions). It is clear that thick edges will form paths. Possible con?gurations of paths around a vertex are shown on ?g. 1.1. a1 b1 c1

a2

b2

c2

Figure 1: The 6 vertices and their weights We will use a1 , a2 , b1 , b2 , c2 , and c2 as names of the vertices. We denote by the same letters Boltzmann weights assign to these vertices. Domain wall boundary conditions are shown on ?g. 2 and 3. For domain wall boundary conditions every path in a 6-vertex con?guration will have one end at the North boundary of the square and the other end at the East boundary of the square. These paths can be regarded as level curves of a height function. The lowest height function is shown on ?g. 3 and the highest height function is shown on ?g. 2.

1.2

Weights

The weight of a state is the product of weights of vertices and the weight of a vertex is determined by rules from ?g. 1.1. The partition function is sum of

3

DW BChigh

. . ?. ?b1 ? ?b1 ? = ?b1 ? ?b1 ? ?b1 c2 ?

b1 b1 b1 b1 c2 b2

b1 b1 b1 c2 b2 b2

b1 b1 c2 b2 b2 b2

b1 c2 b2 b2 b2 b2

c2 b2 b2 b2 b2 b2

?

Figure 2: Domain Wall Boundary Condition High

.. . ...

? ? ? ? ? ?= ? ? ? ?

DW BClow

? ? ? ? ? ? =? ? ? ? ? ? ?

?

..

. c2 a2 a2 a2 a2 a2 a1 c2 a2 a2 a2 a2 a1 a1 c2 a2 a2 a2 a1 a1 a1 c2 a2 a2 a1 a1 a1 a1 c2 a2 a1 a1 a1 a1 a1 c2

?

.. . ..

.

? ? ? ? ? ? ?= ? ? ? ? ? ?

Figure 3: Domain Wall Boundary Condition Low

..

.

4

weights of all con?gurations: Z6v =
states vertices

w(vertex, state)

where w(vertex, state) is the weight of the vertex (see ?g.1.1). The ratio vertices w(vertex, state) Z6v is the probability of the state. This is the Gibbs measure of the 6-vertex model.

1.3

Local moves and the graph of states

Now let us describe local moves in the space of states. Such a move changes the con?guration of arrows at the minimal number of edges near a given vertex and it acts transitively, i.e. any given state of the model can be transformed to any other given state of the model by a sequence of such moves. Such moves are most transparent in terms of height functions. There are two types of local moves: ? The path from ?g. 4 we can move up, i.e. to the path from ?g. 5. We will call this move ?ip up. ? The path from ?g. 5 can be moved down, i.e. to the path from ?g. 4. We will say that this is the ?ip down.

Sa =

b1 c2

a2 b2

=

Figure 4: A local con?guration that may be ?ipped up

Sb =

c2 a2

c1 c2

=

Figure 5: The local con?guration after an up ?ip has occurred Such moves with all possible surrounding con?gurations we will call ?ips up and down. For each ?ippable vertex we introduce e?ective weight as follows: 5

? For a vertex ?ippable up the e?ective weight is the product of weights of all vertices that can be a?ected by the ?ip, i. e. the vertex itself, and the neighboring vertices to the North , the East, and the North-East of it. ? Similarly for a vertex ?ippable down the e?ective weight of it is the product of weights of the vertex itself, and of next neighboring vertices to the South, West, and South-West of it. The e?ective weight is always the product of four factors. The e?ective weight of vertex v in the con?guration S we denote by Wv (S ).

2
2.1

The Markov process
General strategy

Consider the abstract graph with vertices being states of the model and with edges being local moves. This graph is clearly connected. Our goal is to construct a random walk on this graph converging to the probabilistic measure vertices of this graph which is the Gibbs measure of the 6-vertex model with DW boundary conditions. Let us recall some basic facts. Let Γ be a ?nite connected graph and q : V (Γ) → R+ be a probabilistic measure on the set of vertices of Γ. Let M = {p(a → b)}a,b∈V (Γ) be the matrix of the Markov process describing a random walk on Γ. A traveller moves from a to b with the probability p(a → b). The matrix M must satisfy the total probability condition: p(a → b) = 1

b

If in addition it satis?es the detailed balance condition q (a)p(a → b) = q (b)p(b → a) then it is known that the Markov process converges to q . For details about Markov sampling and estimating convergence times, see [13]. Now our goal is to construct such random walk converging to the Gibbs state of the 6-vertex model. At some point the rate of convergence of this Markov process becomes an important issue. To avoid the complicated analysis of estimating mixing times we will modify the algorithm and will use the “coupling from the past” version. This will be explained later.

2.2

The Markov process for the 6-vertex model

We want to construct Markov process which chooses a vertex at random, then with the probability which we will describe below it will either ?ip the con?guration up at this vertex, or will ?ip it down, or will do nothing. The

6

probability of passing from the state Sa to the state Sb in this process can be w written as follows: P (Sa =? Sb ) = 1 # vertices 1 # vertices Pv (Sa =? Sb ) (1)

v

=

δ (Sa , Sb ) +
(v )non?f lip (v )f lip?up only

Pv (Sa =? Sb )+ (2) Pv (Sa =? Sb )

(v )f lip?down only

Pv (Sa =? Sb ) +

(v )bi?f lip

=

1 (# non-?ippable inSa )δSa ,Sb + # vertices

(v )f lip?up only

Pv (Sa =? Sb )+ (3)

(v )f lip?down only

Pv (Sa =? Sb ) +

(v )bi?f lip

Pv (Sa =? Sb )

=

#non-?ippable δSa ,Sb + # vertices #?ippable #vertices 1 #?ippable 1 #?ippable Pv (Sa =? Sb )+ 1 #?ippable

(v )f lip?up only

(v )f lip?down only

Pv (Sa =? Sb ) +

(v )bi?f lip

Pv (Sa =? Sb )

(4) Here #?ippable is the number of ?ippable vertices and #verices is the total number of vertices. Algorithmically, this means that we do the following: 1. With probability P =
#non-?ippable #vertices ,

do nothing (that is, restart the loop.)
#?ippable #vertices ,

+ #bi-?ip 2. With probability P = #?ip-up-only + #?ip-down-only = #vertices continue to the next part.

If the algorithm continues, select a ?ippable vertex with the probability: P (selection) = 7 1 #f lippables (5)

At this selected vertex the con?guration can be either ?ippable only up, or only down, or in both directions. Depending on this proceed according to the following rules: Three possible conditions now exist: 1. The vertex is ?ippable down only. Two options: ? Flip vertex down with probability Pv (Sa =? Sb ) = ρWv (Sb ) ? Stay with probability Pv (stay ) = 1 ? ρWv (Sb ) 2. The vertex is ?ippable up only. Two options: ? Flip vertex up with probability Pv (Sa =? Sb ) = ρWv (Sb ) ? Stay with probability Pv (stay ) = 1 ? ρWv (Sb ) 3. The vertex is ?ippable up and down. Three options: ? Flip vertex down with probability Pv (Sa =? Sb ) = ρWv (Sb ) ? Flip vertex up with probability Pv (Sa =? Sb′ ) = ρWv (Sb′ ) ? Stay with probability Pv (stay ) = 1 ? ρWv (Sb ) ? ρWv (Sb′ )

Here Wv (Sb′ ) Wv (Sb ) are the e?ective weights of the vertex v in the states obtained by ?ipping up or down at this vertex from the state Sa . E?ective weights were described in section 1.3. The parameter ρ is chosen such that all probabilities of transitions should be positive. In other words it should satisfy all conditions ρ < Wv1 (S ′ ) where v is a vertex ?ippable in the state S either up or down, but not bi?ippable and S ′ is the con?guration after the ?ip. At every bi?ippable vertex in the state S we 1 ′ should have ρ < Wv (S ′ )+ W (S ′′ ) where S is the result of the ?ipping S up at v and S ′′ is the result of the ?ip down. This process satis?es the detailed balance condition, and the total probability condition. Since the graph of states with edges being local moves is connected, this process converges to the Gibbs state of the 6-vertex model. The process also depends on the choice of ρ. It slows down when ρ is small.

8

3
3.1
3.1.1

Random states in the 6-vertex model with DW boundary conditions
Phases in the 6-vertex model

One can write weights of the 6-vertex model as a1 b1 c1 = = = E1 + Hx + Hy E1 ? Hx ? Hy ), a2 = exp(? ) T T E2 ? Hx + Hy E2 + Hx ? Hy ), b2 = exp(? ) exp(? T T E3 c2 = exp(? ) T exp(? (6) (7) (8)

Here E1 , E2 , and E3 are energies of the interaction of arrows (or energies associated with the local shape of level curves of the height function) and Hx and Hy are magnetic ?elds. In this interpretation arrows can be regarded as spins interacting with the magnetic ?led such that the energy of a vertical arrow is ±Hx depending on whether the arrow is heading up or down. The energy of a horizontal arrow is ±Hy depending on whether it is oriented left or right. We assigned the energy of an arrow to the energies of adjacent vertices. Notice that since the total number of c1 - and c2 -vertices satisfy the relation n(c1 ) ? n(c2 ) = N the partition function changes by an overall factor only when we change c1 /c2 . The total numbers of a and b vertices satisfy similar relations: n(a1 ) = n(a2 ) and n(b1 ) = n(b2 ). Because of this for the square lattice with DW boundary conditions we can set a1 = a2 = a, b1 = b2 = b, and c1 = c2 = c without loosing generality. 3.1.2 Let us recall the phase diagram of the 6-vertex model [10], [3] with periodic boundary conditions in the absence of magnetic ?elds. The important characteristic of the model is the parameter ?= a 2 + b 2 ? c2 2ab

The phase diagram for the 6-vertex model with periodic boundary conditions in the absence of magnetic ?elds is shown on ?g. 6. There are four phases: 1. Phase I: a > b + c(? > 1). This an ordered phase where there are two possibilities for the ground state. It either consists of a1 -vertices or of a2 vertices. In either case any change in the ground state gives the state with the total number of b and c vertices comparable with the linear size N of 9

Figure 6: The phase diagram for the six-vertex model in terms of the weights of a and b, assuming c = 1.

the system. Thus, as N → ∞ the energy of these two ground states is macroscopically separated from the energy of other states. In other words these are two frozen ground states. 2. Phase II: b > a + c (? > 1). This is an ordered phase with double degeneracy of the ground state. The ?rst possibility is when all vertices are b1 vertices, the second possibility is when all vertices are b2 -vertices. As in case of phase I, this is a frozen phase.
b+c (|?| < 1). This is a disordered phase. Local 3. Phase III: a, b, c < a+2 correlation functions decay as a power of the distance in this phase. These are the values of a, b, c when |?| < 1. In particular, the free fermionic curve ? = 0 lies entirely in this phase. It is shown by the dotted segment of the circle on ?g. 6.

4. Phase IV: c > a + b (? < ?1). This is an ordered phase with so-called antisegnetoellectic ordering (see ?g. 7). The ground state in this case consists of alternating c1 and c2 vertices. It is double degenerate due to the breaking of Z2 -translational symmetry. In this case microscopic deviations from the ground state are possible. There is a ?nite correlation length in the system and local correlation functions decay exponentially. For details about phase transitions, magnetization, and the antiferroelectric phase etc. see [10] and [3].

3.2
3.2.1

The structure of a random state
Free fermionic point

This is the case when ? = 0. It is convenient to parameterize weights in this case as a = ρ cos u, b = ρ sin u, c = ρ. √ When a = b = 1/ 2 this model is equivalent to the domino tiling of the Aztec diamond. The limit shape was computed analytically in [2] and is a circle. 10

Figure 7: Antisegnetoelectric phase - in this phase, zig-zag paths alternate with zig-zags formed by empty edges

The height functions of the average states for several values of the parameter u are shown on ?g. 8. For ? = 0, the limit shapes can be computed explicitly using methods of [11] and they are ellipses, which agrees with ?g. 8.

2a = b

a=b Figure 8: Free fermionic point with ? = 0

a = 2b

3.2.2

Ordered phases

In phase I the a-vertices dominate and the Gibbs state in this case is given by the lowest height function ?g. 3. In phase II the b vertices dominate and according to [10] we should expect that the average state will be the state with the domination of b-vertices. In other words the average state in this case is given by the highest height function ?g. 2. 3.2.3 Disordered phase

In this case it is convenient to use the following parametrization of weights: a = r sin(γ ? u), b = r sin u, c = r sin γ 11

with 0 < γ < π , 0 < u < γ , and r > 0. In this parametrization ? = cos γ . Phase III contains the free fermionic curve ? = 0. Since all this phase is critical one may expect that the nature of the Gibbs states will be similar for all parameters a, b, c in this region. In particular, one can expect the existence of the limit shape as in the case ? = 0. The particular form of the limit shape may vary but the following common features should common for all values of a, b, c in this region: ? The limit shape is a smooth curve having exactly one common point with each side of the square. At this point the limit shape is tangent to the side of the square. ? Inside of the boundary of the limit shape the height function is a smooth function and it has continuous ?rst derivative at the boundary. The second derivative has a discontinuity in the normal direction to the boundary of the limit shape. ? Outside of the boundary of the limit shape the height function is linear. Examples of Gibbs states in the disordered phase are shown on ?g. 9, 10, 11.

2a = b

a=b Figure 9: Disordered phase with γ =
π 4

a = 2b

3.2.4

The antiferroelectric phase

This region c > a + b is the one which is non-critical and which is also not ordered. In the periodic case the ground state has the domination of c-vertices as it is shown on ?g. 12. It is convenient to use the parameterization a = r sinh(η ? u), b = r sinh u, c = r sinh η with 0 < u < η . In this parameterization ? = ? cosh η . 12

2a = b

a=b Figure 10: Disordered phase with γ =
π 5

a = 2b

2a = b

a=b Figure 11: Disordered phase with γ =
π 8

a = 2b

In the case of DW boundary conditions there is a competition between very rigid restrictions on the states near the boundary which allows only a and b vertices near the boundary and the tendency of the system to have as much as possible of c vertices. Numerical simulations show that these competing tendencies resolve in the separation of three phases. It is fairly convincing from the ?g. 12 that the following should take place: ? The system forms a macroscopical droplet of the antiferroelectric phase with the boundary that does not touch the square. The height function in this domain is linear. The boundary of this domain has four cusps pointing towards sides of the square lattice. This phase is noncritical. Correlation functions in this region decay exponentially. ? Near the boundary the system is ordered. This ordered region is bordered by the disordered region where the height function is smooth. The disordered phase is critical. There is a ?nite magnetization, which means there are excitations with linear dispersions and therefore correlation functions 13

decay according to a power law. The boundary between ordered and disordered phases is a smooth curve with the features similar to the |?| < 1 case.

2a = b

a=b

a = 2b

Figure 12: Antiferroelectric phase with ? = ?3

2a = b

a=b

a = 2b

Figure 13: c-density plots of the antiferroelectric phase with ? = ?3

4

Conclusion

We demonstrated that local Markov sampling for the 6-vertex model with domain wall boundary conditions indicates that the system develops a macroscopical droplet of c-vertices when ? < ?1. For these computations it is not essential that the ground state of the 6-vertex model in this phase is doubly degenerate. This degeneracy corresponds to the translation by one step in the North-East direction on ?g. 12. This degeneracy is important in the computation of correlation functions and other observables. The two ground states correspond to the two parts of the graph of states which are connected by a ”very narrow neck” in the limit 14

N → ∞. We will address this problem in the next publication. The existence of the droplet can be seen from results of [14] where a di?erent numerical method was used. It would be interesting to compare the methods. The droplet of c-vertices is similar to the facets in dimer models. The shape of the droplet and of the surrounding critical phase is similar to the corresponding shapes in the dimer model on the square-octagon lattice [11]. The local Markov sampling which we used here is equally e?ective for other boundary conditions in the 6-vertex model. Some of the results for more complicated boundary conditions can be found in [4].

A
A.1

Functions and Implementations
Main loop

The following tasks must be completed by the Main loop function: 1. Import the matrix from a text ?le 2. Build ?ippables list 3. Set weights 4. De?ne ρ: ρ= 1 max{weight combinations for all ?ip types } (9)

5. Loop the following actions, and after a certain de?ned number of successful ?ips, output a ?le with the current matrix (and status of the Markov Chain) in it. (a) Generate a random real, rand, between 0 and 1 i. If
#?ip-up-only + #?ip-down-only + #bi-?ip #vertices

ii. Otherwise, go to (a).

≥ rand, continue to (b).

(b) Select a random ?ippable position with probability 1 P (selection) = #?ip-up-only + #?ip-down-only + #bi-?ip by calling the Get Flippable Position function. (c) Call Get Weight (which is now scaled by the value of ρ, to ensure that it always returns a value less than 1) to get the probability of an up ?ip and/or a down?ip at the ?ippable location chosen. (d) Generate a random real, rand, between 0 and 1. i. For up or down-only ?ips, i? W (Sb ) ≥ rand, execute the ?ip by calling the Execute Flip function, else restart main loop.

15

ii. For positions that can ?ip up and down, i? W (Sb ) ≥ rand, execute the ?ip corresponding to Sb , else i? W (Sb′ ) ≥ rand, execute the ?ip corresponding the Sb′ , else restart main loop. In practice, this means that once a vertex which can be ?ipped either way is chosen, simply divide up the probabilities of each ?ip occurring as discussed earlier.

A.2

Execute Flip
(a) Change the entry in the list of Flippables for the vertex chosen to make a high ?ip impossible. (b) De?ne the following positions: i. ii. iii. iv. v. vi. vii. viii. ix. x. One = the original position = Base Two = (+1, +0) = Right Three = (+1, +1) = Up Right Four = (+0, +1) = Up Left = (-1, +0) Down=(+0, -1) UpLeft = (-1, +1) UpRight = (+1, +1) DownLeft = (-1, -1) DownRight = (+1, -1)

1. If type is high

(c) Replace 4 parts i. Set Contents of Position One = FlipToOne (Position One, High) ii. Set Contents of Position Two = FlipToTwo (Position Two, High) iii. Set Contents of Position Three = FlipToThree (Position Three, High) iv. Set Contents of Position Four = FlipToFour (Position Four, High) (d) If Up, Down, Right, or Left Positions become ?ippable (call Get Is Flippable on each to check), add them to the Flippables List. (e) Call Fix Low End. 2. If type is low (a) Change the entry in the list of Flippables for the entry chosen to make a low ?ip impossible. (b) De?ne the following positions: i. One = (-1, -1) = Down Left ii. Two = (+0, -1) = Down 16

iii. iv. v. vi. vii. viii. ix. x.

Three = the original position = Base Four = (-1, +0) = Left Right = (+1, +0) Up = (+0, +1) UpLeft = (-1, +1) UpRight = (+1, +1) DownLeft = (-1, -1) DownRight = (+1, -1)

(c) Replace 4 parts i. Set Contents of Position One = FlipToOne (Position One, Low) ii. Set Contents of Position Two = FlipToTwo (Position Two, Low) iii. Set Contents of Position Three = FlipToThree (Position Three, Low) iv. Set Contents of Position Four = FlipToFour (Position Four, Low) (d) If Up, Down, Right, or Left Positions become ?ippable (call Get Is Flippable on each to check), add them to the Flippables List. (e) Call Fix High End. A.2.1 Fix High End

1. De?ne the following positions: (a) HighCreateDownLeft = (-1, -1) (b) HighDeleteDown = (+0, -1) (c) HighDeleteDownDownLeft = (-1, -2) (d) HighDeleteDownRighLeft = (-2, -1) (e) HighDeleteLeft = (-1, +0) 2. If HighCreateDownLeft is ?ippable, add it to the High Flippables List. 3. Delete 4 potential ?ippables on the High End. 4. If any of HighDeleteDown, HighDeleteDownDownLeft, HighDeleteDownLeftLeft, or HighDeleteLeft exists in the High Flippables list, remove them from the list. A.2.2 Fix Low End

1. De?ne the following positions: (a) LowCreateUpRight = (+1, +1) (b) LowDeleteUp = (+0, +1) (c) LowDeleteUpUpRight = (+1, +2) 17

(d) LowDeleteUpRighRight = (+2, +1) (e) LowDeleteRight = (+1, +0) 2. If LowCreateUpRight is ?ippable, add it to the Low Flippables List. 3. Delete 4 potential ?ippables on the Low End. 4. If any of LowDeleteUp, LowDeleteUpUpRight, LowDeleteUpRightRight, or LowDeleteRight exists in the Low Flippables list, remove them from the list.

A.3

Get Weight

1. Get contents of four surrounding positions if the ?ip occurred. 2. Multiply weights together corresponding to the contents of the four positions. 3. Multiply (new weight con?guration product) * ρ

A.4

Flip To

FlipTo functions take a position and a type and return what the vertex at the given position would be after the ?ip of the type speci?ed. 1. FlipToOne (Position, Type) (a) If type is high: If vertex was a1 , it will be c1 ; if it was c2 , it will be a2 (b) If type is low: If vertex was c1 , it will be a1 ; if it was a2 , it will be c2 2. FlipToTwo (a) If type is high: If vertex was b2 , it will be c2 ; if it was c1 , it will be b1 (b) If type is low: If vertex was c2 , it will be b2 ; if it was b1 , it will be c1 3. FlipToThree (a) If type is high: If vertex was a2 , it will be c1 ; if it was c2 , it will be a1 (b) If type is low: If vertex was c1 , it will be a2 ; if it was a1 , it will be c2 4. FlipToFour (a) If type is high: If vertex was b1 , it will be c2 ; if it was c1 , it will be b2 18

(b) If type is low: If vertex was c2 , it will be b1 ; if it was b2 , it will be c1

A.5

Get Flip Position

1. Generate a random integer between 1 and the total number of ?ippable positions; that is, the number of up-?ip only plus the number of down-?ip only plus the number of bi-?ips. 2. Choose the corresponding element in the Flippable Positions list to the random number chosen.

A.6

Get Is Flippable

Get Is Flippable should check the status of a position to determine if it is ?ippable. 1. High ?ippables must be a1 or c2 vertices and must have empty upper right corners (upper right corner must be a2 or c2 ). High ?ippable positions must have an x axis coordinate that is less than or equal to the width of the matrix 1 (where 0,0 is the origin) and a y axis coordinate that is less than or equal to the height of the matrix 1. 2. Low ?ippables must be a1 or c1 vertices and must have empty lower left corners (lower left corner must be a2 or c1 ). Low ?ippable positions must have an x axis coordinate that is no less than 1 (where 0,0 is the origin), and a y axis coordinate that is no less than 1.

19

B

Images of the N = 1000 matrix

Figure 14: N = 1000 plot for the antiferroelectric phase with ? = 3 and 2a = b

20

Figure 15: N = 1000 c-vertex density plot for the antiferroelectric phase with ? = 3 and 2a = b

References
[1] H. Cohn, M. Larsen, J. Propp, The shape of a typical boxed plane partition. New York Journal of Mathematics 4 (1998), 137-165 [2] H. Cohn, N. Elkis, J. Propp, Local statistics of random domino tilings of the Aztec diamons, Duke Math. J., 85 (1996), 117-166. [3] R.J.Baxter, Exactly solved models in statistical mechanics Academic Press, San Diego, 1982. 21

[4] K. Eloranta, Diamond Ice, J. Statist. Phys. 96 (1999), no. 5-6, 1091–1109. [5] A. Izergin, Partition function of the 6-vertex model in a ?nite volume, Sov. Phys. Dokl., 32(1987),878-879. [6] R. Kenyon, An introduction to the dimer models.math.PR/0308284 [7] V.Korepin and P. Zinn-Justin, Thermodynamic limit of the Six-Vertex Model with Domain Wall Boundary Conditions.cond-mat/0004250. [8] Korepin, V. E. Calculation of norms of Bethe wave functions. Comm. Math. Phys. 86 (1982), no. 3, 391–418. [9] G. Kuperberg, Another proof of the alternating-sign matrix conjecture, Int. Math. Res. Notes, 3 (1996),139-150. [10] E.Lieb, F.Wu, Two dimensional ferroelectric models, In: Phase transitions and critical phenomena, ed. by C. Domb, and M.S. Green, Academic Press, 1972. [11] R. Kenyon, A. Okounkov, Amoebae.math-ph/0311005. and S.She?eld Dimers and

[12] J. Propp and D. Wilson, Coupling from the past: a user’s guide. Microsurveys in Discrete Probability, (Princeton, NJ, 1997); DIMACS Ser. Discrete Math. Theoret. Comp. Sci., 41, 181-192, AMS, 1998 [13] A. Sinclair, Algorithms for Random Generation and Counting Birkhauser, Boston, 1993. [14] O. F. Syljuasen, M. B. Zvonarev, Directed-loop Monte Carlo simulations of vertex models. cond-mat/0401491. [15] A.V. Razumov, Yu. Stroganov, Combinatorial structure of the ground state of O(1) loop model. math.CO/0104216. [16] Asymptotic combinatorics with applications to mathematical physics (ed. by A. M. Vershik), Springer Lecture Notes in Math.,1815, (2003). [17] Zinn-Justin, P. Six-vertex model with domain wall boundary conditions and one-matrix model. Phys. Rev. E (3) 62 (2000), no. 3, part A, 3411–3418. [18] J.-B. Zuber, On the Counting of Fully Packed Loop Con?gurations. Some new conjectures.math-ph/0309057.

22