The modular S-matrix as order parameter for topological phase transitions

We study topological phase transitions in discrete gauge theories in two spatial dimensions induced by the formation of a Bose condensate. We analyse a general class of euclidean lattice actions for these theories which contain one coupling constant for each conjugacy class of the gauge group. To probe the phase structure we use a complete set of open and closed anyonic string operators. The open strings allow one to determine the particle content of the condensate, whereas the closed strings enable us to determine the matrix elements of the modular $S$-matrix, also in the broken phase. From the measured broken $S$-matrix we may read off the sectors that split or get identified in the broken phase, as well as the sectors that are confined. In this sense the modular $S$-matrix can be employed as a matrix valued non-local order parameter from which the low-energy effective theories that occur in different regions of parameter space can be fully determined. To verify our predictions we studied a non-abelian anyon model based on the quaternion group $H=\bar{D_2}$ of order eight by Monte Carlo simulation. We probe part of the phase diagram for the pure gauge theory and find a variety of phases with magnetic condensates leading to various forms of (partial) confinement in complete agreement with the algebraic breaking analysis. Also the order of various transitions is established.


Introduction
The study of hidden symmetries and phase structure of gauge theories has a rich history. A central motivation in the early works was to understand the confinement phenomenon in non-abelian gauge theories. An important step forward was the lattice action formulation by Wilson [1] and its Hamiltonian version by Kogut and Susskind [2], which allowed for an expansion around the strong coupling limit where confinement is manifest. The question was whether the confinement regime would extend all the way down to zero coupling. Wilson also introduced his celebrated loop operator as a diagnostic for confinement in the pure gauge theory.
Not long thereafter, 't Hooft and Mandelstam [3] suggested confinement in 3+1 dimensions to be a consequence of the dual Meissner effect, caused by magnetic disorder, notably a condensate of topological degrees of freedom, be it monopoles or fluxes. 't Hooft explained the (2 + 1)-dimensional version of confinement in SU (N ) gauge theories as a consequence of the condensation of Z N fluxes corresponding to the centre elements of the gauge group [4]. Polyakov, on the other hand, proved confinement of (2 + 1)-dimensional compact QED due to monopoles [5] and showed furthermore that the finite-temperature deconfinement transition in d = 3 + 1 is due to a Wilson line (string) condensate [6]. The interpretation of such transitions is related to the spectrum of topological defects in the various Higgs phases. For example for U (1) in 3+1 dimensions the occurrence of a phase transition depends on the compactness of the group, i.e. on the presence of monopoles. For non-abelian groups the situation depends on the representation of the Higgs field. In the case of SU (2) for example, one may or may not find a phase transition between the confined and Higgs phases depending 3 on whether the breaking is achieved with two adjoints or a fundamental representation. Shenker and Fradkin [7] constructed the phase diagrams for the SU(2) (lattice) gauge theory for the distinct cases with Higgs fields in the adjoint and the fundamental representation, and proved the absence of a transition between the confined and the Higgs phase in the latter. The difference can indeed be traced back to the topological structure: in the broken phase with a fundamental Higgs no defects exist because the group is broken completely. In the broken phase with one isovector, there is a residual U (1) and one has monopoles, and these will be confined if one adds another generic isovector. However, with the two vectors there is still the Z 2 centre of the group which is unbroken, which implies the existence of Z 2 fluxes in the broken phase.
Lattice formulations of the Z 2 theories date back to the early 1970s [8,9] and were also studied in [7], where the phase diagram of the Z 2 gauge theory with a matter field is computed in detail. A Hamiltonian formulation of Z N gauge theories and their duality properties in two and three dimensions was made in [10], where the same algebraic structure (1) was constructed out of generalized Ising-type lattice models. All authors agree that the transition to the confined phase is due to the condensation of magnetic fluxes. The analysis is facilitated by the introduction of nonlocal order and disorder parameters such as the Wilson and 't Hooft loops and the algebra they form, from which the phase structure could be understood qualitatively. As these operators are nonlocal, the algebraic structure exhibits the nontrivial braid properties that characterize the phase. For Z N theories one obtains an underlying Z el N ⊗ Z mag N symmetry of which the electric, magnetic and dyonic (anyonic) charges form irreducible representations (n, m) labelled by a pair of integers mod N (m, n = 0, 1, 2, . . . , N − 1). If one defines the loop operators (n,m) (C), then 't Hooft derived the crucial algebraic relation [3]: (n,m) (C) (n ,m ) (C ) = (n ,m ) (C ) (n,m) (C) e 2πi(nm +mn )k/N , where k ∈ Z is the linking number of the loops C and C . These loops can be interpreted as linked timelike loops of worldlines obtained after creating and subsequent annihilation of two particle-antiparticle pairs. But they can also be interpreted as two intersecting space-like loops on a torus winding around the two different cycles, one being a Wilson line and the other a magnetic Dirac string. On the other hand, one may think of equation (1) as specifying the braiding relations for the representation theory of a Z N × Z N algebra endowed with a (unique) nontrivial braiding structure that exactly belongs to the quantum double algebra of Z N denoted as D(Z N ). This operation on the torus is a topologically nontrivial vacuum-tovacuum transition (because no charges or fluxes are left) and that shows that there is a periodic vacuum structure in the theory, leading to an N 2 -fold vacuum degeneracy on the torus. This degeneracy in turn is equal to the number of particle species, i.e. the total number of unitary irreducible representations of D(Z N ), which for the case at hand is equal to N 2 . Similar algebras have surfaced in the study of emergent gauge theories from closely related lattice models such as Kitaev's toric code [11] and the Levin-Wen-type spin models [12], where furthermore the vacuum degeneracy (of a gapped phase) on the torus is proposed as a criterion for topological order.
It was also early on pointed out that 'the' Higgs phase structure of a gauge theory, where no massless gauge degrees of freedom (no continuous groups) survive, may be very rich by itself, allowing for distinct phases which exhibited different spectra of purely topological degrees of freedom [13]. By breaking with sufficiently high-dimensional representations of SU(2) for example, one finds that not just the centre, but any discrete subgroup of SU (2) can be selected to survive, also non-abelian groups like D N , or the tetrahedral group, etc. In 4 the corresponding phases, one finds a rich variety of non-abelian fluxes that exhibits highly nontrivial braiding and fusion properties with charges and also among themselves. The effective low-energy theory in such a Higgs phase is also called a discrete gauge theory (DGT); these are indeed gapped phases with only topological degrees of freedom, and therefore by definition phases with topological order. The full (2 + 1)-dimensional description of all the anyonic sectors of discrete gauge theories was given in [14] and revealed an underlying hidden quantum group structure, i.e. the excitations of a DGT with gauge group H were shown to correspond to the irreducible unitary representations of the quantum double D(H ) of H as defined by Roche et al [15]. The most interesting aspect of this perspective is that the ordinary 'electric charge' and topological 'magnetic flux' degrees of freedom are both part of the same representation theory, the magnetic degrees of freedom labelled by the conjugacy classes of H and the electric parts by the representations of the centralizer group in H of the given magnetic flux. This leads to a rather intricate interdependence of electric and magnetic sectors in case the goup H is nonabelian and a different way to look at order and disorder operators. We return to D(H ) and its representations in section 2.2. An important consequence of this underlying symmetry is that it suggested that in the pure DGT a corresponding set of gauge invariant loop operators should exist, thereby generalizing the Wilson/'t Hooft operators to all anyonic species, and indeed such operators have been constructed [16]. With these operators in hand, one can now probe the full phase structure of discrete gauge theories in a lattice formulation and that is what will be done in this paper.
The hidden quantum group symmetry and its representation theory that forms a modular tensor category is a powerful way to characterize the distinct topological phases of pure DGTs and forms the natural basis for the concept of topological symmetry breaking, which is defined as the breaking of quantum group symmetry [17][18][19]. A brief review of this symmetry breaking mechanism is given in section 1.3. In this paper, we study the validity of this mechanism by probing the various ground states through measuring the expectation values of single as well as linked anyonic loop operators, as is explained in section 1.4. In recent years we have witnessed a growing interest in systems that allow for the realization of different topological phases; the examples can be found in the Levin-Wen models [20], the Kitaev honeycomb model [21], discrete gauge theories [17] and, last but not least, in quantum Hall systems [22][23][24][25][26]. Most of these can be understood from the point of view of topological symmetry breaking involving the formation of a Bose condensate. In most instances, a Hamiltonian framework is used, but in this paper we show that to analyse such systems it may be profitable to switch to a euclidean action formulation that allows for the use of Monte Carlo (MC) simulations to determine the phase structure. In this approach, it is easy to directly measure the modular S-matrix in any phase of the system, which explains why we call this matrix an order parameter of such a system. After briefly recalling the basic ingredients of topological order and topological symmetry breaking and settling some notation, we show how open string operators and the modular S-matrix can be used as order parameters and phase indicators for topological symmetry breaking. We then introduce a class of multi-parameter lattice actions for non-abelian discrete gauge theories and verify the theoretical analysis [18,19] in detail by numerical simulations.

