Topological multicritical point in the Toric Code and 3D gauge Higgs Models
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] A p r Topological multicritical point in the Toric Code and 3D gauge Higgs Models
I.S. Tupitsyn, A. Kitaev, N.V. Prokof’ev, and P.C.E. Stamp Pacific Institute of Theoretical Physics, University of British Columbia,6224 Agricultural Road, Vancouver, BC V6T 1Z1, Canada California Institute of Technology, Pasadena, California 91125, USA Department of Physics, University of Massachusetts, Amherst, Massachusetts 01003, USA (Dated: April 20, 2008)We report a new type of multicritical point that arises from competition between the Higgs andconfinement transitions in a Z gauge system. The phase diagram of the 3d gauge Higgs model hasbeen obtained by Monte-Carlo simulation on large (up to 60 ) lattices. We find the transition linescontinue as 2nd-order until merging into a 1st-order line. These findings pose the question of aneffective field theory for a multicritical point involving noncommuting order parameters. A similarphase diagram is predicted for the 2-dimensional quantum toric code model with two external fields, h z and h x ; this problem can be mapped onto an anisotropic 3D gauge Higgs model. PACS numbers:
Introduction.
Topological quantum phases and anyonsare well known in connection with the fractional quantumHall effect, but they are also expected to exist in frus-trated magnets. It has long been proposed that a certainclass of resonating-valence-bond (RVB) [1] phases carries Z -charges and vortices and has a four-fold degenerateground state on a torus [2]. A qualitative understand-ing of this phase can be obtained from a so-called toriccode model (TCM) [3]. The dimer model on the Kagomelattice is mapped onto the TCM exactly [4] while someother models [5, 6] belong to the same universality class.The TCM is defined in terms of spin-1 / d lattice.The Hamiltonian is as follows: H T C = − J x X s A s − J z X p B p , (1)where A s = Q j ∈ s σ xj and B p = Q j ∈ p σ zj are products ofspin operators ( σ αj are the Pauli matrixes) on the bondsincident to a lattice site s and on the boundary of a pla-quette p , respectively. The ground state corresponds toeigenvalues A s = 1, B p = 1 for all s and p . On a surfaceof genus g , it is 4 g -fold degenerate. Elementary excita-tions are characterized by eigenvalues A s = − Z charge on site s ) and B p = − Z vortex on plaquette p ); all excitations are gapped. Each type of quasiparticleis bosonic, but due to nontrivial mutual braiding, theymust be jointly regarded as Abelian anyons .The Hamiltonian (1) has special properties related toits exact solvability: the two-point correlator vanishesand the quasiparticles have flat dispersion. These fea-tures do not survive a small generic perturbation, whilethe topological character of the ground state and the any-onic quasiparticle statistics are robust. Yet a sufficientlystrong field can polarize the spins, driving a transitionto the topologically trivial phase. Trebst at al [7] stud-ied a perturbation of the form − h P b σ zb and solved theproblem by reducing it to the 2 d transverse-field Ising model, which is mapped to an anisotropic 3 d classicalIsing model. In this paper we consider a more generalHamiltonian: H Q = H T C − h x X b σ xb − h z X b σ zb , (2)where b runs over the bonds of a square lattice and H T C isgiven by Eq. (1). Note that the fields h x and h z inducedifferent types of phase transition. The term with h z creates virtual pairs of Z charges, which condense whenthe field strength exceeds a certain threshold. This phe-nomenon may be described as a Higgs transition, or asvortex confinement. By duality, the field h x causes thecondensation of vortices and charge confinement. Thecompetition of the two terms should result in an inter-esting phase diagram, which is the subject of this paper.We approach the problem by reducing the quantumHamiltonian to a classical anisotropic Z gauge HiggsHamiltonian on a three-dimensional cubic lattice; seeEq. (5) below. We expect the phase diagram to be qual-itatively similar to that for the isotropic case , i.e., model M , as defined by Wegner [8]. Monte-Carlo simulationshave been performed for the latter model because it ismore amenable to numerics. Some properties of thephase diagram in the isotropic case were predicted byFradkin and Shenker [9]. In particular, the topologicalphase is bounded by second-order lines corresponding tocharge condensation (for h x ≪ h z ∼ J x , J z ) and vor-tex condensation (for h z ≪ h x ∼ J x , J z ), but the twocondensate phases are continuously connected. For thequantum Hamiltonian (2), a connecting path is realizedby increasing h z so as to polarize the spins in the z di-rection, rotating the field in the xz -plane, and decreasingit again. However, the two phase transitions are clearlydifferent, therefore the corresponding lines cannot joinsmoothly.A previous numerical study involving 10 sites byJongeward, Stack, and Jayaprakash [10] showed the twolines merge into a first-order line that partially separatesthe charge and vortex condensates, and suggested the2nd-order lines might become 1st-order before merging.As discussed below, our data for larger systems (up to60 sites) do not confirm this conjecture. Thus, thepoint where all three lines meet is likely to be a noveltype of multicritical point. Note that each of the 2nd-order transitions can be characterized by an Ising-typeorder parameter, i.e., the amplitude of the correspondingcondensate. The two parameters must somehow coexistthough they are incompatible at the classical level dueto the nontrivial braiding between charges and vortices.(ie., the kinetic terms describing the charge and vortextransport do not commute.) The reduction to a classical problem.
The Hamilto-nian (2) is not gauge-invariant, but it can be mapped to a Z gauge theory by introducing a dummy spin variable µ s (“matter field”) at every site, but only considering states | Ψ i such that µ xs | Ψ i = | Ψ i . This constraint is a proto-type of the gauge-invariance condition, µ xs A s | Ψ i = | Ψ i ,which says that flipping the spin on a site and all incidentbonds does not change the state. To turn one constraintinto the other, we apply the transformation: σ zuv → µ zu σ zuv µ zv , σ xuv → σ xuv ,µ zs → µ zs , µ xs → µ xs A s , (3)where σ uv belongs to the bond connecting sites u and v .Then, the Hamiltonian becomes: H = − J x X s µ xs − J z X p B p − h x X b σ xb − h z X uv µ zu σ zuv µ zv . (4)Note that in the first term we have replaced A s by µ s using the gauge-invariance condition, µ xs A s ≡ τ = β/n , and approximatethe quantum partition function Z = Tr[exp( − βH ) P ] byTr[exp( − ∆ τ H x ) P exp( − ∆ τ H z )] n , where P is the pro-jector onto the gauge-invariant subspace and H x , H z arethe terms in the quantum Hamiltonian that depend on σ xb , µ xs and σ zb , µ zs , respectively. This expression can bewritten as a classical partition function on a cubic lat-tice. The classical variables σ b,t , µ s,t = ±
1, in each timeslice t , correspond to σ zb and µ zs respectively. But whenwe change from 2 d to 3 d , new vertical bonds (along thetime direction) appear. The classical spins on the verticalbonds between two time slices correspond to a choice ofterm in the expansion of P = Q s (cid:0) (1 + µ xs A s ) (cid:1) . Thuswe arrive at this classical Hamiltonian: H C = − X uv λ || , ⊥ bond µ u σ uv µ v − X p λ || , ⊥ pl Y j ∈ p σ j ; (5) λ || bond = −
12 ln tanh ˜ J x − vertical bonds; (6a) λ ⊥ bond = ˜ h z − horizontal bonds; (6b) λ || pl = −
12 ln tanh ˜ h x − vertical plaquettes; (6c) λ ⊥ pl = ˜ J z − horizontal plaquettes , (6d)where ˜ J x = J x ∆ τ , ˜ J z = J z ∆ τ , ˜ h x = h x ∆ τ , ˜ h z = h z ∆ τ .This model is an anisotropic generalization of the Z gauge Higgs model [9].As a final step, we eliminate the redundancy by fixing µ s . This only changes the classical partition function bya constant factor since Hamiltonian (5) can be written interms of the gauge-invariant variables S uv = µ u σ uv µ v : e H C = − X b λ || , ⊥ bond S b − X p λ || , ⊥ pl Y j ∈ p S j . (7)More detailed calculations show that the quantum andclassical partition functions are related by Z = (cid:16) sinh(2 ˜ J x ) (cid:17) k/ (cid:16) sinh(2˜ h x ) (cid:17) m/ e Z C , (8)where k and m are the number of vertical bonds andplaquettes, respectively.Of course, Eq. (8) holds only in the limit ∆ τ → J x , ˜ J z , ˜ h x , ˜ h z , even thoughthe corresponding quantum problem may not be defined.In the isotropic case, two parameters will suffice: e H C = − λ bond X b S b − λ pl X p Y j ∈ p S j , (9)where λ bond = ˜ h z , λ pl = − ln tanh ˜ h x . This model isequivalent to the isotropic Z gauge Higgs model [9] andin what follows we compute its phase diagram. Phase diagram in the isotropic case. At λ bond = 0 allconfigurations, including ground states, have the samedegeneracy factor 2 N . The actual physical variablesin this limit are plaquette numbers N p = Q j ∈ p S j , andthe model itself is dual to the 3D classical Ising model(Eq. (9) is also known as the 3D Ising gauge theory [8]).Using high-accuracy results of Ref. [12] for the criticalpoint and the duality relation λ pl = − / J/T ),where J is the Ising exchange coupling, we obtain λ ( c )pl =0 . λ bond and λ pl the model is self-dual [13], i.e. it maps to itself under the coupling con-stant transformation λ bond , pl → − / λ pl , bond ).This means that the phase diagram has a symmetry, or self-duality , line defined by λ bond = − / λ pl ).Under the duality mapping ( λ bond = 0 , λ pl =0 . → ( λ bond = 0 . , λ pl = ∞ ), which givesus two Ising-type critical points on the phase diagram.To calculate the rest of the phase diagram we per-formed Monte Carlo simulations using standard single-spin flip updates, supplemented by rare (once per N up-dates) flips of all spins belonging to bonds cut by planesoriented along any one of the crystal axes, or along any ofthe diagonals to these axes. There are 9 N possible planessatisfying this condition, and we select any of them atrandom. The plaquette energy (second term in (9)) isconserved by this update. To determine the 2nd-ordercritical lines, we employed a standard finite-size scalinganalysis of the specific heat C v , for linear system sizes N = 24 , , , and 60 (ie., for 3 N spins). First-ordercritical points were identified and located using energydistributions. These distributions are bi-modal (have twomaxima) for the first-order transitions and single-modalotherwise. We thermalized our samples for up to 10 MCsweeps (one sweep having 3 N elementary updates). Thedata were accumulated for ∼ × MC sweeps.The resulting phase diagram is presented in Fig.1. Thefirst-order transition coinciding with the self-duality linewas observed for 0 . > λ bond > . N = 60. Theinset of Fig.2 shows the evolution of the energy distri-bution function along the self-dual line. Even when thebi-modal structure is observed it is extremely weak, de-veloping only for large N , and the distribution can besampled in the minimum without flat-histogram or sim-ilar reweighting techniques. l bond l pl l bond l pl (I)(II) (III) (b) (a)(1)(3)(2) FIG. 1: (color online). The phase diagram of the Hamilto-nian (9). Circles correspond to the second-order transitions(open and filled symbols are related by the duality transfor-mation). Filled squares describe the first-order self-dual tran-sitions. Bold and dashed lines are used to guide an eye andcorrespond to the 1st- and 2nd-order transitions, respectively.The phases are: (I) - topological phase; (II) - topologically dis-ordered phase; (III) - magnetically ordered phase. In inset (a)we show the region where all phases meet each other. In inset(b) we show three alternative ways of connecting the lines.
As noted above, these results conflict with previous MC simulations in Ref. [10], who suggest the 1st-orderline splits into two 1st-order lines. The inset (a) of Fig.1shows a closeup of the controversial region. Though wewere able to resolve critical points with an accuracy of atleast three digits, we observed no splitting of the self-dual1st-order line into two 1st-order transitions. We also findno evidence for tri-critical points on the Ising-type linesas long as we can resolve two separate transitions. Thereremains a tiny parameter range between the apparentdisappearance of the bi-modal distribution on the self-dual line (this disappearance probably due to our limitedsystem size) and two resolved 2nd-order transitions.
N=48 N=600.22620.22660.22700.22740.0000.0050.0100.0150.0200.0250.0300.035 | E m a x ( ) - E m a x ( ) | -0.82 -0.80 -0.78 -0.76Energy0.000.040.080.12 N u m be r l bond=0.22700.226350.2266Energy distribution, N=60E max(1) E max(2) along self-duality line l bond FIG. 2: (color online). The distance between two maximain bi-modal energy distributions along the self-duality line for N = 60 as a function of λ bond . The inset shows examples ofthe energy distributions at various values of λ bond . To probe the behavior in this tiny parameter range weneed a different approach. We therefore scanned energydistributions at 30 points ( N = 48) along the line perpen-dicular to the self-duality line right in the questionable re-gion (short solid line in the inset (a)). If the 1st-order linewere to split above the scan, the third maximum wouldhave to emerge in the energy distribution right betweenthe two maxima we observe on the self-dual line - imply-ing that the energy maxima on the self-dual line could notmerge smoothly, and right below the split, three maximawould have to be seen in the energy distribution. How-ever all distributions along the scan were found to haveonly one peak. It is also clear from the main part ofFig.2 that on the self-duality line, the energy maximaapproach each other and merge continuously as λ bond in-creases. The curves presented in Fig.2 follow a powerlaw near the vanishing point, with corresponding criticalexponent ∼ . before or after the point where two 2nd-order lines touch the self-dualline (cases (2) and (3) in the inset (b), Fig.1). Unfor-tunately our data cannot distinguish between the alter-natives because the 2nd-order lines seem to touch at ex-tremely small (possibly zero) angle. Formally, option (2)fits the data best. Theoretically, the last two scenariosare less demanding since they fit the existing theory ofphase transitions (our data suggest that the 2nd-ordertransitions cannot merge into a single smooth curve andform a kink at the self-dual line). We are not aware ofany effective theory leading to the first scenario. Phases.
Using the two correspondence equations˜ h z = −
12 ln tanh( ˜ J x ) = λ bond ; (10a)˜ J z = −
12 ln tanh(˜ h x ) = λ pl , (10b)we can reformulate the phase diagram Fig.1 in terms ofthe renormalized parameters ˜ J x , ˜ J z , ˜ h x and ˜ h z of theTCM. The resulting phase diagram in terms of the ex-ternal fields is presented in Fig.3. Let us go through thephases in this Figure. hx hz self-duality line (I) (III) ~~ (II) FIG. 3: (color online). The phase diagram Fig.1 in terms ofthe renormalized external fields, ˜ h x = h x ∆ τ and ˜ h z = h z ∆ τ ,of the model (2). The phases (I), (II) and (II) are the sameas in Fig.1. The phase (I) corresponds to the topological phase ofthe model (2) (the “free charge” phase of the isotropic Z Higgs model (HM), [9]). In this phase the systemtends to have all B p = 1 and A s = 1 and a realizationof such a state is obviously not unique. The plaque-ttes with B p = − A s = − z -direction. However, h σ z i alsohas nonzero value everywhere in the phase diagram. The true order parameter may be written as h µ z i using thegauge-symmetrized Hamiltonian (4). A non-zero valueof this parameter results in the confinement of magneticvortices (no free vortices) and the condensation of elec-tric charges. In the HM this is the “Higgs” phase. Thephase (II) is characterized by a dual order parameterrelated to h σ x i , which can be defined by rewriting theHamiltonian in different variables. Its nonzero value re-sults in non-conservation of total magnetic “charge” andcondensation of magnetic vortices while electric chargesare confined (no free charges). This phase correspondsto the “confinement” phase of the HM. The transitionbetween the phases (II) and (III) is accompanied by asharp change in the number of vortices and charges, cor-responding to a “liquid-gas” type transition. The self-duality symmetry reflects the symmetry between chargesand vortices. Summary.
The topological phase of the toric codemodel (the “free charge” phase of the 3d gauge Higgsmodel) remains stable in a rather wide range of fieldsand breaks down via two Ising type transitions whosecritical lines meet with the 1st-order one correspondingto a liquid-gas type transition. The 1st-order line eithermeets with two 2nd-order lines in one multicritical point,or terminates before or after the point where two 2nd-order lines touch the self-duality line. The constructionof an effective field theory for this multicritical region isan interesting open problem.We thank E. Fradkin, B. Svistunov, S. Trebst, M.Troyer, I. Affleck, and K. Shtengel for discussions. Weare also indebted to M. Berciu and J. Heyl whose researchclusters were used to perform our MC simulations. [1] P.W.Anderson,
Science , 1196, (1987).[2] N. Read, B. Chakraborty, Phys. Rev. B , 7133 (1989).[3] A.Yu. Kitaev, Ann. Phys. (N.Y.) 303, 2 (2003).[4] G. Misguich, D. Serban, V. Pasquier, Phys. Rev. Lett. , 137202 (2002); cond-mat/0204428 .[5] N. Read, S. Sachdev, Phys. Rev. Lett. , 1773 (1991).[6] R. Moessner, S. L. Sondhi, Phys. Rev. Lett. , 1881(2001); cond-mat/0007378 .[7] S. Trebst, P. Werner, M. Troyer, K. Shtengel, and Ch.Nayak, PRL , 070602 (2007).[8] F.J. Wegner, J. of Math. Phys. , 2259 (1971).[9] E. Fradkin and S. Shenker, Phys. Rev. D , 3682 (1979).[10] G.A. Jongeward, J.D. Stack and C. Jayaprakash, Phys.Rev. D , 3360 (1980).[11] M. Suzuki, Progress on Theor. Phys. , 1454 (1976).[12] R. Gupta and P. Tamayo, The critical exponents for the3d Ising model , US-Japan Bilateral Seminar - Maui, Au-gust 28-31, 1996.[13] R. Balian, J.M. Drouffe, and C. Itzykson, Phys. Rev. D11