Generalized modular transformations in 3+1D topologically ordered phases and triple linking invariant of loop braiding

In topologically ordered quantum states of matter in 2+1D (space-time dimensions), the braiding statistics of anyonic quasiparticle excitations is a fundamental characterizing property which is directly related to global transformations of the ground-state wavefunctions on a torus (the modular transformations). On the other hand, there are theoretical descriptions of various topologically ordered states in 3+1D, which exhibit both point-like and loop-like excitations, but systematic understanding of the fundamental physical distinctions between phases, and how these distinctions are connected to quantum statistics of excitations, is still lacking. One main result of this work is that the three-dimensional generalization of modular transformations, when applied to topologically ordered ground states, is directly related to a certain braiding process of loop-like excitations. This specific braiding surprisingly involves three loops simultaneously, and can distinguish different topologically ordered states. Our second main result is the identification of the three-loop braiding as a process in which the worldsheets of the three loops have a non-trivial triple linking number, which is a topological invariant characterizing closed two-dimensional surfaces in four dimensions. In this work we consider realizations of topological order in 3+1D using cohomological gauge theory in which the loops have Abelian statistics, and explicitly demonstrate our results on examples with $Z_2\times Z_2$ topological order.

Topologically ordered quantum phases of matter in 2+1D have been intriguing since their discovery decades ago (see 1 and references therein), due to exotic properties such as fractionalized quasiparticles with anyonic quantum braiding statistics. 2,3 Early on it was realized that in such phases the topological degeneracy of the ground state on the torus corresponds to the number of types of particle excitations (superselection sectors). 4 Furthermore, it was shown that the matrix of Berry's phases experienced by the ground states under the modular transformations of the torus, the S and T transformations (Fig. 1a), are directly related to the quantum statistics of the quasiparticles. 2 In fact, to date the most fundamental conjecture remains that the matrices of S, T contain complete information about a topological order. 2 Therefore one can view the modular S, T matrices as the "non-local order parameters" in a topologically ordered phase. 5 However, in three spatial dimensions some fundamental questions are yet completely unresolved: Is there a physical way to characterize different topological orders in 3+1D? Can braiding of excitations help us in the characterization? Clearly the problem is much more complex, since in 3d there are generically both point-like and looplike excitations, and their geometric interplay is rich. If some type of braiding can help us characterize the topological order in 3+1D, what is the topological property of that braiding process that is relevant?
Motivated by the fundamental role of modular transformations of the torus in 2+1D systems, our approach to these questions is based on considering the analogous transformations on the three-torus (e.g., a cube with periodic boundary conditions). The modular transfor- mations S, T on the torus generate the group SL(2, Z), which represents the different classes of continuous transformations on the torus. 6 In 3+1D quantum states, the analogue is the three-torus, which also has just two associated transformations S, T , generators of SL(3, Z) group, 7,8 namely a 120 • rotation through a diagonal of the periodic cube and a shear, respectively (Fig. 1b). Very recently it has been conjectured that exactly these kinds of transformations can be used to characterize topological order in any dimension. 8 One way to study topologically ordered states is using the exactly solvable models of discrete gauge theories introduced by Dijkgraaf and Witten (DW). 9,10 Although these theories in 2+1D do not provide an exhaustive classification of all possible topological orders, 11 they describe a physically interesting set of states. Most importantly for this work, such cohomological gauge theories with gauge group G are naturally defined in any spatial dimension, allowing us to study 3+1D topological orders. They also host both point-like and loop-like excitations, namely gauge charges and flux-loops, respectively. For simplicity, we restrict to the case of Abelian groups G, and then additionally to cases where loops have only Abelian braiding.
Our approach to the 3+1D problem is based on generalizing some aspects of the 2+1D case, which we review here. It is well known that topological operators, which describe tunneling of quasiparticle across the periodic 2d system, obey a non-trivial algebra. A certain product of these operators gives the identity operator times a complex number which equals an S matrix element, and also equals the quasiparticle statistics. 3, 12 To generalize to three dimensions, we need the topological content of this relation, which is revealed using a picture. Fig. 2c depicts the expectation value of the relevant product of tunneling operators as a time sequence of events where particleantiparticle pairs tunnel across the periodic system. Note that the initial and final state in the picture are the same. The four worldline segments in this process, belonging to two quasiparticle types, can be connected to reveal two linked worldlines, Fig. 2d. The picture shows how the S a) b) c) d) FIG.
2. An S matrix element and braiding in 2+1D. (a) Matrix element equals an overlap of two MES, shown as a time sequence. Any MES is defined by action of particle tunneling operator along x on the appropriate reference state defined by x direction. The MES at later time (blue) has been acted on by S modular transformation, so the tunneling is along y, and the final reference state is defined by y direction. (b) Embedding the space-time process of (a) in three dimensions shows that the two worldlines are linked. Time grows in the radial direction as shown, so the space-time of (a) spans the volume of a toroidal slab. (c) Alternatively, the matrix element equals a product of tunneling operators, also presented as a time sequence. Arrows mark the action of tunneling operator and its inverse. Due to taking the expectation value of this operator product, the initial and final state are the same, in contrast to panel (a). (d) Connecting the worldlines from (c) one obtains two worldlines that are linked.
matrix element and the braiding statistics can be seen in the linking of worldlines. One can detect the linking of worldlines in a purely algebraic way, without drawing the two figures, by calculating the linking topological invariant of the two worldlines. This invariant simply counts the number of links in the one-dimensional worldlines living in the three-dimensional spacetime. There is another way to relate an S matrix element to a topological property of the quasiparticle braiding process in 2+1D, which will be important for our generalization to 3+1D. Recently, the connection between modular transformations on the ground state manifold in 2+1D and the statistics of quasiparticles was further exposed by the introduction of minimum entropy states (MES), a special choice of basis in ground state manifold. 13 Namely, an S matrix element is related to an overlap of two MESs. 13 This relation is due to the fact that every MES can be created by the action of a quasiparticle tunneling operator on the appropriate reference state. 13 A pictorial representation again uncovers the topological content in this understanding: Fig. 2a depicts the overlap of two MESs as a time sequence of applying quasiparticle tunneling operators. Note that here the initial and final states in the picture are different reference states, being rotated by 90 • . In this case we immediately obtain the two worldlines, since each MES in the overlap contributes one. The picture reveals that the two worldlines in this process are linked when the 2+1D space-time pro-cess is embedded in three-dimensional space as in Fig.  2b. This approach again connects the measurement of braiding statistics (S matrix element) in 2+1D to the linking of particle worldlines.
In this paper, we will generalize both above approaches to 3+1D cohomological gauge theory, and obtain some surprising relations between the matrix elements of the three-torus S, T transformations and the braiding of excitations. Noticing that the line, which represents a quasiparticle tunneling operator across a periodic direction on the two-torus, becomes a membrane, which represents the tunneling of a flux-loop across two periodic directions in the three-torus, we will construct the appropriate membrane operators as well as the MES on the threetorus. We will show that in 3+1D topological order, the non-trivial matrix elements of the S, T transformations in the MES basis are due to a non-trivial algebra of the topological membrane operators, analogously to the 2+1D case.
We will next generalize to 3+1D the two 2+1D approaches described above. Namely, we will study the space-time process representing the non-trivial product of membrane operators, i.e., the loop tunneling operators, as well as the space-time process of overlapping two MESs. Both calculations give the same value, which equals a matrix element of the three-torus S transformation. This value also represents a topological quantum phase accrued during some loop space-time process, and the main question becomes: What is the nature of this loop process? In the analogous situation in 2+1D, we described that the particle process is simple braiding, which is characterized by the linking of worldlines.
Most strikingly, we will argue that the obtained S matrix elements in 3+1D relate to certain braiding processes involving three loops simultaneously. This is surprising since there is a simple, seemingly fundamental, braiding process of just two loops, where one loop traces out a torus enclosing the other loop, which is relevant in other physical contexts. 14,15 To uncover the topological underpinning of the two studied 3+1D space-time processes, we will consider the topological invariants of the loop worldsheets in both of them. Although we cannot draw the pictures of twodimensional worldsheets living in the four-dimensional space-time, in analogy with Fig. 2, we will find through calculation that the loop worldsheets in both space-time processes are described by the exact same non-trivial values of the triple linking number. 16 The triple linking number (TLN) is an integer invariant of closed surfaces in four dimensions, and can be seen as the 3+1D generalization of the linking number of closed lines in three dimensions which was relevant for particles in the 2+1D case.
The TLN is obviously the fundamental, truly threedimensional descriptor of the relevant three-loop braiding process which is revealed in the context of this paper. To capture its meaning in a tangible way, in this paper we will also construct a movie of a process involving three loops, shown in Fig. 3, such that the three loop world- sheets in space-time exhibit exactly the same values of TLN as obtained in the two processes above. This movie shows the physical braiding of three flux-loops resulting in a topological quantum phase which measures the properties of the underlying 3+1D topological order (at least for cohomological gauge theory), and equals a matrix element of the three-dimensional S modular transformation. This paper is organized as follows. In Section II, we define the exactly solvable models in 3+1D, which are classified by cohomology group and can be viewed as extension of Dijkgraaf-Witten theory to 3+1D. In the following section, we put these models on three-torus, and find the ground state manifolds. Particularly, we find a MES basis, which is useful for interpretation. Further, we construct membrane operators defined as operators mapping between MES. We work out the modular transformations on MES basis in Section IV. We find modular transformations to be directly related to braiding statistics of flux-loops and particles. We show this by both geometric and algebraic methods. In the last Section we solve these models for some illuminating examples.