Topological quantum field theory basics
In this section, we set the stage and fix the notation for the rest of this paper. We study phases of systems that are described by a topological quantum field theory (TQFT) in 2 + 1 dimensions. 5 We label the different sectors or (anyonic) particle species by a, b, c, . . .. The two interactions between two particles in a TQFT are fusion and braiding.
Fusion. We describe fusion by the rule where the integer multiplicities N ab c give the number of times c appears in the fusion product of a and b. The fusion algebra is associative and commutative, and has a unique identity element denoted as '1' that represents the vacuum. Each sector a has a unique conjugateā (representing the corresponding anti-anyon) with the property that their fusion product contains the identity: Braiding. The particles in a (2 + 1)-dimensional TQFT can have fractional spin and statistics. Rotating a particle a by 2π (also called twisting) multiplies the state vector by a phase equal to the spin factor θ a |a twist → θ a |a , generalizing the usual +1 (−1) known from bosons (fermions) in 3 + 1 dimensions. Adiabatically moving a particle a around another particle b in a channel c is called a braiding and can have a nontrivial effect on the state vector of the system, given by θ c /θ a θ b .
Quantum dimensions. The quantum dimensions d a of particle species a are another set of important quantities in a TQFT. These numbers satisfy the fusion rules (2), i.e. d a d b = c N ab c d c . The quantum dimension of an anyonic species is a measure for the effective number of degrees of freedom, corresponding to the internal Hilbert space of the corresponding particle type. The Hilbert space dimension of a system with N identical particles of type a grows as (d a ) N for N large. In general, the quantum dimensions d a will be real numbers; however, for DGTs they are integers. The total quantum dimension D of the theory is given by and the topological entanglement entropy of the ground state is proportional to log D.

Diagrammatics.
There is a powerful diagrammatic language to express the equations describing the TQFT, which we will use to relate the values of observables as they can be measured in the different phases. In this paper, we will use the notation and definitions given by Bonderson [27]. Particle species are represented by lines, fusion and splitting by vertices. A twist is represented by a left or right twist on a particle line: The evaluation of simple diagrams is rather straightforward, and complicated diagrams can be simplified using braid relations and the so-called F symbols which follow from associativity 6 of the fusion algebra. The simplest examples are the closed loop of type a that evaluates to the quantum dimension d a : whereas the twisted loop is equal to d a θ a : Of particular interest are the generators of the modular group, S ab and T ab = e −2πi(c/24) θ a δ a,b , where the c is the central charge of the theory, not to be confused with a particle type. As we mentioned before, the importance of the rather abstract diagrammatic notation is that the diagrams directly correspond to observables in our euclidean lattice gauge theory formulation. In the euclidean three-dimensional (3D) formulation of topological theories the values these diagrams have correspond to the vacuum expectation values of the corresponding anyon loop operators; for example, in the unbroken phase one may measure where the left-hand side (lhs) is now defined as the value of the path integral with the nonlocal loop operator for particle species a inserted and the right-hand side (rhs) is obtained if we are probing the system in the unbroken phase governed with the ground state denoted as 0 and governed by the algebra A. We use the subscript 0 because the value of the same diagram may be different if it is evaluated in a different phase with a ground state that we will denote by ; in the remainder of the paper we will therefore always use brackets with a subscript.

Topological symmetry breaking
In this section, we briefly recall topological symmetry breaking, the phenomenon that a phase transition to another topological phase occurs due to a Bose condensate [17,18]. The analogy with ordinary symmetry breaking is clear if one thinks of the particle as representations of some quantum group, and assumes that a bosonic degree of freedom i.e. with θ c = 1-fundamental or composite-condenses. The breaking can then be analyzed, either from the quantum group (Hopf algebra) point of view or from the dual or representation theory point of view [19]. Let us illustrate this with an example of ordinary group breaking. Suppose we have a gauge group SU (3) and a Higgs triplet that acquires a vacuum expectation value = (1, 0, 0), then the SU(2) subgroup working on the last two entries will leave invariant. Equivalently, this SU(2) subgroup may be characterized by the way the SU(3) triplet decomposes under the SU(2) action as 3 → 2 + 1 where the singlet on the right corresponds exactly to the new SU(2) invariant ground state. In that sense one may select a specific residual gauge symmetry by choosing an appropriate Higgs representation which has a singlet under that residual group in 7 its branching. For example, if we want to break an SU(3) group to the SO(3) subgroup which is characterized by the branching rule 3 (as well as3) → 3, then we may choose the Higgs field to be in the 6D irrep of SU(3), because then 3 × 3 =3 + 6 → 3 × 3 = 1 + 3 + 5, from which follows that 6 → 5 + 1 and again the singlet on the right corresponds to the SO(3) invariant vacuum state .
In the case of general quantum groups it is this branching rule approach which is most natural and powerful in the context of TQFT because the fusion algebra corresponds to the representation ring of the quantum group. A general treatment with ample examples can be found in [19]. Let us point out some essential features of this procedure that one has to keep in mind. As the quantum group centralizes the chiral algebra in the operator algebra of a CFT, one expects that reducing the quantum group will correspond to enlarging the chiral algebra, and this turns out to be the case. In contrast to ordinary group breaking, the topological symmetry breaking procedure involves two steps: firstly, the condensate reduces the unbroken fusion algebra (also called a braided modular tensor category) A to an intermediate algebra denoted by T . This algebra, however, may contain representations that braid nontrivially with the condensed state, i.e. with the new vacuum and if that is the case, these representations will be confined and will be expelled from the bulk to the boundary of the sample. Confinement implies that in the bulk only the unconfined sectors survive as particles and these are characterized by some subalgebra U ⊂ T . Let us briefly describe the two steps separately.
From A to T . Assuming that a certain bosonic irrep c will condense due to some underlying interaction in the system implies that c will be identified with the vacuum of T . For our purposes, a boson is a sector with trivial (integer) spin, although in fact in the context of 2 + 1 dimensions one has to also require that fusion of this field with itself has a channel with trivial braiding.
The definition of the new vacuum requires a redefinition of fields. Firstly, fields in A that appear in the orbit under fusion with the condensed field c are identified in T , so, if c × a = b then a, b → a . Secondly, if a field b forms a fixed point under fusion with the condensate c, then the field will split at least into two parts: b → i b i . The identifications and splittings of representations can be summarized by a rectangular matrix n t a that specifies the 'branching' or 'restriction' of fields a from A to T with fields t, r, s, . . .: This branching matrix is a rectangular matrix (the number of particle types in the A and T theories is not equal in general) of positive integers. We will also consider the transpose of this matrix denoted as n a t , which specifies the 'lift' of the fields t ∈ T to fields a ∈ A: t → a n a t a = a∈t a.
One may now derive the fusion rules T from the fusion algebra (2). Because of the identifications, it is often the case that the intermediate algebra T , despite being a consistent fusion algebra, is not necessarily braided; in more technical terms, it satisfies the 'pentagon' equation but not the 'hexagon' equation. The physical interpretation of this fact is that the sectors in T do not yet constitute the low-energy effective theory. This is so because sectors t that have an ambiguous spin factor, meaning that not all θ a of the lift a ∈ t are equal, will be connected to a domain wall and hence are confined in the new vacuum. The confined excitations will be expelled to the edges of the system or have to form hadronic composites that are not confined.
Yet the T algebra plays an important role: in [24] for example, it was shown that the T algebra governs the edge/interface degrees of freedom in the broken phase.
From T to U. Some of the sectors in T will survive in the bulk and some will be confined. The physical mechanism behind confinement in (2 + 1)-dimensional topological field theories is nontrivial braiding with the condensate. The vacuum state or order parameter should be single valued if carried adiabatically around a localized particle like excitation. If it is not single valued that would lead to a physical string or 'domain wall' extending from the particle that carries a constant energy per unit length. The unconfined algebra U consists of the representations in T minus the confined ones; it is this algebra that governs the low-energy effective bulk theory. The confined representations can be determined in the following way. First we define the 'lift' of a representation in T as the set of representations b ∈ A that are restricted to t. Now, if all of the representations in the lift of t braid trivially with the lift of the vacuum, the sector t is part of U. Otherwise, it is confined. One may prove that the U algebra closes on itself with consistent fusion rules, while consistent braiding is achieved by assigning the (identical) spin factors of the parent sectors of the unbroken theory to the U fields.
Let us finally mention a useful quantity, the so-called quantum embedding index q defined in [28]; it is a real number characterizing the topological symmetry breaking. This quantity is defined as where the index a runs over the sectors of the unbroken phase A, which correspond to the lift of any sector u or t of the algebra U or T ; the n a u is the lift of sectors u to their parents a and d a is the quantum dimension of the representation a. Observe that this expression is independent of the particular sector u, which is a non-trivial result explained in [28].
Choosing for u the new vacuum, we have d u = 1 and obtain that q just equals the total quantum dimension of the lift of the U (or T ) vacuum in the unbroken A theory. The quantum embedding index is the analogue for the embedding index defined by Dynkin for the embedding of ordinary groups. As an aside we mention that the change in topological entanglement entropy of the disc changes also by log(D A /D U ) = log q in a transition from an A to a U phase [28].
Let us, to conclude this subsection on topological symmetry breaking, illustrate the procedure with a very straightforward example, namely the breaking of the quantum group A = SU(2) 4 . It has five irreps labelled by = 0, . . . , 4 with spin factors θ a = 1, 1 8 , 3 8 , 5 8 , 1. The = 4 is the only boson and we assume it to condense. The lift of the new vacuum corresponds to the = 0 + 4 of A, and hence the embedding index q = d 0 + d 4 = 1 + 1 = 2. The 1 and 3 reps of A are identified, but because they have different spin factors, the corresponding T representation will be confined. In U we are therefore left with the = 2 rep., which splits because it is a fixed point under fusion with the condensate as 4 × 2 = 2. We write 2 → 2 1 + 2 2 . The values for the spin and the quantum dimensions and the fusion rules for these representations fully determine the unconfined quantum group to be U = SU(3) 1 . We recall that the nomenclature of the groups is linked to the chiral algebra; it is therefore not surprising that the SU(2) 4 quantum group breaks to the smaller quantum group SU(3) 1 , which is related to a larger chiral algebra. For the chiral algebras one has the conjugate embedding SU(2) 4 ⊂ SU(3) 1 which is a conformal embedding. This conformal embedding in turn is induced by the SO(3) ⊂ SU(3) embedding mentioned at the beginning of this subsection.

