Classical height models with topological order
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] M a r Classical height models with topological order
Christopher L. Henley
Dept. of Physics, Cornell University, Ithaca, NY 14853-2501, USA
I discuss a family of statistical-mechanics models in which (some classes of) elements of a finite group G occupy the (directed) edges of a lattice; the product around any plaquette is constrained to be the group identity e . Such a model may possess topological order, i.e. its equilibrium ensemble has distinct, symmetry-relatedthermodynamic components that cannot be distinguished by any local order parameter. In particular, if G is anon-abelian group, the topological order may be non-abelian. Criteria are given for the viability of particularmodels, in particular for Monte Carlo updates. PACS numbers: 75.10.Hk,05.50.+q, 2.20.Hj
I. INTRODUCTION “Topological order” [1, 2] in a system means it has an emer-gent ground state degeneracy (in the thermodynamic limit),but (in contrast to symmetry-breaking), no local order param-eter operator can distinguish the states. Topological order hasattracted great interest over the last 20 years, since (i) it can-not (by definition) be captured by the Landau order-parameterparadigm and is hence exotic from the viewpoint of traditionalsolid-state theory; [2]; (ii) it is associated with “fractional-ized” excitations; (iii) it is proposed to implement qubits bythe ground-state degeneracy, the coherence of which is robustagainst environmental perturbations [3, 4]; (iv) the formula-tion in terms of ground-state degeneracy makes it attractivefor numerical exploration by exact diagonalization [5].The best-known examples of topological order arequantum-mechanical: the quantum Hall fluids) and latticemodels based on Z , the simplest group, such as Kitaev’s toriccode [4]. Indeed, Wen [1, 2] once emphasized quantum me-chanics as a defining property of topological order. But we canseparate these notions: topological order (as defined above) ismeaningful in a purely classical model (as developed in thispaper) or in a quantum-mechanical model at T > [9] (so thate.g. its renormalization-group fixed points represent classicalbehaviors). Indeed, I would suggest that the subject of topo-logical order skipped over more elementary examples, owingto historical accident. Compare with the history of sponta-neous symmetry breaking: theorists understood classical crit-icality long before quantum criticality [10], and we approachquantum critical properties in light of their similarities or dif-ferences from the classical case.Analogously, it is hoped that, in the case of topological or-der, classical models will (at the least) be a pedagogical aid,and that behaviors evidenced in classical models may providea framework for conjectures about the quantum models. Dis-entangling classical notions and inherently quantum mechan-ical ones might lead to clearer (or at least different) thinking.Also, the framework in this paper naturally draws us to facehitherto unfamiliar groups – e.g. the group A (see Sec. XX)– and it might inspire the construction of quantum-mechanicalmodels involving these groups.The explicit notion of “classical topological order” was in-troduced and highlighted in [9], in particular the points that(i) it is characterized by ergodicity breaking, (ii) can be im- plemented by hard constraints, and (iii) must have a discreteclassical dynamics – all of which applies to the models in thispaper. However, much of Ref. [9] was framed in terms ofthe relation of the classical model to a quantum model, e.g.by taking the quantum model to a temperature at which quan-tum coherence is no longer important, [40] or by removingsome Hamiltonian terms. In the present work, the model isformulated from the start as an ensemble of classical statisti-cal mechanics, without concern for the existence of a quantumcounterpart. That will permit consideration of a richer set ofmodels (i.e. discrete non-Abelian groups G ), for which wemight not know how to concoct a good quantized version.I will consider models based on either abelian or non-abelian discrete groups. It should be noted that abeliannessin this paper has a different significance than in the quantumcontext. In the latter case, the group in question is the Berryphase (or its generalization to a unitary matrix) induced byevolution of the wavefunction from one state to an equivalentone. A specific case is the statistics of quasiparticles whoseworld lines braid around each other. Quantum-mechanicalnon-abelianness may be realized in models based on abelian groups such as Z . Most of the fractional quantum Hall flu-ids are abelian, but non-abelian statistics is more suitable forquantum computation [4, 6]. Proposed realizations of non-abelian topological order in this sense are formulated in latticemodels as sums over loop coverings [7, 8].The primary focus here is not on ideas that point to analyticsolutions or to connections with the existing literature of topo-logical order. Instead, this is meant as a generic blueprint fornumerical studies. For example, the classification of modelsin Sec. III (as summarized in the tables) is motivated by theneed to select a good one for simulations, and the quantitiesdefined in Sec. V are all measurable in Monte Carlo simla-tions. However, the actual simulation results will be left tosubsequent papers [11]. A. Height models
I shall realize topological order by generalizing the conceptof “height models”. Their defining property [12–20] is theexistence of a mapping directly from any allowed spin con-figuration { σ i } to a configuration of heights h ( r ) , wherein h ( r ) − h ( r ′ ) , for adjacent sites r and r ′ , is a function of thespin variables in the neighborhood (normally, either two spinson the sites r and r ′ , or else on one spin on site i at the mid-point of the r - r ′ bond: spins and heights commonly live ondifferent lattices).The “spins” in the model could be any discrete degree offreedom – e.g. dimer coverings. For { h ( r ) } to be well de-fined, it is necessary (and sufficient) that the sum of heightdifferences is zero round any allowed plaquette configuration.Thus, a (local) spin-constraint is assumed that excludes (atleast) the configurations without well-defined heights, yet stillallows a nonzero entropy of states S in the thermodynamiclimit. In the simplest cases, h ( r ) is integer-valued.Such models may have “rough” phases, in which thecoarse-grained h ( r ) behaves as a Gaussian free field, i.e. theeffective free energy of long-wavelength gradients is F = 12 Z d r |∇ h | . (1.1)Via the apparatus of the Coulomb gas formalism [21, 22], thisimplies the spin variables have power-law correlations (criti-cal state); the topological defects may have unbinding tran-sitions like the Kosterlitz-Thouless transition. Indeed, in twodimensions most critical states can be addressed as “heightmodels” [21] and this formalism provides an alternative routeto computing exact critical exponents [17], besides conformalfield theory.Note that on a topologically non-trivial space (such as thetorus periodic boundary conditions), there are nontrivial loops ℓ , such that the net height difference (or “winding number”) w ℓ added around such a loop is nonzero. It is easy to see thisis a topological invariant, in that w ℓ is unchanged if the loop isshifted and deformed (so long as it stays topologically equiva-lent to the original one.) Thus, if ℓ , ℓ , ... are the fundamentalloops, the configuration space divides up into sectors labeledby ( w ℓ , w ℓ , ... ) . Here, and also for the discrete-group heightmodels introduced in the paper, a sector is each set of config-urations which can be connected to each other by a successionof local spin changes (i.e. “updates” in the terminology usedlater).Point defects may also be admitted and loops around themmay also have a nontrivial w ℓ (in which case they are topo-logical defects). Clearly, the winding number w ℓ in a heightmodel is analogous to t ∈ T in a topologically ordered model;we could almost say this is a special case of topological orderin which T is Z , here meaning the (infinite discrete) group ofintegers under addition.In place of a Landau order parameter, the (near) degeneratestates in a topologically ordered system may instead be distin-guished by a global loop operator, acting around a topologi-cally nontrivial loop ℓ . Just as a Landau order parameter formsa group representation of a broken symmetry, and labels thesymmetry-broken states, in topological order the global loopoperator ought to form a faithful representation of the “topo-logical group T ” whose elements label the distinct states. Inour models, the definition of this loop operator is trivial andtransparent: it is just the generalized “height difference”. B. Outline of the paper
In this paper, I first (Sec. II) generalize the height-modelidea to the case where the height variable belongs to a discrete(finite) abelian or non-abelian group, thus defining a familyof classical models which (in many cases) has a topologicalorder. The models are defined by a lattice, a group, and theselected subset of group elements which are permitted valuesfor the ‘spins” of the model; the spins sit on the bonds. Non-Abelianness of the group has interesting consequences: forexample, a collection of defects no longer has a unique netcharge (Sec. II C)In Sec. III and Sec. IV, I survey the various combinationsfor the smallest non-Abelian groups, using the crude Paulingapproximation as a figure of merit to identify the most attrac-tive models for Monte Carlo simulation, using single-site up-dates. Furthermore, Sec. V suggests what quantities are inter-esting to measure in such a simulation; however, no simula-tion results are reported in this paper. But some first analyticresults are included in Sec. VI, based on transfer matrices andhence implicitly one dimensional: the main point is to shedlight on how the size dependence or defect pair correlationdepends on the group elements labeling the topological sectoror the defects.Finally, the conclusion (Sec. VII) reflects on which topolog-ical behaviors are inherently quantum mechanical, and whichare not (in that the same behavior can be found in classicalmodels). Furthermore, applications are suggested, either tosimulating systems with vacancy disorder, or to constructingquantum versions of the models in this family.
II. DEFINITIONS AND TOPOLOGICAL BEHAVIORS
This paper is meant to introduce (and compare) a wholefamily of models. In this section, I define the general rules forthis family, and then describe the most promising examples.(In the next section, I shall exhibit consequences for MonteCarlo updating of such models.)
A. Model definition
Let us take a “lattice” (not necessarily Bravais, e.g. honey-comb) of sites r . The spins sit on the bonds of this lattice, andtake values in the discrete group G ; they are σ ( r , r ′ ) ∈ G (2.1)where ( r , r ′ ) labels a bond of nearest neighbors. The bondsare directed; if we reverse the direction we view a bond, thespin on it turns into its inverse: σ ( r ′ , r ) ≡ σ ( r , r ′ ) − . (2.2)Then each configuration of the spins induces a configurationof “heights” h ( r ) ∈ G , defined by h ( r ) = σ ( r , r ′ ) ∗ h ( r ′ ) , (2.3)where “ ∗ ” represents the group multiplication. Of course, h ( r ) is only defined modulo a global multiplication by someelement τ , h ′ ( r ) = h ( r ) τ ; to make it well-defined, we couldarbitarily require (say) h (0) ≡ e , where e is the group iden-tity. One could then explicitly construct h ( r ) at the neighborsof site , and iteratively at their neighbors, etc.; the result isindependent of which bonds ( r , r ′ ) are used for this, if andonly if a plaquette constraint is satisfied [Eq. (2.5), below].It will be useful throughout to define the line or loop prod-uct of p spins: γ ( ℓ ) ≡ σ ( r , r p ) ∗ σ ( r p , r p − ) ∗ ... ∗ σ ( r , r ) (2.4)Here the loop ℓ is a string of bonds ( r k , r k +1 ) connecting end-to-end, for k = 0 , ...p ; it is a loop when r p = r .
1. Plaquette constraint
Two constraints are imposed on the spin configurations.The first is the plaquette constraint : we require the loop prod-uct around any elementary plaquette, to be the identity, γ ( ℓ plaq ) = e (2.5)This is necessary (and sufficient) for h ( r ) to be well-defined.One can define variants of any model by relaxing the pla-quette constraint to allow a small number of defect plaquettes,around which the loop product is not the identity. Supposethat defects cost an energy ∆ (possibly depending on the kindof defect): then the basic, defect-free version of the model canbe viewed as a Boltzmann ensemble in the limit T / ∆ → .On the other hand, if we imagine there were a spin-spin inter-action that breaks the degeneracy of the states satisfying theplaquette constraint, then the basic version of the model is the T /J → ∞ limit.Using condition (2.5) and induction (adding one plaquetteat a time to the loop), the loop product must be γ ( ℓ ) = e forany finite loop ℓ that is contractible to a point in small steps.
2. Spin constraint
The second constraint is the spin constraint : choose a “spinsubset”
S ⊂ G such that σ ( r , r ′ ) ∈ S (2.6)everywhere. Hence, the choice of S is a major part of amodel’s definition; in Sec. III C, below, I will discuss otherdesirable features of S . The spin constraint is retained even inversions of the model with defect plaquettes.The constraints should respect the group and lattice sym-metries . Implementing the group symmetry means requiring σ ∈ S ⇒ u ∗ σ ∗ u − ∈ S (2.7)for any conjugating element u . Thus, S must be one of thegroup’s conjugacy classes – the simplest case – or a union of such classes. (Some non-abelian groups, and all abelianones, have “outer” automorphisms, symmetries which cannotimplemented by conjugacy within G : we may also wish toimplement those symmetries, too).To implement lattice symmetry, one asks σ ∈ S ⇒ σ − ∈ S (2.8)so the model respects inversion (around the bond’s midpoint).[41]Let us further require e / ∈ S (2.9)thus a uniform height configuration is disallowed. Finally, andtrivially, S generates the full group G . (2.10)(If not, I could have redefined G as the subgroup generated by S .)Without the spin constraint, the models would be identicalto the lattice gauge models of Douc¸ot and Ioffe [23]. In such amodel, only gauge-invariant (i.e. loop) quantities can havenonzero expectations; other correlation functions are zero,even at the nearest distance. In contrast, these group-heightmodels (like the original height models) have nontrivial finite-size effects and local correlations: in particular, there are me-diated interactions between topological defects.Furthermore, the spin constraint allows the possibility of along-range ordered phase, particularly if we assign differentBoltzmann weights to different configurations, and we mayfind phase transitions as those parameters aare varied. Partialorderings are also possible, and transitions might occur be-tween different topological orders It will be easier to explorethe phenomenology of such transitions in the classical realm. B. Topological sectors and topological order
For these models, a “sector” means simply the configura-tions that can be accessed by a succession of local updates.(Here “local update” means an operation that turns one validconfiguration to another by changing spins in a small neigh-borhood of some site, as might be deployed for Monte Carlosimulation. A “nonlocal” rearrangement, as developed inSec. IV, means the cluster of updated sites can be arbitrarilylarge, and in particular could include a topologically nontriv-ial chain of sites that spans the periodic boundary conditions.)By this definition, “sectors” are well defined in any finite sys-tem larger than the maximum update cluster. In other modelsand with other definitions, passing between sectors is be abso-lutely forbidden, so that sectors (more exactly “components”,like the up and down ordered phases in ferromagnet) are emer-gent in the thermodynamic limit.In these models, sectors can be labeled by loop products γ ( ℓ ) . As noted after (2.5), products around topologically triv-ial loops must give the identity, but others – e.g. through theperiodic boundary conditions of a torus system – in general donot. Such loop products are not changed by local updates, andtherefore must take the same value for all states in a sector.So we call sectors “topological” when they are distinguished(and necessarily disconnected) by having different values ofthe loop product(s). The definition of topological order is that– in the thermodynamic limit – different topological sectorsall become equivalent, in that they cannot be distinguished byany local expectations. Furthermore, just as the ground stateenergies in different sectors should become equal in the caseof quantum topological order, the free energies should becomeequal in our models.Let ℓ , ℓ , ...ℓ g be the basic independent loops, where g isthe genus; then the loop products { γ i ≡ γ ( ℓ i ) } [all takenfrom the same origin] label the possible topological sectorsof configuration space. An interesting question is how manydistinct sectors there are, [24] given the system’s genus g .
1. Invariance
Before counting sectors, we need to explore invarianceproperties of the sector labels, in case some labels are equiv-alent to others. First, there is a sort of gauge freedom: if wehad evaluated these loops starting from r instead of from 0,then γ ℓ → γ ′ ℓ = γ r ∗ γ ℓ ∗ γ − r , (2.11)where γ r is the line product along any path from the originto r ; notice that this same element conjugates all the distinctloops. However, just because our labeling fails to distinguishtwo sectors does not conclusively show they are the same.A better criterion for counting sectors as equivalent isthat one can be turned into another by local updates. Keepthe same origin, but perform a single-site update [seeEq. eqreq:sigma-new-inner, below] hitting on the origin ver-tex, one gets another conjugacy γ ℓ → γ ′ ℓ = τ ∗ γ ℓ ∗ τ − , (2.12)where τ is now the updating multiplier.When the group G is abelian the sector labels are invariantwith respect to how we take the loop and unchanged by localupdates. We can have an independent and invariant loop prod-uct γ i for every topologically independent loop, so the numberof sectors is n G g where g is the system’s genus ( g = 2 fortorus), and n G is the number of elements in the group.
2. Sector counting in non-abelian case
On the other hand, in the non-abelian case, a loop productis invariant only up to a conjugacy, so we have fewer sectors.Furthermore, the allowed values of distinct loop products arenot independent. Consider a square lattice model in a rect-angular system cell L x × L y , with periodic boundary condi-tions; let ( γ x , γ y ) be the loop products along straight lines ofbonds running from the origin site (0 , , in the x or y direc-tions respectively. The loop running from (0 , to ( L x , to γ γ γ γ γ γ γ γ γ γ γ γ a) b) c)P P P FIG. 1: The composite defect charge of a pair may be changed bysliding a third defect between the two. (a). Two defects (stars) havecharges γ and γ , thus a loop enclosing both contains charge γ γ .A third defect with charge γ is being moved from the right. The loopproduct around defects 3 and 2 is γ γ These charges are definedusing the reference point P , and multiplications are from the right.(b). After defect 3 is moved between defects 1 and 2, the productaround defects 2 and 3 is still γ γ . (c). After defect 3 has passedto the other side, the product around defect 2 has the new value γ ′ ;the product around 2 and 3 is γ γ ′ , which is unchanged from (b) andtherefore equal to γ γ . Hence γ ′ = γ γ γ − , and the combinedcharge of defects 1 and 2 is now γ ′ γ = γ γ γ − γ , which can bein a different class than γ γ . ( L x , L y ) to (0 , L y ) and back to (0 , contains no defect, soby inductive use of the plaquette contraint its loop product is e .Yet the four segments of this loop are just γ x and γ y , forwardsor backwards, so the loop product is γ − y ∗ γ − x ∗ γ y ∗ γ x = e ; (2.13)in other words, γ x and γ y must commute. So, in effect, we must define an equivalence relation ( γ x , γ y ) ∼ ( γ ′ x , γ ′ y ) whenever the pair satisfies (2.13), andeach topological sector is one equivalence class. In the caseof a larger genus, we extend in the obvious way to longer lists;Some obvious kinds of equivalence classes are:(i) ( e, e ) (ii) ( ω, e ) or ( e, ω ) (iii) ( ω, ω k ) or ( ω k , ω ) for k = 1 , ..., m − , where ω is anelement (not the identity) of order m .(iv) If the group has a nontrivial center G Z , consisting ofelements that commute with all the other elements, then if ( γ x , γ y ) is a sector then ( z x γ x , z y γ y ) is another sector, where z x , z y ∈ G Z .Table I shows the number of classes µ and the sector count µ for some groups of interest. C. Defects and non-abelian effects
We can allow the possibility (dilutely) of a plaquette thatviolates the plaquette constraint. The loop product around itwill be called β . It is analogous to the Burgers vector of adislocation, or of the topological defects in the usual heightmodels based on Z ).In the case of a height model (in a rough phase), topolog-ical defects are like vortices in a two-dimensional Coulombgas[21, 22]. They behave like U (1) electric charges in a two-dimensional universe. Opposite charges feel an attractive log-arithmic potential. In contrast, in the case of topological order,the attractive potential decays exponentially and the defectsare deconfined [27].If the topological order is non-abelian, the non-commutingproperty of defect charges has some interesting consequences.They are not unique to topological order; this also is a longknown property of defects of traditional ordered states asso-ciated with a non-abelian homotopy group [28]. One conse-quence is that a loop product may be changed when a defectof charge β is passed across it, the action being a conjugation: γ ( ℓ ) → βγ ( ℓ ) β − . (2.14)Thus the topological sector might be changed when a defectwanders around the periodic boundary conditions. Also, thenet charge of a defect pair can be changed by passing anotherdefect between the pair. (See Fig. 1)Another consequence of non-abelianness is that given twogiven defects of specified charges, there is more than onepossible value for their combined charge. In the quantum-mechanical approaches to defects in topologically orderedsystems, this same property is also the hallmark of non-abelianness. In that context, the list of allowed combinationsis known as the “fusion rules”, and there are matrices (gener-alizations of Clebsch-Gordan coefficients) which tell how toform the appropriate linear combinations.A third consequence is pertinent to simulations and the def-inition of topological sectors in the presence of defects. Ifone creates a defect pair and moves one defect around theboundary conditions, it may recombine with the original de-fect into a single defect, rather than annihilate. The fact thattwo defects may not be able to re-annihilate is very similar tothe “blocking” idea of Ref. [24] (for quasiparticles in a non-Abelian quantum Hall state).The single defect state satisfies the generalization of (2.13),namely γ − y ∗ γ − x ∗ γ y ∗ γ x = β. (2.15)This commutation is a “group commutator”. Such a single-defect state may conveniently allow numerical measurementsof the creation free energy of a single defect. Of course, sucha state is never possible for abelian defects; in that case, asystem with periodic boundary conditions must have eitherno defects, or at least two of them. III. POSSIBLE MODELS
In this section, I survey specific models, emphasizing thecriteria which would make some of them particularly attrac-tive for future investigations. To summarize Sec. II A: modelsin this paper are specified by (i) the group G (values of height)(ii) the spin subset S (values of spins) (iii) the lattice whosebonds the spins sit on.Therefore, the models will be named in the form“ G ( m ) latt”. Here “ G ” is the groups name, ( m ) is the or-der of the elements in the selected conjugacy class (usuallythat is unambiguous), and “latt” abbreviates the lattice (“tri”, a a aaaab b bbb bb bc ccc c ccb b b ccc aaa b a γ y aproduct γ x b a b bb a aa b cc c a ac ccb b babac a bc b c ca b ccca a (a) (b) FIG. 2: [Color online.] (a) Example configuration of the (abeliangroup) model Z × Z (2) sq . Each square lattice edge is occupiedby a group element a , b , or c ; directions are unneeded since each ofthese is its own inverse. The loop products γ x and γ y around the pe-riodic boundary conditions are also shown, which define the topolog-ical sector. For an abelian group, their values are independent of thestarting points. (b) Example configuration of the (non-abelian group)model S (2 , tri. Labels a, b, c denote the elements (23),(13),(12),while the arrow denotes a cyclic exchange (132). “sq”, or “hc” for triangular, square, and honeycomb). Thus“ S (2 , tri” means that G is the permutations of three ob-jects, S contains all three pair exchanges, as well as the twocyclic permutations (i.e. every group element except for e ),and “tri” means the spins sit on the edges of the triangularlattice. An variant nomenclature is sometimes convenient, inwhich the “ ( m ) ” in the label gets replaced by “ { σ , σ , ... } ”:the set { σ , σ , ... } is simply the listing of the selected ele-ments. A. Groups
Table I lists the groups and spin subsets I shall be interestedin.For future reference, I mention the automorphism group A G of a group G , which is simply its symmetry group. Each a ∈ A G is a permutation of the group elements preserving itsstructure, a ( gg ′ ) = a ( g ) a ( g ′ ) . When G is non-abelian, thereis a subset of the automorphism group called the inner auto-morphisms , defined as the conjugations, a τ ( g ) ≡ τ gτ − . Ob-viously a τ a τ ′ = a τ ∗ τ ′ , so the inner automorphism subgroupis isomorphic to G / G Z , where G Z (the center subgroup) con-sists of the elements that commute with everything. But manygroups have additional outer automorphisms that are not con-jugations; in particular, all automorphisms of an abelian groupare outer.
1. Abelian groups
We start by considering discrete abelian groups in this fam-ily of models. The smallest of them, Z , does not work since S can have only one element. The next simplest cases arecyclic groups Z q , i.e. the integers modulo q under addi-tion, although these often turn out to be height models (seeSec. III B 2 below).Beyond that we go to direct products of cyclic groups, in-deed any abelian group can be represented thus. If S was also TABLE I: Groups and spin subsets. µ is the number of conjugacyclasses, and µ the number of topological sectors on a torus; n G , n S ,and n E respectively are the number of elements in the group G , thenumber in the selected subset S , and the number of even elements.The effective bond probability p b is given by formula (4.4).group+tag n G µ µ n S n E p b Z × Z (2) 4 2 16 3 – 1/3 S (2) 6 3 8 3 3 0 S (2,3) 5 – 1/5 Q (2) 8 3 10 6 – 2/7 D { m, m ′ } A (3) 12 3 8 8 4/11 A (2) 60 4 20 15 – 45/59 A (3) 20 – 40/59 A (5) 12 – 48/59 taken to be a direct product, of course the model would reduceto a superposition of non-interacting models, one for each fac-tor. However, there are many attractive examples in which S is not a direct product, in particular Z × Z (see Sec. III B 1).
2. Non-abelian groups
The smallest non-abelian group is S , the permutationson three objects (also isomorphic to the dihedral group D ).Here, S may be taken as the class of all pair permutations, themodel S (2) , or as all permutations except the identity, that is S (2 , in our notation.Each of the next two smallest non-abelian groups has eightelements. One of these is the 8-element quaternion group Q ,i.e. the unit elements {± , ± i , ± j , ± k } from the quaternionring. Here S must be the class of the six elements not equal to ± .The other eight-element non-abelian group is D , the sym-metry group of the square lattice.The “alternating groups” A and A are especially attrac-tive for our purposes due to their high symmetry (so we canchoose S to be a single class containing a sizeable fractionof all the group elements). They consist of the even permuta-tions of four and five elements. Note that A and A are alsothe point groups of the (proper) rotations of a regular tetra-hedron and a regular icosahedron, respectively. Being sub-groups of SO (3) , these groups might in some sense serve as adiscretization of it [29], just as clock models are a discretiza-tion of the XY model. That would be interesting as a way tomake a connection to topological models (or gauge theories)defined in terms of Lie (i.e. continuous) groups.Finally, A is the smallest non-abelian simple group, mean-ing it has no normal subgroups; as we shall see in a moment(Sec. III B), normal subgroups are an annoyance since theytend to make the behavior more trivial than would be expectedfor the group G . B. Example models
Next I shall survey the simplest examples. Most of themreduce, in some fashion, to previously known models; that isan advantage for computational studies, since old results canbe used as checks. In several cases, the models in our familyhave “accidental” topological order, i.e. beyond the group G ;In particular, some of them have height representations.The group and subgroup involved in our spin constraint arefinite, and so is each plaquette; thus it can happen that theallowed configurations satisfy stronger constraints than thosethey were designed to fulfill. The first five subheadings belowall, in one sense or another, reduce to known models. Z × Z and the 3-coloring model For a first example, let
G ∼ = Z × Z , an abelian group.Besides the identity, this group has three equivalent elements a , b , c ; each has order two, and the product of any two givesthe third. If we treat these as a class (although they are notconjugate, since the group is abelian), then we must choosethat class to be the spins, Z × Z (2) . Since a = a − , etc., wecan depict the spins using three (undirected) “colors” of theedges. On the square lattice this gives perhaps our simplestexample (Figure 2).What about the triangular lattice case [model Z × Z (2) tri ]? The plaquette constraint is simply that each trian-gle has one edge of each color: the “three-coloring model”. (Itis usually represented on the edges of the dual [honeycomb],where the constraint says each vertex has three colors; in ei-ther case, the spins live on kagome lattice vertices, and theconfigurations are also the ground states of the 3-state Pottsantiferromagnet on that lattice.) This model is known to havea Z × Z height representation [17, 26], in addition to the finite-group height field h ( r ) defined by (2.3).
2. 6-vertex model
For another example, take G to be Z q , with q > , and letthe lattice { r } be the square lattice. Choose S = { +1 , − } .(The two elements are not the same class; they are related onlyby an outer automorphism.) Then the sum of spins around aplaquette can be zero (mod q ) only if it is just zero, i.e. thereare exactly two +1 and two − in the loop. If we express thesespins on the dual (also square) lattice, as an arrow pointingoutwards (resp. inwards) wherever σ = +1 (resp. − ) asthe loop is traversed counterclockwise, we see these just arethe configurations of the six-vertex model – which also has ainteger-valued height field. [42]Since Z or Z m are abelian groups, { +1 , − } is merely an outer class.
3. Groups with an even subgroup
An even subgroup E (with n E ≡ n G / ) has G / E ∼ = Z .That is, any product of an even number of elements lies in E .Say that the spin subset S consists of odd elements. (If it con-sisted of even elements, we would generate at most E .) Noticethat such a model cannot use the triangular lattice, since theplaquette rule cannot be satisfied (the plaquette product mustbe odd, while e is an even element).Now if the simulation cell has even dimensions, the possi-ble topological products γ ( ℓ i ) must lie in E . (Even if the cellhas an odd dimension, the possible values of γ ( ℓ i ) still corre-spond 1-to-1 with elements of E .) Thus, the topological sectorlabels can only belong to the subgroup E .For example, in permutation groups the even subgroup con-sists of even permutations. In the case of S , the even permu-tations are just (123) and its powers, so E ∼ = Z . Conse-quently the model S (2) sq can have only abelian topologicalbehavior. In our list, D is another group that contains an evensubgroup (the proper rotations).
4. Groups with a center
What if G has a non-trivial center G Z ? (The center is sub-group of elements that commute with every other element).For example, if we adopt the group Q of unit axis quater-nions, (which order 8) then Q Z = { +1 , − } ∼ = Z and Q/Q Z ∼ = Z × Z . Thus, the model Q (2) tri projects ontoconfigurations of the the three-coloring model. (Just map ± i → a, ± j → b, ± k → c .) The group D also has a cen-ter (two-fold rotations, i.e. inversions, commute with every-thing.) C. Criteria for models: estimates of entropy
To estimate at once the viability of many different models, Ishall use very crude estimates of the entropy and (in Sec. IV B)updatability. Say the lattice has coordination z and the duallattice has coordination z d , i.e. the number of sides of eachplaquette. (These numbers are related by /z + 1 /z d = 1 / .)Also, say the group has n G elements, of which n S are in theselected subset S . These three parameters — n G , n S , and z (or z d ) – contain much of what we need to characterize thepossible models. See Table II for the parameters related tolattice geometry, and Table I for those related to the groupsand the spin subsets.I will use a Pauling estimate for the entropy. There are N z/ edges and hence n S Nz/ ways of placing spins inde-pendently chosen from S , in a hypothetical ensemble that doesnot (yet) enforce the plaquette constraint. If we knew the frac-tion of all these states that do obey the plaquette constraint, wewould have the total count of allowed states and thus the de-sired entropy.Pauling’s approximation is to pretend the event of satisfy-ing the plaquette constraint is uncorrelated between plaque-ttes. So let f e be the chance that a given plaquette has plaque- TABLE II: Lattices (asterisk denotes an average over two kinds ofplaquettes). Here “ σ -phase lattice” denotes the lattice (3 , , , and “square-octagon” lattice denotes (4 , ) . The columns give thecoordination numbers z of the lattice and z d of the dual lattice, fol-lowed by the lattice’s bond and site percolation thresholds p cb and p cs . (Many more significant digits are known [34].)lattice tag z z d p cb p cs triangular tri 6 3 0.347 0.500 σ -phase – 5 3.333 ∗ ∗ ∗ tte product equal to e . Then in this approximation, the prob-ability is f Nz/z d e that the plaquette constraint is satisfied onall N z/z d plaquettes of the whole system. Thus the Paulingestimate of the ensemble entropy is e NS Pauling = (cid:16) n S z/ f z/z d e (cid:17) N . (3.1)The condition we must satisfy, in order to have an ensembleat all, is S Pauling > , i.e. n S > /f /z d e ., (3.2)If z d is not too small, we may estimate that the plaquetteproduct is equally likely to be any group element, hence very crudely f e ≈ /n G . (3.3)Less crudely, one can work out the the actual probabilities thatthe product of z d random group elements from the allowed set S will give the identity, and these are the f e values in Table III.Comparison of those f e values with n G in Table I shows thatusually, (3.3) is not bad. When the product of two elementsof S is particularly likely to fall into one class, the true f e de-viates more from (3.3), either on the low side (e.g. the model A (2) tri) or the high side (e.g. A (5) sq). The extreme case isif the group contains even and odd elements, and S consists ofodd elements (the S (2) or D ( m, m ′ ) examples in Table III).In that case we should replace (3.3) by f e ≈ /n E = 2 /nG if z d is even, but f e = 0 if z d is odd.There is just one entry in Table I showing a negative Paul-ing entropy S Pauling < , namely A (2) tri. It is convenient toexplain this case in the language of (proper) rotation group ofan icosahedron, which is isomorphic to A . The only way thatthree twofold elements can multiply to give the identity is mu-tually when the two fold axes are mutually orthogonal. Sincethe triangles share edges, any valid global configuration mustuse that same triad in every triangle; this entails a fivefold( S ) global symmetry breaking, since the fifteen twofold axesof icosahedral symmetry break up into five disjoint orthogo-nal triads. Indeed the three used elements form a subgroupisomorphic to Z × Z so we are back at the three-coloringmodel for which S Pauling > . The Pauling approximationgave zero entropy only because it did not take account of thesymmetry-breaking and attempted to mix domains with in-compatible symmetry breakings.My purpose here is not to obtain quantitative estimates ofthe model’s entropy, although the Pauling estimate is some-times surprisingly accurate. Rather, I want to compare thesevalues between different models as a figure of merit, to aidus in guessing which models are the most interesting or themost tractable. To this end, the figures of merit are shown inTable III.To satisfy Eq. 3.2, the three parameters get pushed in thefollowing directions, but there are considerations limitingeach of the three.(1) We want n G as small as possible; however, there arenot so many small, discrete, non-abelian groups: onlythree have n G ≤ , namely S (permutations of threeobjects), Q (quaternion group), or D (point group of asquare).(2) We want larger n S , meaning the model is lessconstrained (and more tractable). In the limiting case n S = n G , the model is just a pure gauge theory, whichis trivial apart from its global topological properties.On the other hand, a sufficiently large n S requires in-cluding more than one conjugacy class in S , so that thespins can have inequivalent “flavors”. That is estheti-cally undesirable: a generic model (with unequal statis-tical weights) needs more parameters, and it is harder toimagine how such a model could be realized physically.(3) We want large z d , as in the honeycomb lattice.However, it is esthetically harder to implement a prod-uct constraint in a physical model. (When the prod-uct string is short, there are only a few symmetry-inequivalent cases for it, and it is easier to concoct aHamiltonian term which does not reference the groupmultiplication, but which has those cases as its energyminimum.)To satisfy (3.2) with a large group but S consisting of justone conjugacy class, the group must have high symmetry.E.g., the alternating group A (the proper icosahedral rota-tions) has n G = 60 and contains conjugacy classes with 12,15, or 20 elements, which using (2) would need z d > . , . , or . respectively. IV. MONTE CARLO UPDATING
For us, one essential criterion of a model is the possibil-ity of Monte Carlo simulation. I limit consideration to theequal-weighted ensemble, in which every allowed configura-tion has the same weight. Then detailed balance is satisfiedif the forwards and backwards rate constants are the same forany update move. But what is the minimum sufficient updatemove? For the six-vertex model it sufficed to reverse the ar-rows on the four edges of one plaquette, which changes the height field on one site, a purely local update. For the three-coloring model, the minimal update involves switching two“colors” (e.g. a ↔ b ) along a loop , a nonlocal update move.What happens generically for our family of models? A. Cluster update move
The update move is simplest described in terms of theheight function (defined in (2.3)). First pick at random a groupelement τ = e and a starting site r . Say D is the domain be-ing touched by the update. (It will be explained in a momentwhat determines D ). Then I prescribe that the update premul-tiplies the heights in this domain by τ , so as to “shift” them: h ′ ( r ) = ( τ ∗ h ( r ) for r ∈ D ; h ( r ) for r / ∈ D . (4.1)This induces the following update of the spin configura-tion: [43] σ ′ ( r ′ , r ) = τ ∗ σ ( r ′ , r )) ∗ τ − for r , r ′ ∈ D ; τ ∗ σ ( r ′ , r )) for r ′ ∈ D , r ′ / ∈ D ; σ ( r ′ , r )) ∗ τ − for r ′ / ∈ D , r ′ ∈ D ; σ ( r ′ , r ) for r , r ′ / ∈ D . (4.2)I call this a “gaugelike” transformation [33]: it has the sameform as a gauge transformation would, but it is valid onlywhen an additional spin constraint is satisfied too.If both endpoints of the bond are in D , then σ ′ is conjugateto σ and must be legal (since we include whole conjugacyclasses in S ).On the other hand, where the ( r , r ′ ) bond crosses the do-main boundary ∂ D , the spin constraint is nontrivial to satisfy.Let’s place an arrow along the edge from r to r ′ if and only if σ ( r , r ′ ) ∗ τ − / ∈ S . (4.3)In other words, there is an arrow from r to r ′ whenever in-cluding r in D forces us to include r ′ as well. This arrow isnot bidirectional. (it is in the case τ = e ). Thus, we mighthave any of four possibilities (no arrows, arrows both way, orarrows one way) along each bond.Then the update rule is to construct the arrowed-percolationcluster consisting of site r , with the rule that site r ′ is in-cluded if site r is included and there is an arrow r to r ′ . Thisis the smallest possible updated domain containing r . Ofcourse, we do not actually need to construct all the arrows; in-stead, we grow the cluster from the initial site, and constructarrows only from sites already in the cluster.Notice that (only) in the case G is abelian , (4.2) reducesto σ ′ = σ throughout the interior of D . In other words, theupdate only changes spins along the boundary ∂ D and thus isa loop update. In a non-abelian model, however, the update isgenerally a cluster update.In some models (see next subsection) there is a strongchance to hit a system-spanning cluster, including most of thesites, which tends to be inefficient. (Updating all the sites isequivalent to no update). To avoid this, a limiting size s max for the update cluster D should be set; if this limit is reached,we cancel the tentative move and start over, choosing a newrandom r and τ . B. Numerical criteria for cluster updates
Notice that in growing a cluster from r , we never caredabout the reverse arrows. Therefore, we obtain the same clus-ters as we would in an ordinary (not arrowed) percolationproblem, if the occupied bond probability p b is identified withthe probability of an arrow in a pre-selected direction; thatprobability is simply p b ≈ − n S − n G − , (4.4)if we chose the candidate updating factor τ at random. Theseprobabilities are shown in Table III. In a case that any groupelement τ works as a multiplier on any bond, I would write p b = 0 in Table III, rather than use (4.4). In such a case, ourmodel is locally trivial: exactly n G N configurations may beaccessed, simply by applying one arbitrary group element atevery site. In other words, there is a (locally) 1-to-1 mappingto the trivial model in which every site has an independentdegree of freedom. That model is just the gauge model, whichwas studied previously in Ref. [23].In the spirit of the Pauling approximation, let us now pre-tend that arrows on different bonds are uncorrelated : withinthat assumption, we must obtain the same cluster distributionas in the (thoroughly studied) problem of uncorrelated perco-lation on these lattices. It follows that the updating behaviortends to depend on the relation of p b to the critical percolationfraction p cb . On the one hand, if p b > p cb , then the clustergrows without limit, including a nonzero fraction of the wholesystem; in that case, the update move certainly is not efficient.On the other hand, if p b /p cb is too small, we never get a clus-ter at all, or else a single-site update (next subsection) wouldsuffice. The “interesting” case when a cluster update is neces-sary and helpful, would be for p b /p cb close to or slightly lessthan unity. C. Single-site updates A single-site update is the case that the updated cluster D is just one site, thus only the z spins around it are updated.When p b is much less than p cb , most clusters are small, andthe probability P of a single-site update is appreciable. If P is large enough, it might be ergodic to use only single-siteupdates (i.e. to pick s max = 1 ), in which case we can omitthe cluster-growing algorithm. I shall concentrate on thesecases, which are the easiest to simulate (and also the likeliestto extend to quantum models).To estimate P , I pretend the z bonds around a site are inde-pendently occupied by randomly chosen elements of S . Then P = 1 − Y α (1 − q zα ) n α (4.5) TABLE III: Entropy and updatability parameter estimates for se-lected models. Formulas from eqs. (3.1), (4.4), and (4.5). Note a : inthese cases, any even element τ can always update, but no odd τ canever update.Model name f e exp( S Pauling ) p b /p cb P Z × Z (2)tri 2/9 / ... 0.241 Z × Z (2)sq 1/3 / S (2)sq 1/3 a S (2)hc 1/3 a S (2,3)tri 4/25 / S (2,3)sq 21/125 / Q (2)sq 7/54 / D { m, m ′ } sq 1/4 a D { m, m ′ } hc 1/2 √ a A (3)tri 1/16 A (3)sq 3/32 A (2)tri 2/225 / A (3)tri 7/400 / A (5)tri 5/144 / A (2)sq 71/3375 / A (3)sq 147/8000 / A (5)sq 53/1728 / where α indexes the group class (with n α elements) that τ might be in (excluding the identity), and I defined q α to be thefraction of times τ ∗ σ ∈ S given that τ falls in class α .To digest the implications of (4.5), let’s make an evencruder version of the estimate, replacing q α by it’s averageover all τ ′ s , namely q α → − p b : I get − [1 − (1 − p b ) z ] n G − which is a lower bound on P as given by (4.5). Evidently,to have a high single-site success rate, we want (i) n G as largeas possible, (ii) z as small as possible, and (iii) p b as smallas possible; in light of (4.4), the third criterion amounts towanting n S /n G as large as possible. Those are the same threeconsiderations given in Sec. III C as favoring a large entropy.I include these estimates in Table III, particularly focus-ing on the models using group A . We see from Table IIIthat P is large enough in many cases that we can rely onsingle-site updates. However, whenever P gets close to1, our model is “too easy” in some sense – it is practically agauge model, with only mild constraints eliminating some ofthe configurations.The entry P = 1 for A (3) is delusory. This comes be-cause q α = 1 for a certain class of update multipliers, namelythe order-2 class (double pairwise exchanges). If we limitedourselves to this class, indeed every update would be success-ful, but (it can be checked) the move would not be ergodic(does not access the whole ensemble). A similar situation ap-plies in the cases of S (2) or D { m, m ′ } : any τ from theeven subgroup is always accepted, while an odd τ is never ac-cepted; but single-site updates based on the even subgroup donot access the whole ensemble.To implement an actual simulation, one would not wantto choose τ at random, but biased towards the group classes0with a larger q α (the success fraction looking at just an iso-lated bond). In particular, one would omit group classes with q α = 0 ; if the group contains even/odd elements and S in-cludes only one parity of element, then q α = 0 for every classof odd elements. The values of p b and P in Table III for S (2) and D ( m, m ′ ) were computed assuming τ ∈ E .Another way to implement a single-site update is, afterchoosing a random vertex r , to examine the local environmentof its z bonds, find the entire list of τ ’s which can update it,and choose randomly from this list. Typically, configurationdependent choices like this are avoided in Monte Carlo algo-rithms because they tend to violate detailed balance. In thepresent case, however, it can be checked that the number ofpossible τ ′ s is always the same in the old and new configu-ration, i.e. the rate is the same for the forward and backwardstep, which is sufficient to ensure detailed balance (and anequal ensemble weight for every configuration). D. Criteria for initial conditions
In height models, certain special states (e.g. the “columnar”arrangement of dimers on the square lattice) were “ideal” inhaving the maximum number of possible update moves. (Fora model requiring loop updates, we might replace that cri-terion by “having the shortest typical loops.”) Certain otherstates (e.g. the “herringbone” packing of dimers) were inert,in that no finite updates are possible (in the thermodynamiclimit). These states, in a height model, correspond respec-tively to a zero coarse-grained gradient of the height variable,or the maximum gradient.In a non-abelian height model, the coarse-grained heightgradient is undefined, but one can still construct “ideal” and“anti ideal” states. It is recommended that simulation runs bestarted in both kinds of state, being in some sense oppositeextremes of the configuration space. A diagnostic for equi-libration is then whether the expectations from the two startsconverge to the same values.More exactly, rather than a single domain of anti-ideal state,one should divide the system into two domains. Then, updatesare initially possible along the domains’ border, using loopswhich extend across the system. Gradually, a larger fraction ofthe system’s area become updatable, and the loops get smaller.On the other hand, starting from an ideal state, the loops areinitially small and get larger. Thus, tracking the loop distri-bution is an obvious diagnostic to test for convergence to thesame equilibrium state.
V. POSSIBLE MEASUREMENTS IN SIMULATIONS
In this section, I sketch how one might confirm the topolog-ical order numerically, or measure other interesting quantities,given a working Monte Carlo simulation.
A. Correlation functions
Correlation functions are an obvious starting point. Ofcourse, a topological order state has exponentially decayingcorrelations, so this serves primarily as a negative test: wecheck that the system is not a height model in disguise (seeSec. III B), which would have power-law correlations, and thatit doesn’t have long-range order (which can emerge even inequal-weighted entropic ensembles, or because the definingconstraints are too restrictive). Correlations are also of inter-est near a critical point where long-range or quasi-long-rangeorder emerges.In models with vector spins s i , one was accustomed to eval-uating the expectation of s i · s j , or occasionally its secondmoment. It may not be immediately obvious what to mea-sure now. One can, of course, simply tabulate frequencies ofdifferent combinations, e.g. (for the “height difference”) howoften γ → r belongs to each conjugacy class. It is preferable,though, to reduce the measurements to a single (meaningful)number, and the appropriate generalization of the dot productis the trace of the matrices in the right group representation.Thus we are led to use a character function χ ( x ) , where x is any group element; this is always the same within eachconjugacy class of the group. I divide the actual character bythe dimension of the representation, so that χ ( e ) ≡ for anyrepresentation, and | χ ( x ) | ≤ for any element. Presumably,the best choice of representation is the one that has the largestpositive χ ( σ ) for spins (for σ ∈ S ). This corresponds concep-tually to using a distance metric, within the group G , countingmany multiplications by some element of S are needed to takeyou from element to the other one.
1. Height difference correlation
In the old “height models” (sketched in Sec I A), a naturalmeasure of fluctuations was h| h )0) − h ( r ) | i . The natural gen-eralization of this for the present models with finite (possiblynon-abelian) groups is C h ( r ) ≡ h χ ( γ ( ℓ → r )) i . (5.1)Of course, the product γ ( ℓ → r ) is independent of which pathis taken from to r – provided the path does not wrap aroundthe periodic boundary conditions.As just noted, choosing χ ( . ) so that χ ( σ ) is as close toone as possible, provides that C h ( r ) does express how fastthe group element wanders from the identity under repeatedcompositions; that is the choice likeliest to give a monotonicdecay with distance. If γ ( ℓ → r )) is equally likely to be anygroup element – which one expects large r – then it followsthat C ( r ) = 0 .
2. Spin correlations
Similarly, we can compute G ij ≡ h χ ( σ i ∗ σ − j ) i . (5.2)1 B. Defects
It is easy to augment the simulation to allow a defect pla-quette where the plaquette constraint is violated. The same(single-site) update rules will work correctly next to the de-fect, but they cannot change its position. To make a defectmobile, one can add additional update rules specific to thedefect, by (say) arbitrarily choosing one bond of the plaque-tte and changing it to make the plaquette’s loop product be e (which, of course, the loop product not be e for the pla-quette on the other side of that bond, unless that was also adefect plaquette and this is the annihilation event.) The sim-ulation would normally be run with a constraint or bound onthe number of defects.The idea is to create a pair of defects, by hand, and thenevaluate expectations depending on them. The first thing tomeasure is the distribution P d ( R ) of defect separations R . Inthe case of topological order, we expect deconfinement, mean-ing P d ( R ) → const for R > ξ , where ξ is a (not very large)correlation length. One can define an effective (entropic) po-tential V ( R ) by P d ( R ) ∝ exp( V ( R )); (5.3)physically, V ( R ) is the difference in entropy due to plac-ing the defects near to each other. In the case of a heightmodel, V ( R ) ∝ ln | R | , and d P ( R ) decays to zero as a powerlaw. [44]In fact, since there are various flavors of defect labeled bydifferent group elements b , one really needs to write the effec-tive potential as E = U ( b ) + U ( b ′ ) + V b,b ′ ,c ( R ) (5.4)where b and b ′ are the respective defect charges, and c is thenet charge of the combined defect. Here, U ( b ) and U ( b ′ ) are“core energies” of these respective defects; these, and usually the inter-defect potential, are functions only of the conjugacyclasses of b , b ′ , and/or c . Implicit in the form (5.4) is that theexponential confinement length probably depends on all of b , b ′ , and c .Measuring how the effective potential depends on class ismore physical, since (i) it decides whether a defect is stableagainst decays into other defects (ii) measurements of defectbehavior (in simulations or in real systems, were any to bediscovered) might be used to discover the universality class ofthe topological order, if that were not known. I conjecture thatthe dependence on b , b ′ , and c is also described by a characterfunction; it would be interesting to see if that can be exploredanalytically in some model.Incidentally, since (with the appopriate boundary condi-tions) we can have a single defect in our system cell, that givesadditional opportunities to evaluate e.g. the core energy U ( b ) without the complication of a second defect.Finally, if the single-site updates of Sec. IV C are not feasi-ble, defects provide a less elaborate alternative update move,in place of the cluster update of Sec. IV A. Namely, we cre-ate a pair of defects and allow them to random-walk until theyannihilate. (However, if their paths differ by a loop around the periodic boundary conditions, they may be unable to an-nihilate.) Many Monte Carlo schemes [31, 32] are based on asimilar process. C. Topological sectors
The tests of topological order outlined up to here have beennegative; none of them catpures the positive property of topo-logical order, which is the degeneracy of topological sectors.This can be measured in a classical simulation, if we use a(necessarily nonlocal) update which can change sectors, whilesatisfying the detailed balance condition. Either the clusterupdate of Sec. IV A or the defect-pair update just outlined inSec. V B will suffice.From the relative fraction of time spent in different topo-logical sectors, we can infer a free energy F L ( γ x , γ y ) , where ( γ x , γ y ) are the loop products characterizing the sector, and L refers to the system size. This is a finite size effect, since (bydefinition of topological order) the difference between sectorsvanishes in the thermodynamic limit; F L is expected to decayexponentially as a function of L .In a similar fashion, if we allow transitions between stateswith and without a defect as part of the dynamics, we canevaluate the defect core energy U ( b ) . Of course, F L ( γ x , γ y ) is very analogous to U ( b ) , since b is a loop product encirclingthe puncture where a defect sits, just as γ x is from the loopproduct encircling the system. [45] VI. TRANSFER MATRIX AND ANALYTIC APPROACHES
In a quantum mechanical models with topological order,the energy differences between different topological sectorsdecays with system size as exp (const L ), and the correlationsof a defect pair decay as exp( − R/ξ ) . Up to now, I have as-sumed without justification that this would carry over to thepresent classical models.This section finally examines the basis of exponential be-havior. I turn here to an analytic treatment using the (practi-cally) one-dimensional framework of transfer matrices. Firstof all, this sheds some light on why the finite-size depen-dences, as well as the defect-defect interaction, are exponen-tially decaying with distance. More specifically, they clar-ify the pattern of how sector-weight splittings or defect-pairdistributions relate to the group’s representations and symme-tries. A. One-dimensional model
Imagine the most trivial system which can have topologi-cal sectors: the one-dimensional version of the discrete-groupheight models. There can be no plaquette constraint. Our en-semble simply consists of chains of length L – with periodicboundary conditions – having a group element σ i placed oneach link, the only constraint being that σ ∈ S . All ( n S ) L sequences are equally likely.2If we let γ ( x ) ≡ σ x ∗ σ x − ∗ ... ∗ σ . (6.1)then the topological sectors are labeled by γ ( L ) . Define the( n G × n G dimensional) transfer matrix T in the standard fash-ion: let T γ ′ ,γ be the number of ways to get γ ( x + 1) = γ ′ from γ ( x ) = γ . Then ( T L ) γ,e is the partition function (thetotal number of states) for the sector with γ ( L ) = γ . Notethat T commutes with permutations that implement the sym-metry operations (automorphisms) of the group G ; hence, theeigenvalues/eigenvectors of T are classified by the represen-tations of the automorphism group (mentioned in Sec. III A).The transfer matrix has eigenvalues { Λ k } and correspond-ing eigenvectors { v k,m } ; the index m labels each family ofsymmetry-related eigenvectors belonging to the same (degen-erate) eigenvalue Λ k .The restricted partition function for topological sector γ is X k,m [ v k,m ] γ [ v k,m ] e Λ Lk . (6.2)Hence, in any sector the overall (entropic) free energy perunit length is ln Λ , where Λ is the largest eigenvalue,and L -dependent corrections depend on some larger eigen-value Λ s . For this trivial one-dimensional model, v =(1 , , ..., , / √ n G . More generally v must be totally sym-metric under all automorphisms of G , i.e. it belongs to thetrivial representation. Indeed, v s must also belong to the triv-ial representation, since [ v k,m ] e in (6.2) is independent of m ,but P m [ v k,m ] γ = 0 for any other representation. We let Λ s be next largest (necessarily nondegenerate) eigenvalue of thefully symmetric representation, after Λ .Hence, P ( γ ) P ( e ) ≈ c γ (Λ s / Λ ) L c e (Λ / Λ ) L . (6.3)where c g ≡ [ v s ] g [ v s ] e [ v ] g [ v ] e (6.4)where [ v ] g [ v ] e = 1 /n G , for this one-dimensional model. Itfollows from (6.3) that ln P ( γ ) P ( e ) ≈ ( c γ − c e ) e − L/ξ (6.5)where exp( − /ξ ) ≡ | Λ s / Λ | . Often Λ < ; in this case,we must add a factor ( − L on the right-hand-side of (6.5).Furthermore, at short L , we may see subdominant terms withshorter decay lengths ξ etc., deriving from other eigenvectorsof T .
1. Example
A useful example is any group G when S = G \ e ,i.e. every element but the identity is allowed. In this case T = (1 , , ..., ⊗ (1 , , ..., − I . Thus Λ = n G − and Λ = Λ = ... = − . Thus exp( − /ξ ) = 1 / ( n G − andthe deviations have alternating signs, i.e. the ( − L factor isneeded in (6.5). For the model Z × Z (2) , the matrix is T = . (6.6)
2. Symmetry-class reduced matrix
We can classify eigenvectors as “symmetric” or “asymmet-ric” according to what representation of the automorphismgroup they transform under. “Symmetric” eigenvectors areinvariant under group symmetries, while “asymmetric” eigen-vectors represent a bias of the probability distribution favoringcertain local patterns over other (symmetry-related) ones.Since the sector probability ratio (6.3) is the same for allsymmetry-related γ , I believe not only v − but also v mustbe totally symmetric. That affords a considerable simplifica-tion, for we can replace T by its projection e T onto the groupelement symmetry classes. (Such a class consists of elementsthat map to each under under some automorphism, so theseare at least as large as the conjugacy classes.) Whereas thedimension of T was the number of group elements n G , thedimension of e T is the number of group classes: e T ji tells thenumber of times that σγ belongs to class j , if γ belongs toclass i and σ runs over all n S elements in S .For example, in the case of the group A , the matrix is re-duced from × to × , with entries for elements oforder one (identity), two, three, and five. (There are two con-jugacy classes with order five, but they are equivalent by anouter automorphism.) For the model A (3) , we get e T = . (6.7)This matrix is similar to a symmetric matrix D − / e T D / ;here D = diag(1 , , , for this group, or in general isthe diagonal matrix with entries being the count of each class. B. Sector probabilities in two dimensions?
A two-dimensional finite-group height model is also de-scribed by a transfer matrix T . However, now the vector that T acts on represents all possible path products γ x,y taken to apoint ( x, y ) , and thus is ( n G ) W dimensional, where W is thewidth of the strip (in the y direction; iteration still runs in the x direction). We must replace c γ → c γ ( W ) and ξ → ξ ( W ) inEq. (6.5). Conceivably ξ ( W ) → as W → ∞ , as is very wellknown in gapless systems, so the form of exp( − L/ξ ( W )) does not prove exponential decay in d = 2 .3Nevertheless, we can make a plausible guess to obtain afitting form for comparison with numerics. Since all correla-tions are expected to be rapidly decaying, a strip of width W is like W/w independent, one-dimensional strips of width w in parallel. But all these strips are constrained to have thesame, or equivalent, sector label γ .) The consequence is that P ( γ ) P ( e ) ≈ " c γ (Λ / Λ ) L c e (Λ / Λ ) L W/w (6.8)so ln P ( γ ) P ( e ) ≈ ( C γ − C e ) W e − L/ξ (6.9)in place of (6.5), with C γ ≈ c γ /w and ξ ≈ ξ independentof W .My chief motivation for introducing the transfer-matrix for-malism is separate from such guesses about the W scaling,and is much better founded. Namely, the eigenvectors forthe corrections to P ( γ ) are representations of the automor-phism group. Furthermore, which representation goes withthe longest correlations is probably the same as in the one-dimensional case. What really matters here is that our choiceof a selected set S defines a sort of metric on G : the distancefrom g to g ′ is the number of times you need to multiply byan element of S to get from g to g ′ . Then, the first nontriv-ial eigenvector v is the mode that is slowest varying on G according to this metric (apart from v which is uniform).The one-dimensional correlation length ξ can be computedfor any combination of G and S and can serve as another “fig-ure of merit” for a group. That is, in light of the previous para-graph’s argument, it should be roughly related to the true sec-tor probability decay length ξ for the two-dimensional model(and likely related to the defect-defect decay length as well). C. Defect separations and W = 2 transfer matrix Whereas the one-dimensional model already seems to cap-ture the essence of how sector probabilities depend on systemsize and sector label, it does not admit topological defects andhence sheds no light on the parallel question of how p ( R ) fora defect pair decays with separation or depends on the respec-tive defect charges.Clearly, p ( R ) must be associated somehow with the eigen-vectors and eigenvalues of the two-dimensional transfer ma-trix, since all possible correlation information is expressed init. But it is not self-evident just what kind of distortion of theensemble is being propagated, or what sort of subdominanteigenvector: the symmetric kind (which governed the sectorprobabilities) or the asymmetric kind.I will work out here a toy calculation, again using a transfermatrix, of the correlation decay due to asymmetric eigenvec-tors. I believe they are the ones that matter for the case ofan abelian group. In that case, the charge of a defect is aparticular element: the loop product around the defect givesthat same result, no matter how big the loop, and only another defect with the inverse of that charge can cancel it. In thenon-abelian case, however, these properties would seem to bedefined only modulo conjugacy classes. Therefore, the picturepresented here is only asserted to go with abelian groups.The simplest property that could influence or be influencedby a defect’s presence is the correlation of two adjacent spinson the same plaquette, i.e. sitting on bonds that make a 90 ◦ angle. A simple example is the model Z × Z (2) sq, in which n S = 3 elements are allowed – all except the identity. Theseelements are { a, b, c } . Consider a plaquette with the spinson two edges specified and the remaining two spins to be as-signed (there are unconstrained ways to do so). When ad-jacent edges on a plaquette have the same element, there arethree ways to satisfy the plaquette constraint, but only twoways if the given adjacent edges are different. On the otherhand, if we want to make a defect plaquette, there are six wayswhen the given spins are the same but seven ways when theyare different.Let’s set up a W = 2 strip, the narrowest kind that cancapture defect correlations. This transfer matrix, unlike theprevious one, refers to the actual spin configurations in eachvertical pair of bonds; we add up all the possible horizontalbonds. I assume the upper row of plaquettes are constrainedto be have identity product around the plaquette. Plaquettes inthe lower row can have any product – defects are permitted –with a weight θ for the identity or θ a , θ b , θ c for the respectivedefect charges a, b, c . We imagine the limit in which θ a,b,c aresmall and ask for the corresponding defect correlations.Although T has × = 81 elements, in fact there areonly ten distinct kinds by symmetry, as given in Table IV;“no.” represents the number of times each kind occurs in thematrix. The factors θ ≈ and θ σ ≪ are omitted in the ta-ble. To compute the matrix elements, note that when σ = σ ′ in the upper plaquette, the central horizontal bond may be anyelement [three possibilities] but if σ = σ ′ , the central bondmay be only σ or σ ′ [two possibilities]. There are alwaysthree possibilities for the lower horizontal bond, so the table’srows add up to 9 or 6 depending whether or not σ = σ ′ .The probability to find a defect of charge β ′ at separation R , given there is a defect of charge β at the origin, is then p ( R ) = Tr (cid:16) [ T (0) ] L − R − T ( β ′ ) [ T (0) ] R − T ( β ) (cid:17) Tr (cid:16) [ T (0) ] L − T ( β ) (cid:17) (6.10)For a large power M , we can replace [ T (0) ] M → (Λ ) M v ⊗ v + (Λ ) M v ⊗ v . Here v and v are the eigenvectors be-longing to the maximum and next-largest eigenvalues of T (0) .Assuming L ≫ R ≫ , we get p ( R ) = Λ L − R − P k =0 , P m ( v , T β ′ v k,m )( v k,m T β v )Λ R − k,m Λ L − ( v , T β v ) (6.11) = p ( β ′ ) h c ( β ′ ) c ( β ) (cid:16) Λ Λ (cid:17) R i (6.12)where p ( β ′ ) = ( v , T ( β ′ ) v )Λ (6.13)4 TABLE IV: Example for group Z × Z (2): Transfer matrix elements T ( β ) σ ′ ,σ ′ ; σ ,σ no. ( σ , σ ) ( σ ′ , σ ′ ) T (0) T ( a ) T ( b ) T ( c ) ( aa ) ( aa ) ( aa ) ( ab ) ( aa ) ( ba ) ( aa ) ( bb ) ( aa ) ( bc ) ( ab ) ( ab ) ( ab ) ( ba ) ( ab ) ( ac ) ( ba ) ( ca ) ( ab ) ( ca ) — remember ( v , T (0) v ) = Λ ] — and c ( β ) ≡ (cid:16) Λ Λ (cid:17) / ( v , T ( β ) v )( v , T ( β ) v ) . (6.14)Please remember, the eigenvector called v here is asym-metric, and is thus not the same as the symmetric eigenvectorcalled v s in Sec. VI A. We see that asymptotically, ln p ( R ) ∝ c ( β ) c ( β ′ ) e − R/ξ (6.15)Notice first that the decay length ξ is independent of the defectcharges, but different defect charges have different projectionsonto this eigenmode. As is clear from the derivation, a moregeneral form could be written, including subdominant contri-butions: ln p ( R ) ∝ X k c k ( β ) c k ( β ′ ) e − R/ξ k (6.16)where ξ > ξ > ... . The later terms could be importantcorrections to include in fits at short R , particularly when thesmaller ξ k ’s happen to be associated with larger coefficients c k ( β ) . Also, if c ( β ) = 0 for certain defects, their asymptoticinteraction gets carried by the first mode that has nonzero pro-jections onto both defects.The formulas basically apply to any width of strip. (If de-fects are allowed in more than one horizontal row of plaque-ttes, then the defect distribution is no longer a function justof R but also of the two y coordinates; the only modificationnecessary is that T ( β ) → T ( β ; y ) , labeled not only by the de-fect’s flavor but by its y coordinate.) I would speculate thatthe W = 3 strip, with defects in the central row, may be agood approximation in practice, although of course there isno control parameter to make small. The basis for this is sim-ply the notion that, when we have rapid exponential decays,these are associated in the ensemble with strings connectingthe defects; any influence carried by a less direct chain wouldbe exponentially smaller in correspondence with the longerlength. D. Other approaches to p ( R ) in d = 2 I conjecture there is an alternative approach which is morecongenial to d = 2 . Namely, in the vicinity of a defect, theprobabilities of local patterns have small deviations from thebulk values, which could be represented by operators O k ( r ) and small conjugate fields h k ( r ) . That is, adding a Hamil-tonian P r h k ( r ) O k ( r ) (in the absence of the defect) wouldperturb the ensemble the same way as the defect does. Notethat the operator “ O k ( r ) ” is schematic, in the sense that suchoperators probably involve two spins at different r (in light ofthe same logic laid out in the first paragraphs of this subsec-tion). Then possibly some sort of mean-field approximationproduces a difference equation for h k ( r ) , the discrete ana-log of Poisson’s equation ∇ h k ( r ) = h k ( r ) /ξ k , which hassolutions ∼ e − R/ξ k /R . In this approach, we have a sort ofsmall parameter in that the influence of a perturbation decaysas e − R/ξ k , which becomes arbitrary small at sufficiently large R . We can therefore rely on linear response in that regime.One can conceive additional approaches to p ( R ) whichdepend on a genuine small parameter; the difficulty is thatthe actual model families defined in this paper are far fromthat limit. For example, one could expand around the puregauge theory: in place of the spin constraint, there wouldbe no constraint but configurations would have a statisticalweight exp( λ P r , r ′ u ( σ ( r , r ′ )) , where u ( σ ) would penalizeall σ / ∈ S . In the limit λ → , we have a pure gauge the-ory in which all correlation lengths ξ k are zero, so hopefully ξ k would scale as a power of λ . The models under consid-eration are, unfortunately, the case λ = ∞ . Still, since thetopological phases are like the pure gauge models at largescales, they should be adiabatically connected and hence thisapproach should be qualitatively valid. VII. DISCUSSION
I have put forward the notion of purely classical topologicalorder, defined by an ergodicity breaking into sectors depen-dent on the topology, and not distinguishable by thermody-namic expectations of any local operator. A family of explicitmodels has been described, along with a suitable Monte Carlotechnique, and criteria were suggested to pick out the mostpromising cases (in having a nontrivial and updatable ensem-ble of allowed states).A framework was set up (Sec. II A) to define models withthree variable attributes: which group, which class(es) out ofthe group to selected as “spin” variables, and which latticeto place the model on. These are characterized by parame-ters – the sizes n G and n S of the group and the spin subset,the coordination number z of the lattice and z d of its dual– which entered crude formulas that estimate the entropy ofthe model and its updatability under single-site Monte Carlomoves (Sec. III C and IV B), which are the only actual cal-culations in the paper. Groups with normal subgroups tend tobe “less nonabelian ”, thus perhaps less attractive (Sec. III B).Although topological order superficially would appear to beintrinsically featureless, there is sufficient richness of measur-5able functions when one considers the dependence of free en-ergy on topological indices – finite size dependence on sectoror finite distance dependence on defect separations (Sec. V).In trying to connect the classical picture to the quantumtheory of topological order, it is intriguing that a given twodefect charges (see Sec. II C) can combine in more than oneway, in a classical non-abelian model, reminiscent of fusionrules in a quantum model. If further investigation finds that thesector counting gives the same degeneracies in the classical asin the quantum case, one would conclude that this is one ofthe shared properties, not an intrinsically quantum one.The physical manifestations of classical topological orderand/or of non-abelianness are less striking, perhaps, than forthe quantum case. Most prominent is the behavior of topolog-ical defects. Topological order implies deconfinement in theclassical model for nonabelian and abelian cases alike. Non-abelianness (of the group) changes the rules for addition ofdefect charges, and braiding has physical consequences, eventhough there are no Berry phases in a classical model. (Itmust be noted, however, that the same behaviors are seen innon-abelian defects of ordinary long-range order [28] – theyare not inherent to topological order.)Degeneracies of different topological sectors, the definingproperty of topological order, work differently in the non-abelian than in the abelian case: for example, there are fardistinct fewer sectors in the non-abelian case (Sec. II B). A. Quantum mechanics
Several central concepts of topological order are inherentlyquantum-mechanical and thus have no counterpart in classi-cal topological order. They are mainly related to phases inwavefunctions and braiding of worldlines in 2+1 dimensions,namely anyon and mutual statistics. Most real or imagined ex-periments relating to topological excitations have involved in-terference phenomena (e.g. tunneling in various geometries ofquantum Hall fluids) and thus probe the quantum-mechanicalaspects of topological order.Another feature missing in the classical models is the dualdefect or quasiparticle (such as the vison [39]), which is adistortion of the phase factors in the many-body wavefunction.A final attribute of topological orders is the nontrivialcounting statistics of the excited states made by several quasi-particles, which is quantum mechanical in that it concerns thelinear dimension of a Hilbert space. One cannot rule out theappearance of similar concepts in classical stat mech – there,too, the partition function contains combinatorial factors forthe placement of defects, after the other degrees of freedomehave been integrated out. However, I am not aware of a clas-sical situation in which such a nontrivial counting actuallyemerges. B. Construction of quantum models?
Any of the classical height models with topological ordermay be converted into a simular quantum model if we can endow it “flipping” move, just as classical dimer (and other)models get converted into quantum dimer models using theRokhsar-Kivelson (RK) prescription [37]. A barrier to thisis that the only generally guaranteed “flip” move is a clusterupdate, as explained in Sec. IV A.Fortunately, whenever the single-site update (Sec. IV C)suffices, we can define a quantum model with a simple “flip-ping” term in the Hamiltonian, usually parametrized by anamplitude t , as well as a “potential” term of strength V = t that penalizes each flippable place. At the RK point V = t ,the ground state wavefunction is a superposition of all config-urations in the same topological sector, with the same (equal)weighting as in the classical ensemble, and (mutually inac-cessible) topological sectors are trivially degenerate. One isalso free to set V = 0 – obtaining a simpler model in whichflippable sites are so favored that an ordered state is likely tobe the outcome – or to vary V /t with the hope of crossing aphase transition.The above recipe is incomplete, in that there are many pos-sible choices of update (labeled by the multiplier τ of Sec. IV),and presumably all or many should be included in “flipping”term of the quantum Hamiltonian, which requires a prescrip-tion for the relative magnitudes of coefficient to put for eachclass of τ , as well as the relative phase factors. Presumably, aproper choice is taking the same phase factors for every term,i.e. the Hamiltonian transforms by the fully symmetric (triv-ial) representation of the automorphism group of G . Alterna-tively, in lucky cases, one might select a site-dependent pat-tern of τ i ’s so as to link the group symmetry to the lattice sym-metry, in the spirit of Kitaev’s honeycomb model [4]. Anotheroption is to include a second, quantum fluctuating field of τ ′ s which are used for the update. If the τ ′ s are derived from asecond set of “spins” also having the gauge-like structure ofa finite-group height model, we might be able to realize dual(“magnetic”) defects having mutual statistics with the σ -spintype (“electric”) defects described in this paper. C. Possible simulations
As laid out in Sec. V, several quantities e.g. correlationfunctions can be measured in classical non-abelian heightmodels as a test-bed, whereas the analogous calculation mightbe very challenging computationally in a quantum mechanicalmodel. Of course, the answers need not be the same, but thequestions may be much clearer once the classical results arein hand. First, one can create defect pairs and evaluate the his-togram of their separations, which will reveal whether or notthey are deconfined. Secondly, one can evaluate the probabil-ities of different topological sectors, which is the direct test oftopological order.Furthermore, if we generalize the models to include classi-cal Hamiltonians (so as to weight configurations according tothe Boltzmann distribution), phase transitions can be studied.Just as a standard height model may have a “smooth” phase,in which one or more height components becomes locked, itseems conceivable that a discrete-group height model mighthave a phase in which the loop products can take values in a6subgroup G ′ ⊂ G . If so, one might encounter critical pointsseparating different topological phases, and characterize thecritical exponents. D. Dilution and effective interactions of local degrees offreedom?
A classical model might be a helpful too for investigatingthe consequences of dilution disorder in a model with topolog-ical order. Each diluted site or plaquette is like a very smallhole cut in the system, thereby increasing the genus and thenumber of topological sectors. If the hole were big, these sec-tors would be truly degenerate (by the definition of topologicalorder), however this degeneracy is broken since the holes aresmall.The values of γ ∗ on each dilution site are local pseudospins,which are expected to have (exponentially decaying) interac-tions mediated by the fluid in between them, like the emergentspin-1/2 degrees of freedom in diluted spin-1 antiferromag-netic chains [35]. The ground state of such a system could beconstructed by a renormalization group that iteratively com- bines the most strongly coupled pair of pseudospins into a sin-gle effective pseudospin, as was originally done for the (expo-nentially decaying) antiferromagnetic coupling of the charge-bound electron spins in P-doped Si [36].In the non-abelian case, at least, we do not know for surewhether these interactions lead to an inert singlet phase (as inthe antiferromagnet) or could give a state with some kind oforder among the pseudospins. Thus the system with defectswould not be a topological liquid, and this would be a novelscenario of how order can emerge due to disorder. [38] Acknowledgments
I thank M. Troyer for suggesting the problem; also L. Ioffe,A. Kitaev, D. A. Ivanov, Simon Trebst, C. Castelnovo, C. Cha-mon, S. Papanikolaou, and R. Lamberty for comments anddiscussion. ALso, I thank J. Papaioannou and R. Maimon forpreliminary work on the simulation algorithm. This work wassupported by NSF Grant No. DMR-1005466. [1] X.-G. Wen, Phys. Rev. B44, 2664-72 (1991).[2] X.-G. Wen,
Quantum Field Theory of Many-body Systems (Ox-ford University Press, Oxford, 2004), chapters 8 – 10.[3] L. B. Ioffe, M. V. Feigel’man, A. Ioselevich, D. Ivanov, M.Troyer, and G. Blatter, Nature 415, 503-506 (2002).[4] A. Kitaev, Ann. Phys. 303, 2-30 (2003).[5] G. Misguich and C. Lhuillier, “Two-dimensional quantumantiferromagnets”, in
Frustrated spin systems , ed. H. T.Diep (World Scientific, Singapore, 2005); see ArXiV/cond-mat/0310405.[6] M. Freedman, A. Kitaev, M. J. Larsen, and Z. Wang, Bull.Amer. Math. Soc. 40, 31 (2003),[7] M. A. Levin and X.-G. Wen, Phys. Rev. B71, 045110 (2005),[8] M. Freedman, C. Nayak, and K. Shtengel, Phys. Rev. Lett. 94,066401 (2005),[9] C. Castelnovo and C. Chamon, Phys. Rev. B 76, 174416 (2007).[10] S. Sachdev,
Quantum Phase Transitions (Cambridge UniversityPress, 1999).[11] R. Z. Lamberty, S. Papanikolaou, and C. L. Henley, unpub-lished.[12] H. van Beijeren, Phys. Rev. Lett. 38, 993 (1977),[13] H. W. J. Bl¨ote and H. J. Hilhorst, J. Phys. A 15, L631 (1982)[14] B. Nienhuis, H. J. Hilhorst, and H. W. J. Bl¨ote J. Phys. A 17,3559 (1994).[15] W. Zheng and S. Sachdev, Phys. Rev. B 40, 2704 (1989): L. S.Levitov, Phys. Rev. Lett. 64, 92-5 (1990).[16] Jan´e Kondev and C. L. Henley, Phys. Rev. B 52, 6628 (1995).[17] J. Kondev and C. L. Henley, Nuc. Phys. B 464, 540 (1996).[18] C. Zeng and C. L. Henley Phys. Rev. B 55, 14935 (1997).[19] R. Raghavan, C. L. Henley, and S. L. Arouh, J. Stat. Phys. 86,517-550 (Feb. 1997).[20] J. K. Burton and C. L. Henley, J. Phys. A 30, 8385-8413 (1997).[21] B. Nienhuis, in
Phase Transitions and Critical Phenomena ,edited by C. Domb and J.L. Lebowitz (Academic, London,1987), Vol. 11. [22] D. R. Nelson, in