II. COHOMOLOGICAL GAUGE THEORY IN 3+1D
In this section, we define the cohomological gauge theory for a general manifold in 3+1D, based on the Dijkgraaf-Witten (DW) topological invariant. The theory is topological and defined by a discrete gauge group FIG. 4. The 4-cocycle ω assigns a U (1) complex number ω ε (g54, g43, g32, g21) to a 4-simplex, where ε is the chirality of the 4-simplex, defined as ε = sgn[det ( 12,23,34,45)]. The dashed lines represent that the vertex 5 has a different coordinate in the fourth dimension (time) with respect to the other vertices.
G. However, there are distinct topologically ordered states for a fixed G, and in 3+1D they are classified by the fourth cohomology group of G with coefficients in U (1), namely H 4 (G, U (1)). In Appendix A we give a brief review of cohomology concepts relevant for the rest of the paper, while referring the reader to, e.g., Refs.17 and 18, for more details.
In this paper we will work in 3+1D, and therefore the theory will be defined using the 4-cocycle (sometimes we call it simply cocycle) ω, for which the cocycle condition becomes: ω(g 2 , g 3 , g 4 , g 5 ) · ω(g 1 , g 2 · g 3 , g 4 , g 5 ) · ω(g 1 , g 2 , g 3 , g 4 · g 5 ) (1) =ω(g 1 · g 2 , g 3 , g 4 , g 5 ) · ω(g 1 , g 2 , g 3 · g 4 , g 5 ) · ω(g 1 , g 2 , g 3 , g 4 ), where ω ∈ H 4 (G, U (1)), and g i ∈ G. In this paper, we will use the "canonical" 4-cocycle, meaning that ω(g 1 , g 2 , g 3 , g 4 ) = 1 if any of g 1 , g 2 , g 3 , g 4 is equal to 1 (the identity element of group G). The gauge theory is now defined by using ω to construct topological invariants of a 4D manifold. For a given 4D manifold M without boundary, one can triangulate it using a finite number of 4-simplices. 19,20 The 4-simplex is a higher dimensional analogue of a regular polyhedron, and can be constructed from a tetrahedron by adding a fifth vertex and moving it into the fourth dimension so that all the edges from it to the four original vertices are the same length as the tetrahedron's edges. The triangulation basically requires completely filling the manifold M with 4-simplices without overlaps. The vertices of this triangulation are then ordered arbitrarily, and the ordering is represented by assigning arrows going from the lower to the higher ordered vertex on each edge, Fig. 4. Let us denote a 4-simplex of the triangulation, together with the ordering of its vertices, by σ I , where I = 1, 2, . . . , S labels 4-simplices and S is the total number of 4-simplices in M . Next, one defines a coloring ϕ of all the edges in the triangulation, by assigning group element to them. Let us denote the group element assigned to the bond connecting vertices j and i as g ij , following the ordering from j to i: j → i; we then automatically assign g ji = g −1 ij . In addition, the three assigned group elements for any given face must satisfy the constraint g ij · g jk · g ki = 1, and i, j, k are the three vertices of the face. This constraint is the "zero-flux rule".
With these definitions, one can assign a U (1) phase to every 4-simplex by computing ω ε (g 54 , g 43 , g 32 , g 21 ), where ε = sgn[det ( 12,23,34,45))] determines the chirality of the simplex, as shown in Fig. 4. 21 For a given coloring ϕ and simplex σ I , we label this U (1) phase as W (σ I , ϕ) ε(σ I ) . Finally, one can compute the product of all W for the simplices: For a given coloring ϕ, we will have one such product. The key result 9 is that the complex number where |G| is the number of elements in group G, and V is the number of vertices in the triangulation, is a topological invariant of the manifold M . More precisely, Z M does not depend on the triangulation and the ordering of vertices (while different colorings are already summed over), owing to the cocycle condition in Eq. (1). One can further show that equivalent cocycles (i.e., cocycles differing by a coboundary) give the same value of Z M . 9 The topological invariant Z M is exactly the partition function of the cohomological gauge theory, which is a topological quantum field theory for discrete gauge group G in 3+1D. It is the higher dimensional version of the DW theory 9,22 , and it only depends on inequivalent elements in H 4 (G, U (1)).