Observables
Our objective is to verify the theoretical predictions of the topological symmetry breaking scheme in a class of euclidean gauge theories that are expected to exhibit transitions between different topological phases. We will numerically evaluate the expectation values of various topological diagrams using MC simulations, and in this section we calculate the predicted outcomes of a variety of possible measurements from theory. The strategy has two steps: (i) determination of the condensate (including the measurement of the embedding index q) by evaluating the basic nonlocal open string order parameters, given by equation (33); (ii) measuring the so-called broken modular S-matrix and from that construct the S-matrix of the U phase. We also will see that the condensate fixes the branching and lift matrices and having determined those we can also predict the outcome of measurements of other topological diagrams corresponding to the lifts of U fields to A fields.

Determination of the condensate and the embedding index q.
We measure the open string operators in the model. Note that in our pictorial representation time flows upward, so a vertical line physically represents the creation, propagation and annihilation of a single particle. For the particular case of a DGT, which we study in this work, these lines have a realization as operators on a spacetime lattice, see equation (33).
If the symmetry is unbroken we will have for any nontrivial field a that because the diagram represents the creation and subsequent annihilation of a single a-particle. However, in the broken situation the expectation value will be nonzero for all fields φ i ∈ A in the condensate, which we denote by . So writing we obtain that, in general, This in turn implies that it is simple to measure q as

Determination of confinement and other topological data of the broken phase.
Once we have determined the components of the vacuum we can determine the lifts of the t fields, simply by studying the fusion rules of × a = t , where t denotes the lifts of those t fields which contain a, i.e. for which n a t = 1. Having obtained the lifts of the t fields the next step is to make the measurement determining whether a given t field is confined. This involves the measurement of the index η, or simply: Alternatively, one can measure certain closed a i loop operators that are also defined for fields a that split under branching and that will be defined later, for which it holds that: It follows that, from these measurements, the fields that are confined can be determined, but also the quantum dimensions d u and twists θ u of the unbroken U theory are obtained.

The broken modular S-and T-matrices.
Instead of the fusion coefficients N ab c , an alternative specification of a (modular) topological field theory is by its representation of the modular group S L(2, Z) generated by the Sand T-matrices with C being the charge conjugation matrix. The corresponding matrix elements can be expressed in the fusion coefficients and spin factors: where D is the total quantum dimension and the constant c is the conformal central charge of the corresponding conformal field theory. We recall that the central charge of a DGT is zero, so in that case the T -matrix is just the diagonal matrix containing the spin factors. The great advantage of switching to the modular data is that unlike the fusion coefficients these generators can be directly measured using the anyon loop operators that arise naturally in a 3D euclidean formulation of the theory. We will evaluate the expectation value of these S-matrices numerically in our lattice formulation of multiparameter discrete gauge theories later on. The measured Sand T -matrix elements do not satisfy the relations (15) directly; however, using the measurements the full Sand T -matrices of the U theory, which do satisfy the modular group relations, can be constructed. In the unbroken theory, the measured S-matrix elements S ab correspond to the expectation values of the Hopf link with one loop coloured with representation a and the other with representation b: where S ab is the S-matrix of the unbroken A theory. We can, however, also determine the modular S-matrix of the residual U theory S uv directly from measurements if we take the splittings of certain fields a ⇒ {a i } into account appropriately. We will show how to do this later for the DGTs in detail and give a more general mathematical treatment of this elsewhere [29]. Then we will arrive at an explicit formula and algorithm to determine S uv : This expression involves not only the branching (lift) matrix n a i u , but also what we will call the broken S-matrix defined asS a i b j = S a i b j , which, because of the splitting, clearly involves a larger size matrix than the modular S-matrix of the original A phase. From the broken S-matrix we may directly read off S uv , the S-matrix of the effective low-energy TQFT governed by U. An important observation is that the values of the S-matrix elements in a broken phase will be different from the ones in the unbroken phase, for example because of the contribution of the vacuum exchange diagramS depicted below, in which the condensed particle is exchanged giving a nonzero contribution in the broken phase, while it would give a vanishing contribution in the unbroken phase: In the explicit calculations later on we show that this vacuum exchange diagram leads to a change in the S-matrix which depends on the subscripts introduced above. It turns out that it is also possible to calculate the broken S-matrix from first principles; this will be discussed in a forthcoming paper [29]. As is to be expected, one finds identical rows and columns in the broken S-matrix, for components that are identified, whereas the entries for confined fields will be zero. With this prescription the formalism outlined above is applicable in any phase of the theory including the unbroken one where there is no splitting and the vacuum exchange diagram gives a vanishing contribution. The measured T -matrix, on the other hand, is given by again with T ab 0 = T ab . After measuring or calculating the Sand T -matrices in a given phase, we can reconstruct the fusion coefficients with the help of the Verlinde formula [30] To conclude, we have in this section summarized the basic features of a TQFT and considered some aspects of topological phase transitions induced by a Bose condensate; furthermore, we explained how the measurement of the L-, Sand T -operators in the broken phase fully determine the quantum group of a (broken) topological phase. The general scheme to analyse the breaking pattern of some multiparameter TQFT is to first use the open string operators to probe which fields are condensed in the various regions of parameter space. In a given broken phase we can subsequently compute/measure what we will call the broken S-matrixS a i b j , where, as mentioned, the subscript labels the splitting of the corresponding A field. From the broken S-matrix, we can read off the S-matrix of the U theory. In the remainder of the paper, we will explicitly execute this program for discrete gauge theories.