A. Exactly solvable models
We define our exactly solvable models in 3+1D as Hamiltonian versions of the cohomological gauge theory. We consider space triangulated using a tetrahedron lattice with oriented edges (bonds), where these orientations are compatible with some ordering of lattice sites, and assign an element g ij ∈ G to each oriented edge j → i, according to the above discussion.
An arbitrary quantum state in the Hilbert space H of our model is then labeled by |a = |{g ij } . The building block for the Hamiltonian is the operatorB s p labeled by a group element s ∈ G, and a "plaquette" p containing all 3-simplices (tetrahedra) that share the vertex i. The plaquette operator acts on group elements on the edges that share i. To define its action, we introduce an additional edge rising into the fourth dimension, connecting i to an auxiliary vertex i . To edge i → i we assign the element s ∈ G. The group elements are changed as and these new values are represented on auxiliary edges i → j and k → i . Further, the non-zero matrix elements ofB s p , namely B s p = f(s)|B s p |i , are assigned the following quantum amplitude where the 4-simplices σ I are built by triangulating the 4D volume formed by the tetrahedra in the plaquette p and the auxiliary edges.
It is important to note that the zero-flux rule is by construction satisfied on all faces (triangles) of 4-simplices, if it is satisfied in the tetrahedra of p, and this must be imposed for the B s p to be well-defined. We can then define the plaquette operatorsB p as having matrix elements TheB p are projectors, which can be easily checked using the cocycle property to show f|B s pB s p |i = B s·s p , which then implies f|B pBp |i = B p . Similarly, it can be shown that the plaquette operators commute, [B p , B p ] = 0, ∀p, p .
Let us also introduce the operator Q t , which projects flux in a triangle t to zero, i.e., it enforces the zero-flux rule. Then the Hamiltonian takes the form where the label t ∈ p enumerates all the triangles making up the plaquette p. As mentioned above, the factor t∈p Q t is actually crucial to ensure that H is welldefined. Further, it is easy to see that plaquette operator termB p t∈p Q t actually commutes with the projectors Q t . Since all the terms in H commute with each other, the model is exactly solvable. Let us briefly mention the connection of the Hamiltonian formulation to the gauge theory, which is exhibited in the ground state manifold. Since all the terms in H are projectors, the ground state manifold is the image of the projector P = pB p t∈p Q t . On the other hand, P is exactly the projector defining the cohomological gauge theory on the 4D manifold having two copies of our spatial manifold M as boundaries (see Ref. 18 for details). The ground state sector of H, to which P projects with eigenvalue 1, is also the ground state sector of the cohomological gauge theory 9 defined on M .
. This phase can also be written as

B. Geometrical reduction of 4-cocycles
In this section we present some cohomology equations for reducing the 4-cocycle to lower order cocycles, and explain their geometric meaning. These equations crucially simplify all following calculations. From now on, we will focus on Abelian groups G for convenience.
First, let us consider a triangulated 4D manifold in Fig.5, with the shown coloring. (Note that some edges needed for full 4D triangulation are omitted, but coloring and ordering are fully defined.) The U (1) phase calculated from all the 4-simplices spanning this 4D volume, with the 4-cocycle ω given, equals β ε s (c, b, a), with: and ε = sgn[det( a, b, c, s)]. Using the 4-cocycle condition for ω, it is straightforward to show that β s is a 3-cocycle. This shows that lifting all vertices of a tetrahedron produces a quantum phase which is only a 3-cocycle, for any given ω.
Another quantity that appears naturally from a cubic geometry is γ a,b , whose geometric meaning is shown in Fig. 6. It is defined from the 3-cocycle β a as: It is straightforward to show that δγ a,b (c, d, e) = 1, namely, γ a,b is a 2-cocycle (see Appendix A). Further, from Eq. (7) and Eq. (8), . This equality follows also from the geometry in Fig. 6.

III. GROUND STATE ON THREE-TORUS AND MEMBRANE OPERATORS
A. Exact models on three-torus We now put our model, Eq. (6), on the three-torus in 3+1D. It is important to note that the exactly solvable model has correlation length zero. Therefore, we can consider the simplest triangulation of a three-torus shown in Fig.7. All eight cube vertices are identical due to periodic boundary conditions. It is triangulated by six tetrahedrons. There are three independent edges, which are assigned group elements a, b, c ∈ G, with G a finite group. Edges with the same direction share the same group element value. The corresponding quantum state is labeled by |a, b, c . We also require G to be Abelian for simplicity.
Since there is only one vertex, we denote the plaqutte operatorB p simply asB, which equals 1 |G| s∈GB s . The action ofB s on state |a, b, c iŝ We can directly write down the above result due to the observation that the 4D graph we obtain by acting witĥ B s is in fact made out of two copies of Fig. 6. 23 Notice that the U (1) phase obtained by action ofB s is a fully antisymmetric function of a, b, c, s, as can be seen both geometrically and algebraically.

B. MES as ground state basis
Let us first briefly review topological order in 2+1D. It is partially characterized by ground state degeneracy on torus. 3 One can understand this degeneracy by applying Wilson loop operators of distinct topological excitations winding around one of noncontractible loops on the torus. From this point of view, one can see that ground state degeneracy equals the number of distinct topological superselection sectors.
Nonchiral topological order is fully determined by braiding statistics and topological spin of its topological excitations. 4 Remarkably, one can read the information about excitations from ground state by using modular transformations 2,13 , namely, by considering the S, T matrices of the S, T transformation in the ground state manifold. Dimension of S, T equals the number of topological sectors. In a proper ground state basis, we can obtain the "canonical form" of S, T matrices, for which the entries of S matrix are the braiding statistics and the diagonal elements of T are the topological spins of quasiparticles. Ground state basis for canonical S, T matrices is formed by minimal entropy state (MES). 13 We can extend these concepts to 3+1D. However, there is a major difference in this case: Topological excitations can be flux loops in 3+1D. Without loss of generality, we only consider the MES in z direction.
Inspired by the case of 2+1D cohomological gauge theories discussed in Ref. 24, we have found the MES in z direction as whereχ a,b λ is a one-dimensional projective representation. Here, λ labels different projective representations of the group G (see Appendix A), and the 2-cocycle γ from Eq. (8) plays the role of factor-system of these projective representations: We will only consider the case of Abelian (onedimensional) projective representations χ a,b in this paper. This assumption implies that the 2-cocycle γ a,b is a 2-coboundary. We believe this is related to the physical assumption of Abelian statistics of loops. Note that our assumption about γ a,b implies that γ a,b (c, s) = γ a,b (s, c), so that the phase factor in Eq. (9) is just identity, and states |a, b, c for any a, b, c ∈ G are in the ground state manifold. Firstly, we verify that this state is indeed in the ground state manifold. Acting with projection operatorB on the state, we get where the second row uses Eq. (9), and in the third row we used γ a,b (c, s) = γ a,b (s, c) which follows from the above mentioned assumptions. Next, we prove that this state is indeed an MES in z direction. Let us retriangulate the three-torus, so that it has two unit-cells in z direction. The ground state defined on this two unit-cell system can be evolved from that in one unit-cell, as shown in Fig. 8: As seen from the above, |a, b, λ defined on two unit-cells can be written as a direct product state. So, entanglement entropy of this state in z direction is zero, which must be minimum. We therefore conclude that this state is indeed an MES in z direction. Note that the entanglement entropy for this small system size vanishes due to a cancellation between the "area law" contribution and the topological contribution (see Ref. 13 and references therein), which will not occur for larger system sizes.
Similarly, it is easy to write down the MES in x and y directions: |a, ν, c = 1 whose properties can be derived in the same way as above.

C. Membrane operator
Although we constructed the MES in 3+1D, the physical picture is still unclear. Recall that in 2+1D all MES can be obtained from inserting ribbon operators (Wilson loop operators) into "trivial" MES, which corresponds to topological trivial sector. In the following we will show that membrane operators are the relevant operators for such a procedure in 3+1D.
Let us start with the MES in z direction, |a, b, λ . Characteristically in discrete gauge theory, we can interpret a group element as a label of flux-loop. Consequently, a membrane, which is the 3+1D analogue of the Dirac string, is also labeled by a group element. Further, a group representation labels a particle. 10 Then |a, b, λ can be viewed as state with membrane a in yz plane and membrane b in zx plane, as well as string λ (world-line of particle) in z direction. So, it is natural to define trivial MES as where e ∈ G is identity element. Here 1 means the trivial linear representation. The central question becomes: What are the operators that send one MES to another? It is natural to assume that these operators correspond to membrane insertion in yz and zx plane, as well as string insertion in z direction. Besides, we expect that a string in x(y) direction can measure a membrane in yz(xz) plane while membrane in xy plane will measure strings in z direction.
Following this intuition, we define membrane insertion operators in yz, zx, xy planes, respectively, as shown in Fig. 9: where u, v, w label the spatial planes of the membranes. Further, we can define where we interpret F (z) u,λ as inserting membrane u (in yz plane) and string λ in z direction, and interpret G (z) v,λ as inserting membrane v (in zx plane) and string λ in z direction. To confirm this, we act with these operators on state |e, e, 1 , getting It is not hard to obtain the "fusion rule" of membranes and strings, namely where the representation λ 3 is determined as the one in which every element g ∈ G is represented by the product of numbers which represent g in the λ 1 and λ 2 representations (note that the Abelian representations considered here are always one-dimensional, i.e., just numbers). The fusion rules follow from the properties of the 2-cocycle γ, namely, assumeχ a,b µ is a projective representation with factor system γ a,b , Then it follows that where the representations µ 1 , µ 2 , µ 3 are related in the same way as λ 1 , λ 2 , λ 3 just above. Similarly to the above derivations, we can define where H (x) (H (y) ) creates membrane in xy plane and string in x(y) direction. Acting with these operators on MES in z direction, we get It is then natural to interpret H (x) (H (y) ) as operator that measures strings in z direction and membrane in yz(zx) plane.
We will also write down the remaining two operators that send MES to MES for later convenience:

IV. TOPOLOGICAL OBSERVABLES AND THEIR PHYSICAL INTERPRETATION
A. S and T matrices from modular transformations In this section, we will calculate the Berry phase of ground states obtained during modular transformations. The derivation is largely a higher dimensional generalization of 2+1D case in Ref. 24.
In real space, we can write the modular transformations, Fig. 1b The question is what is the action of S and T on our exact models? We follow the strategy of Ref. 24, but generalize it to 3+1D. We consider a T 3 × [0, 1] manifold (T 3 is three-torus), and put the initial ground state at T 3 ×0, final state at T 3 × 1. Then we carefully triangulate the 4D manifold T 3 × [0, 1] and compute the quantum amplitude from the initial to the final state. After lengthy but straightforward calculations, we find: −1 c, a)|a, b, a −1 c . Now we act by T 31 on MES in z direction: We can see that |a, b, λ is indeed an eigenstate of T matrix. We can also get S matrix element in z-direction MES basis: Taking into account our assumption that γ a,b is a 2coboundary, the projective representationχ can be rewritten asχ ab µ (g) = ε a,b (g)·χ µ (g), where χ µ (g) is an ordinary linear representation of G, and ε a,b is a 1-cocycle for which γ a,b = δε a,b (see Appendix A). Then we get a factorized form: While the physical meaning of this element is not so clear for general case, it is instructive to see the simple case where a = b = e. Then the 1-cocycle part of Eq.(31) is trivial and only χ λ (b )/χ λ (a) is left. We can interpret this phase as Aharonov-Bohm phase of particles going around a flux loop in three dimensions, namely, particle λ sees flux loop b and particle λ sees flux loop a. In the following, we will show how the most general form of S matrix element, including the 1-cocycle contribution, can be interpreted as statistics of flux loops as well as particles.

B. Braiding statistics from S matrix and triple linking number
In 2+1D space-time dimensions, the S matrix elements can be directly related to quasiparticle braiding by considering quasiparticle tunneling operators. 3,13 One may ask is it possible to capture the S matrix elements in 3+1D space-time dimensions using a loop braiding process by considering membrane operators? We will show that this is indeed possible.
Let us review the 2+1D case, mentioned in the Introduction. An S matrix element can be expressed as an overlap of two minimum entropy states (MES), where every MES can be created by the action of a quasiparticle tunneling operator on the appropriate reference state. 13 Fig. 2a depicts the MES overlap as a time sequence where application of tunneling operator is represented as a space-time event of particle-antiparticle pair tunneling across the periodic system. We get two worldlines, since both MESs in the overlap contribute one. The occurrence of braiding can be revealed by realizing that the two worldlines in this process are linked when the 2+1D space-time process is embedded in three dimensional space as in Fig. 2b. Therefore an S matrix element describes braiding which can be seen as linking of worldlines.
We now present an analogous interpretation of an S matrix element in 3+1D, using membrane operators and a process involving a triple linking of worldsheets in 3+1D. Starting from the S matrix element in the MES basis in z direction, Eq. (30), it is straightforward to show where the last row follows from the definition of membrane operators. In contrast with the 2+1D case, each MES is created by the action of two membrane operators on the reference state |e, e, 1 . A membrane operator, such as Eq. u,λ describes the tunneling of loop in zy plane. (A membrane also contains a string in its plane, which represents particle tunneling.) The key question now becomes: What is a robust characterization of the process involving the four membrane operators in Eq. (32)? Inspired by the 2+1D case, where embedding the process involving one-dimensional worldlines in three dimensions revealed their linking, we embed the process involving two-dimensional worldsheets into four dimensions. We find that three worldsheets have a non-trivial triple linking number (TLN), which is a topological invariant generalizing the linking number used for worldlines in 2+1D. The details of embedding the 3+1D space-time process into four dimensional space are in Appendix B.
Here we briefly summarize the triple linking number (TLN) invariant.
The triple linking number T lk M N P of three oriented two-dimensional surfaces M, N, P smoothly embedded in four dimensions was defined in Ref. 16 as an analogue of the linking number of classical links. In our case, M, N, P are the three flux-loop worldsheets in 3+1D.
T lk M N P is an integer topological invariant. 25 It can be non-zero only if the surfaces M, N, P are distinct, and the T lk obey the relations T lk M N P + T lk N P M + T lk P M N = 0, (33) T lk M N P + T lk P N M = 0, and are therefore fully determined by two integers. 25 There are different ways to calculate the TLN. 25 We describe the one that is most convenient for the braiding problem: One projects the surfaces M, N, P from 3+1D onto a three-dimensional slice using an arbitrary projection direction, and looks for triple-points, namely, points in the projected manifold where all three projected surfaces intersect. For each triple point s one checks the stacking order of surfaces along the projection vector, and assigns the top surface to I s , the middle to J s , and the bottom to K s . Finally, the sign s is calculated as the handedness of the three I s , J s , K s surface normals at the point s, see Fig. 10. Now let (I, J, K) be a permutation of the three surfaces (M, N, P ). Having the above information, T lk IJK equals the sum of s over the points s for which I s = I, J s = J, K s = K. If no triple point contributes to a certain choice IJK, then T lk IJK = 0, and this has to be consistent with other values of I J K according to relations Eq. (33).
The number of triple points and the stacking order of surfaces both depend on the chosen projection vector in 3+1D, however, the resulting TLN is topologically invariant.
Here we present the TLN result, leaving the calculation details to Appendix B. We find that in the embedded process the worldsheets corresponding to membrane operators F We note that even for general 3+1D topological order, beyond cohomological models, the S-matrix elements are given by MES overlaps, while the direct connection between MES overlap and triple linking of tunneling operator worldsheets, such as discussed above and demonstrated using Eq. (32), remains a general property based on the purely geometrical embedding construction. The only required ingredient is that the MESs of the topological order can be expressed using tunneling operators of loop-like excitations, which we expect to be generally true in 3+1D.