Z 2 gauge theory and topological order: a prelude
Before applying our approach to (non-abelian) discrete gauge theories (DGTs) in general, let us make some connections to previous work from different perspectives. After Wilson's seminal work on euclidean lattice gauge theory for non-abelian Lie groups to study the confinement of quarks [1], a Hamiltonian formalism for the same problem was soon developed [2]. It was clear that because the gauge fields now take values in the group instead of the Lie algebra, one could also study models based on a finite group. These models, which we call DGTs, were mostly studied as approximations to U(1) or SU(N ) theories in times when computers were not as powerful as they are today. They are, however, also interesting in their own right, since they are purely topological: there are no local degrees of freedom, and only the topological (generalized Aharonov-Bohm-type) interactions survive.
This does not automatically mean that all observables are topological quantities: the appearance of virtual flux-antiflux pairs gives small size-dependent corrections to the loop-like observables in these theories. However, since these excitations are gapped, these corrections are exponentially small. We will show below by explicit calculation that as long as one stays away from the critical points, it is justified to think of the observables in these theories as topological quantities.

Hamiltonian formalism.
To connect with other work on topologically ordered systems, let us first go to a Hamiltonian formalism. This is formally done by taking a time slice of the spacetime lattice and taking the limit in which the temporal spacing goes to zero [31]. The Hamiltonian of (2+1)-dimensional Z 2 gauge theory on a square spatial lattice is where the operators P l and Q l act on links, the second term is a sum over the elementary plaquettes of the lattice where p1, . . . , p4 are the links of a single plaquette and λ is the coupling constant. The operators satisfy {Q l , P l } = 0, P 2 l = Q 2 l = 1, which means that a possible representation can be given in terms of the Pauli matrices P l = σ 3 , Q l = σ 1 acting on spin-1 2 bosons living on the links. Note that the algebra above is the same as the Z 2 version of (1), and indeed a closed string of P i operators generates a Wilson loop, whereas a closed string of Q i operators creates a closed Dirac string. Gauge transformations act on the star of four links i1, . . . , i4 adjacent to a site i G i = P i1 P i2 P i3 P i4 , and to build the gauge invariant Hilbert space, one has to implement a Gauss law for physical states |ψ Now we can make the connection with the work by Kitaev [32] and Hastings and Wen [33]. Their models (Toric code, Z 2 string nets) correspond to Hamiltonian Z 2 DGT where the coupling λ = 0 and the gauge constraint (21) is not strictly enforced. Setting λ = 0 makes the theory purely topological, the ground state is an equal weight superposition of all states where C is a closed loop of links and |0 is the state with the property P l |0 = |0 for all links l.
Viewing the link variables as spin-1 2 bosons, this vacuum state corresponds to all the spins being in the up state. Since the expectation value of any loop operator (22) in the ground state is equal to one, and these loops are Wilson loops in the gauge theory language, the theory is topological.
By not enforcing the gauge constraint (21) strictly but adding it as a term to the Hamiltonian, these models allow for massive open strings. Such open strings are not gaugeinvariant at their endpoints and therefore correspond to external charges.

Euclidean formalism.
The Z 2 gauge theory in the euclidean approach, where we discretize both space and time, is described by the action where the sum is again over all plaquettes (now both spatial and temporal) and the U variables are numbers ±1. The partition sum are the quantities of interest here. The gauge invariance, which in the Hamiltonian formulation was enforced by projecting out states from the Hilbert space, is now manifest in the action and the operators. The partition sum is over all gauge field configurations, but since all sums are finite, gauge fixing is not required 4 . If the coupling β is large, the dominant contribution from the partition sum will be from field configurations where all plaquettes UUUU = +1. In the limit β → ∞, this is strictly true, and one is left with a TQFT, as was the case for the Hamiltonian (20) with λ = 0. For β small there is a confining phase; the phase transition is at β = 0.7613 [34].
In most of this work, we study the topological properties of a DGT, for a general group H . To show that for a finite coupling constant β this is a good approximation, let us perturbatively calculate the expectation value of a Wilson loop in this Z 2 theory. The Wilson loop W (C) is the product of U variables around a closed loop C For large β, the action is minimized by configurations for which all plaquettes are +1. The first-order perturbation comes from those configurations in which one link is −1. In three dimensions, this excites four plaquettes, so the Boltzmann weight for such configurations is e −4β smaller than for those with no excited plaquettes.
If the lattice has size N × N × N , there are 3N 3 links. For a contour C of length L, This shows that the corrections to the purely topological result W (C) = 1 are, for β several times larger than the critical point, negligible for simulations of reasonable lattice sizes: for a Wilson loop size 10 × 10, β = 3.0 yields corrections only in the third digit.
Another gauge-invariant quantity is the 't Hooft loop, which lives on a loop C of the dual lattice. Such a loop pierces a number of plaquettes p, and the 't Hooft operator flips the sign of the coupling β → −β for these plaquettes. This forces a Z 2 magnetic flux through these plaquettes. We will define operators generalizing the 't Hooft and Wilson loops for general non-abelian DGTs shortly.

Phase structure.
The action (23) can realize three phases when one also allows for negative coupling. For large positive β, the phase mentioned before is realized, where almost all plaquettes are +1. For large negative β, almost all plaquettes are −1. For small |β|, a confining phase where the magnetic Z 2 flux has condensed is realized.
To study the phase diagram of a (non-abelian) DGT in full, we find it convenient to formulate the action in the class basis, instead of the irrep basis. This means that we do not take the character of the group element of the plaquette product UUUU in the action, but we define delta functions on each class. We will explain in detail how this works in section 2.3. For Z 2 the phase diagram is one dimensional, but the introduction of a second coupling constant will get rid of the need for negative couplings: where δ A (U ) for a group element U and a conjugacy class A gives +1 if U ∈ A and zero otherwise. For non-abelian groups this formulation makes the phase diagram much more intuitive; for Z 2 it is rather artificial.
In figure 1 the phase diagram of the pure Z 2 gauge theory is shown as a function of the conjugacy class couplings β +1 and β −1 . Later in this work, we present similar phase diagrams for theD 2 gauge theory.
It is well known that the inclusion of matter coupled to the gauge fields complicates the phase diagram strongly. The question of whether there exist good order parameters to distinguish the phases in coupled gauge-matter systems is interesting in its own right and highly nontrivial [35], but it is not something we will go into here.

DGT and the quantum double of a finite group
The particles in a DGT, their fusion and braiding properties, spins and so on are all obtained by working out the representation theory of the underlying quantum group, which is the quantum double of the finite discrete subgroup. We will not give a detailed account of the emergence of quantum group symmetry in DGT, as this can be found in the literature [14], but we do present a short summary of the basics to fix the notation and introduce some key concepts required later on.
Consider the following operators acting on states in the Hilbert space of a DGT. Firstly, there is the flux projection operator, denoted by P h , which acts on a state |ψ Secondly, we have the operator g, for each group element g ∈ H , which realizes a global gauge transformation by the element g: where it should be noted that we have not yet modded out by the gauge group to obtain the physical Hilbert space. These operators do not commute, and realize the algebra The set of combined flux projections and gauge transformations {P h g} h,g∈H generates the quantum double D(H ), which is a particular type of algebra called a Hopf algebra. The representation theory of the quantum double D(H ) of a finite group H was first worked out in [15], but here we follow the discussion presented in [36] and follow the conventions of those lecture notes.
Let A be a conjugacy class in H . We will label the elements within A as for a conjugacy class A of order k. In general, the centralizers for the different group elements within a conjugacy class are different, but they are isomorphic to one another. Let A N ⊂ H be the centralizer for the first group element in the conjugacy class A, denoted by A h 1 .
The set A X relates the different group elements within a conjugacy class to the first: This still leaves a lot of freedom, but we fix our convention such that A x h 1 = e, the group identity element. The centralizer A N , being a group, will have different irreps, which we label by α.
The vector space for a representation α is spanned by a basis α v j . The internal Hilbert space corresponding to an irrep of the quantum double that combines magnetic and electric degrees of freedom, V (A,α) , is then spanned by the set of vectors where i runs over the elements of the conjugacy class, i = 1, 2, . . . , dim A and j runs over the basis vectors of the carrier space of α, j = 1, 2, . . . , dim α. These irreducible representations correspond precisely to the particle types a, b, . . . in section 1.2. They obey a set of fusion rules as in (2) and it is possible to calculate the modular S-matrix, F-symbols and so on.
To see that this basis is a natural one to act on with our flux measurements and gauge transformations, consider the action of a pure flux projection P h e The matrix action (A,α) of an irreducible representation (A, α) of some combined projection and gauge transformation P h g: where the elementg is that part of the gauge transformation g that commutes with the flux A h 1 , defined asg and D α (·) i j is the matrix representation of the centralizer. To conclude, we give a simple expression for the modular S-matrix that can be obtained by calculating the trace of the monodromy matrix where [g, h] is the group theoretical commutator: [g, h] = ghg −1 h −1 .

Lattice actions and observables
We discretize the 3D spacetime into a set of sites i, j, . . . using a rectangular lattice. The gauge field U i j , which takes values in the gauge group H , lives on the links i j, jk, . . . connecting sets of neighbouring sites. The links are oriented in the sense that U i j = U −1 ji [1]. We note that the gauge field U i j takes care of the parallel transport of matter fields that are charged under the gauge group from site i to site j. An ordered product of links along a closed loop is gauge invariant up to conjugation by a group element and measures the holonomy of the gauge connection. Gauge transformations are labelled by a group element g i ∈ H and are performed at the sites of the lattice. The gauge field transforms as where the orientation of the links (incoming or outgoing) has to be taken into account as shown in figure 2. The standard form for the lattice gauge field action makes use of the ordered product of links around a plaquette i jkl: which transforms under conjugation by the gauge group, i . The gauge action per plaquette that corresponds to the Yang-Mills form F 2 µν in the continuum limit for H = SU(N ) is given by where χ α is the group character in irrep α and β α is inversely proportional to the square of the coupling constant for irrep α. This is known, for H = SU(2) and the sum over representations limited to the fundamental one, as the Wilson action [1]. The action (29) is the euclidean analogue of the Kitaev Hamiltonian introduced in [32] (which in itself is a variation of the Kogut-Susskind Hamiltonian [2]) where the electric field is represented by the timelike plaquettes and the gauge constraint (21) is not necessary since everything is manifestly gauge invariant.
For SU(N ) gauge theories, one usually only includes the fundamental representation and is thus left with only one coupling constant. This is not necessary however: gauge invariance of the action is ensured by the fact the characters are conjugacy class functions and therefore we will consider actions where the number of independent couplings is equal to the number of conjugacy classes i.e. the number of irreps (for a finite group, these numbers are finite and equal).
For our purposes, namely the study of magnetic condensates in DGTs, equation (29) is not the most convenient to work with. We perform a change of basis in the space of coupling constants to write it as a sum over delta functions on conjugacy classes: δ A (h) = 1 if h ∈ A, and 0 otherwise. In this basis the action becomes This formulation allows us, in particular, to control directly the mass of the different fluxes in the theory, which will ease the search for different vacua in the phase diagram. Increasing the coupling constant for a certain conjugacy class (magnetic flux) A will increase the contribution of configurations carrying many A fluxes to the path integral. Likewise, setting all β A to zero except β e , the coupling constant for the trivial conjugacy class will result in an 'empty' vacuum and therefore an unbroken phase.
To perform the transformation to the conjugacy class basis, we need to make use of the following orthogonality relations valid for all finite groups H : where |H | is the order of the group H , |A| is the order of the conjugacy class A, R is the set of irreps and group integration is defined as Equations (30) and (31) show that the irreducible representations of a group H form an orthonormal set for functions on conjugacy classes of H . We thus expect the conjugacy class delta function to be expressible in terms of characters for some set of constants {c α }. We multiply both sides of this expression by a character of the same group element in another irrep β and perform the integrations by the use of the orthogonality relations (30) and (31) where the slightly abusive notation χ α (A) means the character of any group element of A in the representation α. This shows that which in turn implies that the difference between (29) and (30) is just a change of basis: where C is the set of conjugacy classes and β A (β α ) is given by To probe the physics of the system for a fixed set of values of the coupling constants in the action, we will use a set of order parameters and phase indicators. These order parameters are in one-to-one correspondence with the set of fundamental anyonic excitations of the theory.
Order parameters and phase indicators. We distinguish two different sets of order parameters that are closely related to one another. The first is the set of closed loop operators, which physically correspond to the creation, propagation and annihilation of an anyon-anti-anyon pair in spacetime. The second is the set of open string operators that create, propagate and annihilate a single anyon. In the background of a trivial vacuum, only the loops can have nonzero expectation values, since the creation of a single particle would violate the conservation of the quantum numbers of the vacuum in such a background. This means that the open strings tell us something about possible Bose condensates, whereas the closed loops tell us about the behaviour of external particles put into this background. We define the open string operators only for the purely magnetic sectors, since in this work we only study magnetic condensates 5 . First we will define the loop operators. This set of nonlocal order parameters was introduced in a previous publication [16]. For a full discussion, see that work. Here we recall the essentials and fix the notation. The closed loops are a generalization of the Wilson and 't Hooft loops. They create a particle-antiparticle pair from the vacuum and annihilate them at a later time. These loops allow us to calculate the Aharonov-Bohm-type phases and determine which anyonic excitations will be confined. In SU(N ) gauge theories, the Wilson loop for a free excitation, e.g. in the Higgs phase of SU(2) theory, in general falls off as e −c P , with P being the perimeter of the loop, whereas a confined excitation, such as the 3 charge of an external quark source in pure SU(3) gauge theory describing QCD, falls off as e −c A , with A being the area of the loop.
Because the excitations in a DGT are gapped, numerically we find that the expectation values of loop operators are constant as a function of size. The argument for this behaviour for the Z 2 theory is in section 2.1.2. Although strictly true only in the limit of infinite coupling constant, the gap suppresses the dependence on size so strongly that we will assume that the theory is a purely topological one in the region of coupling constant space that we are interested in.  (33). The ordered product U p of links around a plaquette 'p' needs to be taken with an orientation that has to be constant throughout the loop.
Let us draw a closed loop on the dual lattice; this loop pierces a set of plaquettes C through which we will force magnetic flux. Now draw another loop, this time on the real lattice, such that (i) each point of this loop lies on the corner of a plaquette in C and (ii) the two loops do not link 6 . The combination of the two loops establishes a framing: we have selected a location for the electric charge of a flux-charge composite. This framing also provides us with a point and orientation on each plaquette from which to take the plaquette product U p : for non-abelian groups the product of the four links depends on which corner you start.
To insert a flux h in a plaquette p ∈ C, we have to 'twist' the Boltzmann factor of this particular plaquette by locally changing the action from S(U p ) to S(h −1 U p ): if the minimum of the action was previously obtained for U p = e, it is now shifted to U p = h. We want to perform this twisting procedure for all plaquettes in C. Figure 3 shows the first two plaquettes of such a loop C.
The notion of a group element as a magnetic flux is not gauge-invariant: under a gauge transformation by g, a flux h transforms as gh g −1 . Therefore it is necessary to sum over the group elements h in a conjugacy class A in some way.
One can go about this in two different, and inequivalent, ways.
• The authors of [37] only studied pure magnetic flux loops (without electric charge) and performed a sum over the conjugacy class for each plaquette in C individually. This leaves a gauge-invariant expression, but the loop loses its framing, since a conjugation by the element U 01 maps for a plaquette spanned by group elements U 01 , . . . , U 30 .
• When dealing with nontrivial braiding properties of loop operators, it is necessary to choose a basepoint i 0 in space with respect to which all operators are defined; it provides a calibration that serves as a 'flux bureau of standards', borrowing a term from [38]. This point can be anywhere in spacetime and does not need to be on the loop. We then define a function k i p (h, {U i j }) of h and the gauge field variables {U i j } for the twist element that has to be inserted into the plaquette product for plaquette p, where i p is the corner of the plaquette chosen in the framing With this notation and the above considerations, the anyonic operator (A,α) is given by 7 Here p j iterates over the plaquettes in C and D α is the representation function of the centralizer irrep α of A N . The link U j−1, j neighbours the plaquette p j , and the combination in brackets always takes values in the centralizer subgroup of the conjugacy class A. The exponential of the difference of two actions changes the minimal action configuration to one containing flux h for the plaquette under consideration. The operator in expression (33) is a generalization of the Wilson and 't Hooft loops, and by constructing it we have established the desired one to one correspondence between irreducible representations of the quantum group and loop operators for the pure DGT. If we fill in for A the trivial conjugacy class, the exponent vanishes and the x group elements are equal to the group unit, so after we multiply out the D α -matrices we are left with (e,α) (C) = χ α U 1,2 U 2,3 · · · U n−1,n U n,1 , where the product of U s is an ordered product along the loop on the lattice. On the other hand, if we replace α by the trivial representation, we are left with which is comparable to the order parameter proposed in [37], but the gauge invariance with respect to the transformations (28) is ensured in a different way. We sum over the conjugacy class only once and insert the flux in a gauge invariant way by parallel transporting it along the loop from a fixed basepoint. The operator in [37] sums over the conjugacy class for each individual plaquette. This way also gauge invariance is achieved, but the loop loses its framing and therefore is not suitable for describing true anyonic charges. The open magnetic string operators are a variant of expression (33) where the set of plaquettes C corresponds to an open string on the dual lattice. Looking at the h-forest configurations, it can immediately be seen that such a string, corresponding to the creation and subsequent annihilation of a single particle, has zero expectation value in the trivial vacuum. For these strings to acquire a nonzero expectation value, a vacuum exchange contribution is required, which we will focus on now.
The vacuum exchange contribution. We use the set of operators { (A,α) } to measure the elements of the S-matrix by picking two loops C 1 and C 2 that link each other once: In the trivial vacuum the S (A,α)(B,β) -matrix elements of fluxes g ∈ A and h ∈ B for which gh g −1 h −1 = [g, h] = e evaluate to zero (this is what we measure using the operators (33) and calculate algebraically (27)). If we, however, measure the S-matrix elements of such noncommuting fluxes in a broken vacuum, nonzero matrix elements can appear. This is most easily explained by considering an example. The main contribution to a single loop of pure magnetic flux is of the form depicted in figure 4. This configuration is called the h-forest state in the earlier work [37]. Modulo gauge transformations are the dominant configuration in the trivial vacuum that contributes to a loop of flux labelled by conjugacy class A, where g ∈ A. Expression (33) contains a sum over these group elements within a conjugacy class, but let us for now focus on one of the group elements. Each link in this configuration has the value e, except for the fat links in figure 4, which have the value h. That this configuration leads to a loop or tube of flux is easily seen: within the forest each plaquette has a value e h e h −1 = e, whereas at the edges the value is e h e e = h (depending on the orientation of the plaquette product). This is also the easiest way to see the origin of the Aharonov-Bohm effect on the lattice: an electric charge loop having linking number 1 with the flux loop will have exactly one link with the value h in it; therefore its value will be χ α (h).
Consider now the dominant configuration that contributes to the S-matrix element S (A,1)(B,1) . We again pick two group elements g ∈ A, h ∈ B and draw a similar diagram. This is shown in figure 5. By a similar logic this causes the plaquettes at the boundary of either forest to have the value g, respectively h. Inside the forests most plaquettes still have the value e; however, there are some plaquettes that are different. There is a tube of plaquettes that have the value [g, h], where the two forests intersect. In general, this group theoretical commutator is not equal to the identity element for non-abelian groups. This is the physical reason behind the appearance of zeros in the S-matrix for non-abelian theories. This tube of plaquettes represents a flux [g, h] going from one loop to the other. In the trivial vacuum, this flux will be gapped, so the contribution of this diagram to the path integral expectation value will be negligible.
However, a different situation appears when we are in a vacuum where the flux [g, h] has Bose condensed. We cannot give a single configuration that contributes dominantly to the path integral (there are many), but we can say that configurations like the one in figure 5 are now contributing since the mass for the flux [g, h] has disappeared. Thus we expect that in the measurements there will be cases when zeros in the original S-matrix will obtain a nonzero value in the broken phase.
An auxiliary A N gauge symmetry. The operators (33) are invariant with respect to the local H gauge transformations (28). However, in our formulation of the operators we have tacitly introduced another, auxiliary A N gauge symmetry that is less obvious. A crucial property that allows one to determine the topological symmetry breaking pattern in detail is that the loop operators do transform nontrivially under this symmetry. In a nontrivial ground state, these symmetries may be broken and will therefore lead to the lifting of certain degeneracies related to the splitting of fields in the topological symmetry breaking process. So this hidden symmetry turns out to be a blessing in disguise.
Let us first note that there is no preferred choice for the coordinate system (24) we define for the conjugacy classes. Once a certain choice {x h i } has been made such that will do just as well. In the trivial vacuum, the S-matrix is invariant with respect to this transformation. This is most easily seen by looking at the algebraic expression (27), but it is also confirmed by our measurements of (34). This invariance can be understood on the operator level by multiplying out the representation matrices of the centralizer in equation (33). Generally, this will lead to terms of the form where g is the product of links on the loop and h k = g h i g −1 , implying that indeedg ∈ A N . When the loop is linked with another loop, the element g will in general be in the conjugacy class of the flux of this other loop. Under the transformation (35) of the conjugacy class coordinate system, the above expression will transform as h k = e. In a nontrivial ground state where noncommuting fluxes may have nonzero S-matrix elements due to a vacuum interchange contribution, the transformation (35) may manifest itself in different measured matrix elements. This means that in such cases the entry (A, α) may split into multiple entries {(A i , α i )}. As we are interested in these multiple entries, we will in our calculations always include the nontrivial behaviour of our observables under this auxiliary A N action. This turns out to be one of the two mechanisms responsible for the splitting of irreps of A into multiple irreps of U, the other of which we turn to now.
An auxiliary H/ A N symmetry. There is another symmetry but now on the level of the fusion algebra that turns out to be useful. Suppose in the theory there exists a rule of the form where (e, β) is some 1D purely electric representation. This turns out to be the case whenever the representation (e,β) (·) evaluates to unity for all elements in A N , the normalizer of conjugacy class A. We can prove this using the explicit expression for the fusion coefficients in terms of the quantum double characters: Picking (B, β) = (e, β) and (C, γ ) = (A, α), · · · Tr (A,α) (P h g) * where in the latter line we have made use of the orthogonality of the characters. We assumed that (e,β) (P e g) = 1 for all g ∈ A N . The sum over h is restricted since if h ∈ A the matrix element will be zero and the sum over g is restricted since if g ∈ A N the matrix element will be off-diagonal and thus not contribute to the trace. So we see that the fusion rule (36) leads to a degeneracy in the calculation of S-matrix elements, since by definition However, on the operator level this equality does not hold. Indeed, when we probe the lhs of this equation in a nontrivial vacuum, the result will in general differ from the rhs. In particular, it turns out that the different U representations that lift to the same A representations (A, α) differ Table 1. Character table of the groupD 2 and of Z 4 as a centralizer of the conjugacy class X i .
precisely by such a fusion. So this degeneracy may be lifted in the broken phase and give rise to additional splittings of certain entries (A, α). Consequently, in our numerical calculations we have to explicitly keep track of the presence of such electric representations (e, β), which satisfy (36), and see whether they give rise to additional splittings.
To conclude this section, we remark that we have very explicitly indicated how one gets from the modular S-matrix S ab to the extended or broken S-matrixS a i b j , from which the topological data of the broken U phase can be immediately read off.