C. Braiding statistics from membrane operator algebra
Returning to the 2+1D case, as briefly reviewed in the Introduction, there is another way to relate an S matrix element to a topological property of the quasiparticle braiding process. Namely, the quasiparticle tunneling  11. Movie for process G −1 F −1 H −1 F HG . The worldsheets in this process share the same topological properties as the three-flux-loop braiding process in Fig. 3.   FIG. 12. (color online) Projection of the membrane process movie in Fig. 11 to three-dimensional space at t = −∞. Lines show the pairwise intersections of projected worldsheets. Purple lines: for F and G worldsheets; Orange: for F and H; Green: for G and H. Although there are eight triple points here, the triple linking is still the same as for the three-fluxloop braiding process, Fig. 13. The directions t1,2,3 show the time ordering of contributions to projection from worldsheets G, H, F , so clarify to which T lkIJK some triple point contributes (see after Eq. (33)). For example, at point a, direction of t1, t2, t3 shows that worldsheet projection at this point comes from: G rather than G −1 , H −1 rather than H, F −1 rather than F , respectively. Therefore point a contributes to T lkF HG.
operators obey a non-trivial algebra, and a certain product of these operators gives the identity operator times a complex number which equals an S matrix element. 3,12 Fig. 2c depicts such a product of tunneling operators as a time sequence of events where particle-antiparticle pairs tunnel across the periodic system, in the same fashion as in Fig. 2a. The four worldlines in this process can be connected to reveal two linked worldlines, Fig. 2d. Notice that this procedure can be done in the representation of the system as a parallelepiped with periodic boundary conditions, without need for embedding the 2+1D space time in a three dimensional space. One again confirms that S matrix element describes braiding, which can in turn be seen as linking of worldlines.
Here we present an analogous interpretation of a 3+1D S matrix element. The algebra of membrane operators follows from their definition: Consider the membrane operator product G −1 F −1 H −1 F HG , where the expectation value is obtained in state |a, e, c with a, c arbitrary (note that such a state is in the ground state manifold, see after Eq. (11)). Here we label F = F so the quantum amplitude equals S matrix element v, w, µ|S|u, v, λ up to factor |G|! As in the previous subsection, the membrane operators Eq. (18), (24), (26), are interpreted as events of tunneling a flux-loop-antiflux-loop pair across a plane in the periodic system; for example, from Eq. (18), the F (z) u,λ describes the tunneling of loop in zy plane. The quantum amplitude in Eq. (36) can then be seen as the "movie" in Fig. 11.
Analogously to the 2+1D case, where two worldlines, defined by a particle tunneling operator and its inverse, were connected into a single worldline, we can consider that the events defined by F and F −1 form a single worldsheet, and so on for G and H. Having exactly three worldsheets, we calculate their TLN. Projecting the "movie" onto the three dimensional space slice at time t = −∞, we find eight triple intersection points of the projected worldsheets, Fig. 12. For simplicity of presentation, we offset the spatial position of inserted operator and its inverse, i.e., the membrane is moved slightly between the time of its appearance and disappearance. We checked that this offset does not influence the result. Notice that the orientation of membrane operator and its inverse are opposite, see Appendix B. A straightforward calculation from each triple point gives: a, c, e : T lk F HG = 1, b : T lk HF G = 1, f : T lk GHF = −1, e, g, h : T lk GF H = −1, which is exactly the same triple linking number obtained from S matrix calculation!

Braiding process of three flux loops
In previous subsections, we have interpreted the membrane operators as representing an instantaneous event of creating a loop-anti-loop pair and expanding the loop across a plane in the periodic system until the loops annihilate. However, this kind of worldsheet evolution can be smoothly deformed to represent a more physically clear  33)). For example, at point c, t2 shows that projection of H at this point comes at later times, so after F , while t1 shows that G comes from earlier times, so before both F and H, altogether contributing to T lkHF G. This process has same triple linking number as the one in Fig. 12. process. We therefore make a movie of three-flux-loop braiding process that gives exactly the same nontrivial triple linking number as the membrane process, as shown in Fig. 3. By projecting this braiding movie, we get Fig. 13, in which it is straightforward to measure the TLN: Triple point a gives Tlk GHF = −1, triple point b gives Tlk GF H = −1, triple point c gives Tlk HF G = 1 and triple point d gives Tlk F HG = 1. The obtained values of T lk IJK for the three-flux-loop braiding exactly match the membrane calculation result.