The D(D 2 ) gauge theory
We turn to the particular example we have chosen to work out in detail: a DGT with gauge groupD 2 , also called the quaternion group. The representation theory was worked out in [36]; here we summarize the results that are required to describe the breaking by a Bose condensate.

Algebraic analysis
The groupD 2 contains eight elements that can be represented by the set of 2 × 2-matrices where the σ i , i = 1, 2, 3, are the Pauli spin matrices. We denote the conjugacy classes as e = {1}, e = {−1}, X 1 = {iσ 1 , −iσ 1 }, X 2 = {iσ 2 , −iσ 2 }, X 3 = {iσ 3 , −iσ 3 } and the irreducible group representations as 1, the trivial irrep, J 1 , J 2 , J 3 three 1D irreps and χ the 2D irrep given by (38). The character table is given on the left-hand side of table 1. The centralizer groups for the conjugacy classes e and e are bothD 2 since the elements in these conjugacy classes constitute the centre of the group. The conjugacy classes X i , i = 1, 2, 3, have nontrivial Z 4 centralizer subgroups, of which the character table is given on the rhs of table 1. The irreducible representations of the quantum double are labelled by a combination (A, α) of a conjugacy class A and a centralizer irrep α. The full set of fusion rules for the D(D 2 ) theory is given in appendix. All in all, there are 22 sectors: the trivial flux paired with the five irreps ofD 2 , the e flux paired with the five irreps ofD 2 and the three X i fluxes paired with the four Z 4 irreps. The sectors that involve an X i flux or a χ irrep have quantum dimension 2; the others have unit quantum dimension. One finds that the total quantum dimension for the theory D A = 8.
Breaking: (e, 1) condensate. In this case, the lift of the new vacuum is φ = (e, 1) + (e, 1), which implies that q = d (e,1) + d (ē,1) = 2. To determine the effective low-energy theory we fuse φ with all particle sectors of the theory and look for the irreducible combinations that appear. As before, the notation (A, α) stands for a particle with magnetic flux A and electric charge α.
The lines marked with ( * ) have components on the right-hand side that carry different spin factors, implying that they are confinement in the broken phase. Studying the fusion rules of the surviving combinations of irreps leads to the conclusion that the effective U theory is D(Z 2 ⊗ Z 2 ). We denote the four different irreps and conjugacy classes of the group Z 2 ⊗ Z 2 by the labels ++, +−, −+, −−, the first (second) symbol standing for the first (second) Z 2 . This means that D 2 T = 32 and D 2 U = 16. The branchings of A irreps into the unconfined U theory are and have d t = 2.
Breaking: (X 1 , 0 ) condensate There is an obvious symmetry in the fusion rules between the three (X i , 0 ) particle sectors. We choose to study the case when the (X 1 , 0 ) condenses. This gives for the new vacuum φ = (e, 1) + (e, 1) + (X 1 , 0 ), from which it follows that q = 4 in this case. We now read off the lifts of the T fields on the right: We have used the symbol δ i j that is 1 when i and j are equal and is zero otherwise and η i j that is 1 when i and j are not equal and is zero when i and j are equal. The U theory is D(Z 2 ) Z 2 ⊗ Z 2 . This means that D 2 T = 16 and D 2 U = 4. The lifts of the unconfined fields are (e, 1) + (e, 1) + (X 1 , 0 ) 1 → (+, +), d (+,+) = 1, and of the confined fields (e, J 2 ) + (e, J 2 ) + (X 1 , 2 ) 1 → t 1 , d t 1 = 1,

Measurements by lattice Monte Carlo simulations
The five couplings {β A } for conjugacy class A that appear in the action of the D(D 2 ) theory are inversely proportional to the masses of the fluxes A. For example if we put all couplings to zero except for β e , which we make large (at least as large as 2.0 as we will see shortly), the trivial vacuum is realized: this is the configuration where for all plaquettes U p = e. Deviations from this configuration occur because of quantum fluctuations, but since all excitations are gapped they will be exponentially suppressed. The gap in this vacuum is easily calculated to be of the order of 4β e , since the smallest excitation above the configuration in which all plaquettes are e is the one in which one link has a value h = e. This excites four plaquettes and changes the action (39) by a value of 4β e .

MC considerations.
For the other, nontrivial phases in this theory, the dominant configurations contributing to the path integral are not so readily identified. To gain insight into what configurations contribute we use an MC simulation, in particular a modified heat bath algorithm. Bluntly applying this algorithm to our problem leads to various complications; therefore we briefly point out the method, the complications and how we have resolved them. The procedure starts with some initial configuration of link variables {U } 1 . We then update all links in lexicographic order, a process called a sweep, and arrive at a new configuration {U } 2 . The updating process for each link proceeds as follows. Consider the link U i j . We identify which plaquettes contain this link: in three dimensions, there are four such plaquettes. Now we calculate, for each element g ∈ H , what the sum of the plaquette actions for each of these four plaquettes would be if U i j were to have the value g. This gives a set of numbers {S g 1 , S g 2 , . . . , S g |H | }, where S g k is the sum of the four plaquette actions with U i j equal to g k . We now calculate a localized partition sum Z U i j : which can be used to calculate a set of probabilities { p(g)} g∈H for each group element g However, for our purposes this scheme is troublesome for two reasons: it is tacitly assumed that the presence of the operator O in (40) does not change the value of the minimum of the action S and furthermore the loops of magnetic flux are very nonlocal objects and therefore highly unlikely to appear when using a local updating algorithm. This is illustrated in figure 6. The shift upward of the functional S[U ] is due to the presence of a magnetic flux string and the shift to the left is due to the nonlocality of the magnetic excitations. The latter shift also occurs when a single loop of flux is inserted.
The minimum of the action in the calculation of an S-matrix element (34) is altered by the insertions of the loop operators: the configuration for two noncommuting fluxes carries a string (see the discussion around figure 5) that is massive and thus costs a finite amount of action. There is no way to get rid of this string and therefore the minimum value of the action in the presence of the two loops is shifted. We therefore have to amend the standard MC algorithm. Defining S = S min + δS, without operator insertion, S =S min + δS, with operator insertion, and noting that around the minimum the actions behave identically, implying that δS and δS are the same functions, expression (40) becomes This leads to a modified MC average We now describe two approaches to the second problem in our MC measurements: the low probability that the local updating algorithm will converge to a gauge field configuration containing a (set of) magnetic flux loop(s). We will assume that a single loop of pure magnetic flux is inserted, as nothing substantial will change in the case of multiple loops or the addition of dyonic charge.
The first approach is based on the observation, illustrated in figure 4, that we know the gauge field configuration (up to gauge transformations) that extremizes the action in the trivial vacuum with the insertion of a loop of magnetic flux: the h-forest. We can therefore use this configuration as an ansatz in the MC algorithm. We start with a 'cold lattice', all links U i j = e, except for the h-forest-for these links we set U i j = h. This is an extremum of the action for the action if we set all β A =e = 0 and β e 1. To carry out a measurement at some other value of the coupling constants, we can slowly change the coupling constants towards the desired values, performing a few MC updates after each step.
The second approach is a more physical one. We initialize the lattice directly at the desired point in the coupling constant space. The trick then is not to insert the loop all at once, but to slowly grow it, as illustrated in figure 7. We start by twisting the action for four plaquettes around one link, as shown by the shaded plaquettes in the top left of figure 7. After this, a number of MC updates are performed. Then the set of plaquettes that have a twisted action is changed as in the top right corner of the figure. Again a number of MC updates are performed and so on. We have checked that in the trivial vacuum one obtains the h-forest configuration using this procedure. Both of the methods to insert flux loops have been used by us and we have verified that they lead to completely equivalent results.