V. EXAMPLES
Here we will present the example of G = Z 2 × Z 2 cohomological gauge theories. Since H 4 (G, U (1)) = Z 2 × Z 2 , they can represent different topological orders. This will show how the loop statistics can distinguish different topological orders.
It is convenient to label group G elements a as (a 1 , a 2 ), where a i ∈ {0, 1}. Group multiplication rule a · b is defined as ( a 1 +b 1 , a 2 +b 2 ), where we introduce notation x ≡ x mod 2.
Since the cohomology group is H 4 (Z 2 × Z 2 , U (1)) ∼ = Z 2 × Z 2 , it can be parametrized by 4-cocycles with multiplication rule The explicit form of these 4-cocycles is 26 It is straightforward to check that these ω indeed satisfy the 4-cocycle condition.
One can now work out the induced 3-cocycle β a and 2cocycle γ a,b using their definitions in Eq.(7) and Eq.(8). For induced 3-cocycle, we get It follows that the 3-cocycle β a can be expressed as where P a ij is some integer matrix. According to Ref.27, then the induced 2-cocycle must be a coboundary Altogether, for inequivalent 4-cocycles we get the induced 2-cocycle as Now, we are ready to calculate statistics of loops and particles. We will focus on |G| · w, u, ν|S|v, w, µ =χ v,w In the second equality, we have defined where χ µ (χ ν ) is one-dimensional linear representation of Z 2 × Z 2 . One can easily check the above definition ofχ µ andχ ν is consistent, due to γ a,b being a 2-coboundary. Labeling µ = (µ 1 , µ 2 ) as Z 2 × Z 2 group element, First, let us consider the case w = (0, 0). In this case, only the χ λ factors are non-trivial in second line of Eq. (44), which is interpreted as contribution from Aharonov-Bohm phase of braiding particles around fluxloops. In this case, the phase factor equals e iπ( µ· u− ν· v) , which is independent of choice of cocycle. Namely, statistics between particles and loops cannot distinguish different phases.
Then, we turn to the general case. We get an additional phase factor s l beyond e iπ( µ· u− ν· v) , and the s l factor comes from ε in Eq. (44). In other words, it is present even when µ = ν = 0, i.e., χ representations are trivial, so there are no charged particles. Therefore, s l represents statistics of flux-loops. We list s l obtained from different 4-cocycles as follows • ω 00 : s l = 1.
• ω 01 : s l = e We can see that flux-loop braiding can indeed distinguish different topological orders in 3+1D, recalling here that the membrane operator expression is identified with a particular type of three-flux-loop braiding. In particular, according to previous section we can identify the fluxloops (blue, red, black) in Fig. 3 with fluxes (u, v, w) here, and there are no charges present. Now, we turn to T matrix elementχ u,v λ (u) = u, v, λ|T 31 |u, v, λ . In the same way as above, we get • ω 00 :χ u,v λ (u) = e iπ λ· u .
While the e iπ λ· u can be interpreted as Aharonov-Bohm phase of particles going around loop, the remaining part also encodes information about loop statistics. While we do not have a proof at this time, we believe that this phase is related to the ribbon nature of flux loop, or in other words to a thickness of the membrane.