Results
In this subsection we present the results of our MC simulations. The first quantity we measured was the free energy as a means of mapping out a suitable subspace of the parameter space. It gives us an indication of the validity of our naive intuition about where nontrivial condensates should occur.
Once we have found some region where symmetry breaking occurs, we measure the open string expectation values to determine the respective condensates. After that we measure the unbroken and broken S-matrix elements. Using the straightforward algorithm involving the auxiliary symmetries of our loop operators discussed in section 2.3 allows us to find the branching matrix n a 1 u as well as the S-matrix of the effective U theory in the broken phase. Mapping out the phase diagram. The space of coupling constants in our theory is five dimensional but it is not our goal to analyse it completely. We have restricted our search to some representative regions where nontrivial condensates do indeed occur. To study the location of the corresponding phase transitions we measured the free energy F, which we define as the expectation value of the plaquette action, averaged over the spacetime lattice.
The left plot of figure 8 shows F as a function of (β e and βē) and all other couplings are equal to zero. For small values of all the couplings appearing in the action (39), we are in the completely confining phase of the gauge theory, where all the open string operators of magnetic flux have a nonzero expectation value, and all loop operators carrying electric charge are confined. This corresponds to the plateau in the graph where F is maximal and tends to zero.
The regions where the magnetic flux (e, 1) and (X 1 , 0 ) have condensed can be anticipated on theoretical grounds by realizing that the coupling β A is inversely proportional to the mass of flux A. In fact, when we look at the subgroup K A generated by the elements in conjugacy class A, in particular K e = {1, −1}, and set the couplings for the conjugacy classes containing the elements in K A equal to one another, there is an extra gauge invariance U p → k U p for an element k ∈ A in the plaquette action (39). In particular is invariant with respect to U p → −1U p and These left multiplications are exactly the kind appearing in the definition of the (loop) order parameters (33). Therefore, one can establish, even without reverting to MC measurements, that the above actions, for large values of β, produce the desired flux condensates.
One may verify this reasoning in figures 9 and 10 where we have probed the phase diagram in more detail by measuring the spacetime averaged expectation value of δ A (U p ) for all conjugacy classes A as a function of the relevant coupling parameters β. The red colour indicates high values for the expectation value and we see that for all coupling parameters near zero all fluxes are condensed and thus all charges will be confined. This is what is traditionally called the 'strong coupling' phase (g ∼ 1/β 1). Looking at the colourings for the various operators, one readily identifies the various phases as indicated in the schematics of subfigures (f). For example the symmetry with respect to the diagonal of figures 9(a) and (b) shows that there are 'Ising'-like ordered phases, one with all plaquette values U p = e and the other with all U p =ē. The in between region is the region with theē flux condensate. Note that if theē flux would be the only one that phase would continue all the way to the origin, and we would exactly end up with the Z 2 pure gauge theory phase diagram.
In the region with β e larger than approximately 2.0 and all other couplings near zero, the trivial vacuum is realized. All string operators with nontrivial magnetic flux have zero expectation value there.
As pointed out in previous sections for example in relation (12), there is a very direct way to determine the condensate as well as the quantum embedding index q. This is by measuring the expectation value of the open string for each pure flux A and then summing over all fluxes. In the (e, 1) vacuum we obtain There is one more issue we like to address in our simulations, i.e. to determine the order of the transitions we have identified. A conventional approach is to search for a hysteresis effect across a first-order transition, but because of the relatively modest size of the lattices used this is not an optimal approach. A method that is working much better is to directly probe the system at a given sequence of coupling constants around the transition and to see whether there is a coexistence region where both phases occur in the sampling 8 . To make these measurements we use the parallel tempering method [39] to overcome local minima in the action landscape. The idea behind this method is to initialize a range of lattices simultaneously, all at different couplings along a trajectory in coupling constant space starting in phase 1 and ending in phase 2. The updates of this ensemble then consist of the updates of each of the individual lattices and, occasionally, a swap of two adjacent lattices. The swap between lattices 1 and 2 is accepted with a probability where S 1 (2) means using the action (in particular, the set of couplings) of lattice 1 to evaluate the field configuration of lattice 2 and so on. One can prove that this satisfies detailed balance. In effect, each lattice will perform a random walk through coupling constant space along the chosen trajectory, allowing a 'cold' lattice to thermalize in the 'high-temperature' region, thus overcoming the local minima of the action. We have made measurements for the trajectories indicated by the arrows in figures 9(f) and 10(f). The results of these measurements for the horizontal arrow are given in figure 12 and those for the vertical arrow in figure 13. We find that in the horizontal trajectory the transition from the strongly coupled phase corresponding to the left peak in figure 12 to the trivial phase The result for the vertical trajectory corresponding to the transition from the trivial phase to the broken (X 1 , 0 ) phase is given in figure 13, where we see that the peak shifts continuously, implying that the transition is second order. We can understand this transition as follows. In this region of coupling constant space, all fluxes except the e flux are very heavy. This means that the ground state is essentially that of a Z 2 gauge theory. Since Z 2 gauge theory in three dimensions is the Kramers-Wannier dual to the 3D Ising model, it has the same phase structure [34]. We therefore expect this phase transition to lie in the same universality class.
Measuring the (broken) modular S-matrices. We have measured the (broken) S-matrix elements using the simple algorithm involving the auxiliary symmetries of our loop operators. This allows us to obtain the unbroken S-matrix as well as the branching matrix n a 1 u and the S-matrix of the effective U theory in the various broken phases. Here we exploit the relation (27) for the measurement, and the relation (18): relating S uv to the measured S-matrix in the broken phase. We first measured the unbroken S-matrix in the D(D 2 ) phase and obtain results identical to the matrix calculated using the defining formula (6); the result is given in table 2 and is of course consistent with the matrix obtained from the relation (18) with = 0. The accuracy of the measured matrix elements represented in the table as integers is smaller than 5%. The branching matrices n a u can be obtained from measuring the broken S-matrices. The columns in these matrices correspond to the different U sectors. If we see two rows or columns with different parents a, b in the A theory that are proportional to each other, a and b branch to the same U sector u. Conversely, if different u fields correspond to the same a field that means that the a splits in the broken phase. We have listed the results for the broken S-matrix in the (e, 1) vacuum in table 3. To realize the splittings between the irreducible representations using the auxiliary gauge symmetries alluded to in section 2.3, we found the following construction to suffice.
We see that in table 3 the columns (rows) for the sectors (e, α) and (e, α) for α = 1, J 1 , J 2 , J 3 are identical and thus that the corresponding fields have to be identified. This leaves us with 16 sectors for the broken U theory. Summing the entries as prescribed by formula (3.3) yields exactly the S-matrix of the D(Z 2 ⊗ Z 2 ) theory, which is given in table 4. The S-matrix for the (unbroken) D(D 2 ) theory (up to the normalization factor 1/D A = 1/8) as measured in the trivial vacuum. We put integers in the table as the accuracy is below 5%, i.e. one should read 1 as 1.0 ± 0.05.

Conclusions and outlook
In this paper, we have studied euclidean lattice models for discrete gauge theories. We have introduced a set of multiparameter actions for these theories that display a rich phase structure, and showed in particular that all the allowed condensates of pure magnetic flux are realized in certain well-anticipated regions of coupling constant space. The set of open string operators that we defined form a set of order parameters that allowed us to determine the content of the condensate and to measure the topological symmetry breaking index q.
Once the condensate is identified, we have shown how to unambiguously reconstruct the S-matrix of the low-energy theory in a broken or unbroken phase by measurements of the complete set of braided loop operators, using the anyonic loop operators we proposed in the earlier work [16]. Due to an auxiliary gauge symmetry these operators are particularly well suited for detecting the nontrivial splittings of fields that correspond to fixed points under fusion with the condensate. We found that as expected the excitations that are confined in a broken vacuum give rise to rows and columns of zeros in the broken S-matrix. Our work clearly demonstrates that the euclidean approach allows for a very straightforward method to completely determine the nature of the broken phase.
Our work showed that the reason why the modular S-matrix changes in the broken phase is largely due to the contribution of the so-called vacuum exchange diagram. In an upcoming more theoretical paper [29], we will extend the approach used in this work, the use of observables and in particular the S-matrix to determine the phase structure of a TQFT to a far wider range of theories, in particular the SU(N ) k TQFT arising from the Chern-Simons actions.
It would be interesting to study different models exhibiting different topological phases by somehow formulating them in the euclidean 3D framework; to our knowledge such an approach is unfortunately not yet available for the Chern-Simons theories. One expects that for the Levin-Wen models [40] our approach could be implemented however. Another path is to investigate the phase structure after adding dynamical matter fields that transform nontrivially. It is known that in such situations the Wilson-type criteria break down as the strings can break; this necessitates the development of different diagnostic tools [35,41].