VI. DISCUSSION AND CONCLUSIONS
One of our main results is the construction of MES states on the three-torus for the 3+1D cohomological gauge theory, which can be trivially generalized to arbitrary number of unit-cells. The S, T transformation matrices take a simple form in this basis.
We discussed that the S-matrix elements are directly related to the braiding of loop excitations. The T -matrix elements, which are diagonal in the MES basis, correspond to the generalization of topological spin for loop excitations. Here physically the loop excitations are generally expected to be ribbon excitations with two different loop-edges. We expect that the geometrical interpretation of the T -matrix elements is related to the braiding involving different loop-edges.
Although we use exactly solvable models and 3+1D topological quantum field theories to compute their S, T matrices, these 3+1D S, T matrices are in principle measurable quantities in practical model Hamiltonians. In particular, given a topologically ordered phase in 3+1D with its topologically degenerate ground sector on threetorus T 3 , one can firstly find a MES basis, similarly to the algorithms proposed in 2+1D. 13 For instance, for the S-matrix element between two MES |Ξ i and |Ξ j : S ij , one can perform the following thought numerical measurement. Because the topological properties do not depend of local geometry, we can assume that these ground states live on a cube with periodic boundary conditions. Then one can consider the state rotated by 120 • along the (111) direction of the cube: R 120 • |Ξ i . Because R 120 • |Ξ i and |Ξ j belong to the same topological phase, in the absence of symmetry there should exist a Hamiltonian path H(τ ) (τ ∈ [0, 1]) such that |Ξ j (|Ξ j ) are the ground state of H(0)(H(1)), and the ground state sectors of H(τ ) are adiabatically connected. One can then define a projection operatorP τ into the ground state sector of H(τ ) for any given τ . The many-body quantum amplitude related to the adiabatic time-evolution process of the S-transformation can be computed as Ξ j |P N −1/N · ... ·P 2/N ·P 1/N R 120 • |Ξ i as N → ∞. This computation is a realization of the topological quantum field theory time-evolution.
We expect that this quantum amplitude is related to the S-matrix elements s ij at most by an overall ambiguity U (1) phase e iθ , which is due to the nonuniversal local physics in the time-evolution, and a phase e iφi−iφj which is due to the gauge choice of |Ξ i ,|Ξ j . Even with these ambiguities, such measurements can still be used to extract useful information about the S, T matrices which potentially could fully determine them.
Recently, there has been a lot of progress in relating topologically ordered phases to symmetry protected topological (SPT) and symmetry enriched topological (SET) phases, for example by partially or completely ungauging the gauge group G, i.e., by transformations between global and local symmetries. 18,[28][29][30][31][32][33][34][35] We therefore expect that our work will be useful in characterization of SPT and SET phases too.
Our discussion has been limited to the case in which both group G and its projective representations are Abelian. These constraints are introduced here for simplicity rather than due to difficulty of principle. In the most general case, we expectχ to be the character of a projective representation, while modular transformations can easily be generalized, as shown in Ref. 24 for 2+1D case. In those cases, three-loop braiding may give a unitary matrix which corresponds to non-Abelian statistics of loops. In fact, we expect that the triple linking is relevant for any 3+1D topological order, beyond cohomological models. Namely, as long as the MESs in 3+1D topologically ordered state can be written using tunneling operators of loop-like excitations, which we believe should always be the case, our general geometrical approach based on space-time embedding shows that MES overlaps, which appear as S-matrix elements, will involve the triple linking of worldsheets defined by those tunneling operators.
Finally, let us consider a trivial but ubiquitous example of G = Z 2 . In this case, H 4 (G, U (1)) = Z 1 , so the cocycle can be set to identity map. The braiding phaseχ vw µ (u)/χ wu ν (v) reduces to a linear representation χ µ (u)/χ ν (v), where group elements u, v = 0, 1, and µ, ν = 0, 1 label the representations of Z 2 : The braiding phase therefore equals e iπ(µu−νv) . There is no contribution from flux-loop braiding, since the 1cocycle factors in Eq. (31) are trivial. In summary, the modular S transformation for common Z 2 gauge theory in 3+1D tells us that particles see a flux-loop as a π-flux, and the flux-loops themselves have trivial braiding. Using the MES basis and Eq. (29), (30), we directly obtain the S and T matrices of 3+1D Z 2 theory in their canonical form:  Diag(1, 1, 1, 1, 1, −1, 1, −1), where the MES basis |a, b, λ , with a, b, λ ∈ {0, 1}, is here naturally ordered according to binary numbers with digits abλ. These matrices are consistent with the S and T matrices derived for the same theory in Ref. 8. During the preparation of this manuscript we noticed a very recent paper that also considers aspects of flux-loop braiding in 3+1D. 36 This work was supported by the Alfred P. Sloan foundation and National Science Foundation under Grant No. DMR-1151440.
As t → −∞, r 2 , r 3 → 0, so the space manifold shrinks to a circle in XY plane: X 2 + Y 2 = r 1 [t = −∞] 2 . As t → ∞, r 1 , r 2 , r 3 → ∞, and the space manifold Σ asymptotically approaches the W -axis, which can be viewed as the S 1 (3) circle (γ circle) with r 3 = ∞. Using all time (−∞, ∞) with these asymptotic limits in the embedding is useful for consistently removing ambiguities in the embedding of events corresponding to membrane operators, as will soon become clear.
It is useful to think about this equation as a two-step process in time: from t = −∞ to t = 0 and from t = 0 to t = ∞. At t = 0 − , the system is in state |u, v, λ , obtained by insertion of membranes G u,λ at some time t < 0 into the trivial MES |e, e, 1 . On the other hand, the interval from t = 0 to t = ∞ can be interpreted as the conjugate of an appropriate history from t = −∞ to t = 0. This gives us the bra µ, v, w| by simply inserting membranes G and H into the trivial MES at t = ∞. Joining the two histories at t = 0, the product of bra and ket is obtained. So the quantum amplitude which equals the S matrix element is naturally interpreted as a sequence of space-time events.
Crucially, the coordinates x, y, z used for the membranes are in the original three-dimensional space, and we yet have to consistently identify them with the α, β, γ coordinates which were defined above in the embedding of space into R 4 . It turns out that F creates a sheet covering S 1 (1) ×S 1 (2) , G creates the sheet S 1 (3) ×S 1 (1) and membrane H covers S 1 (2) × S 1 (3) , as we show below. In other words, if we cut open the three-torus T 3 to a cube, then x axis corresponds to γ, y axis corresponds to β and z axis corresponds to α, remembering that F is membrane in the yz plane, G is membrane in the zx plane and H is membrane in xy plane.
To see this, we start from the trivial MES ket |e, e, 1 = 1 √ |G| c |e, e, c , where the three labels inside the ket correspond to the three directions x, y, z. Consider the corresponding limit t → −∞, in which the β, γ circles shrink to points, while the α circle remains at finite radius r 1 . Therefore the product of group elements along the α circle remains unconstrained in this limit. On the other hand, we conclude that the product of group elements along β or γ direction has to be identity e to be shrinkable to a point. More concretely, consider a consistent triangulation of entire R 4 , so that the nested Σ manifolds coming from different times are connected by the consistent triangulation. Then the β and γ circles at times t < 0 are bounding consistently triangulated discs, and as they shrink the zero-flux rule through the discs would force the product of group elements along the circles to identity e. Altogether, only the α coordinate can be identified as z, while β, γ correspond to x, y (not necessarily in that order).
We can perform a similar analysis for bra 1, e, e| and the corresponding limit t → ∞. The finite loop at t = ∞ is S 1 (3) , the γ circle. Then, it is easy to confirm that x corresponds to γ. So, the remaining axis y corresponds to S 1 (2) , i.e., to the β circle. Finally we discuss some conventions we have chosen in the calculations of TLN in both Section I and IV C. As described after Eq. (33), the TLN requires the calculation of normal vectors to worldsheets which are projected from four to three dimensions. The overall sign of T lk IJK therefore depends on precise definition of orientations of worldsheets, but these orientations are not inherent in the membrane operators and we need to choose them consistently. The F membrane operator defines a sheet in the yz plane, and we assign it the ordered pair (z, y). After a projection to three dimensions, the derivatives with respect to z, y define the sheet tangent vectors n 1 , n 2 , respectively, and the normal vector to sheet F is chosen as n 1 × n 2 . For F −1 , the normal is defined with opposite sign. In Section IV C, the total worldsheet composed of F and F −1 therefore has a consistent normal vector throughout. In the exact same fashion, we assign to G sheet the pair (x, z), and to H the (y, x) pair. Finally, we note that the Jacobian of the transformation in Eqs. (B1) has a negative determinant, so the embedding reverses the handedness of the coordinates in the three-dimensional projected space. Due to this, the TLN values calculated from the embedding are additionally multiplied by −1.